再入弹头三维非对称烧蚀外形模拟

周述光，国义军，贺立新，刘骁

1.中国空气动力研究与发展中心 空气动力学国家重点实验室，绵阳 621000 2.中国空气动力研究与发展中心 计算空气动力研究所，绵阳 621000

再入弹头三维非对称烧蚀外形模拟

周述光1,*，国义军2，贺立新2，刘骁2

1.中国空气动力研究与发展中心 空气动力学国家重点实验室，绵阳 621000 2.中国空气动力研究与发展中心 计算空气动力研究所，绵阳 621000

远程战略导弹再入大气层，热环境十分严酷，烧蚀量很大，外形呈现非对称性。烧蚀变形、质量引射和表面烧蚀花纹反作用于气动力、气动热及表面温度场。它们之间的相互作用十分复杂，影响落点精度。依据等效转换，将有迎角、侧滑角的绕流转换成总迎角绕流。对于不同子午面，将其变换到纵向坐标轴沿等效流动方向的气流坐标系，通过等效迎角化为零迎角绕流。采用Monte-Carlo统计方法，生成伪随机数，对弹头表面粗糙度随机抽样，表征不同子午线上转捩初始位置的随机分布。改进了零迎角烧蚀外形计算方法，拓展了其应用范围，并通过与能反映转捩点移动方向的松耦合计算策略相结合，采用数值求解烧蚀外形方程，成功模拟了小倾角再入过程中端头非对称烧蚀外形。

再入弹头；非对称体；热防护；气动热；烧蚀；碳/碳材料；松耦合模拟

传统远程导弹弹头以零迎角惯性再入，再入时间短，热流高，烧蚀严重，烧蚀外形基本上呈现对称性。近年来，新型远程战略导弹为提高机动性，采用持续有迎角再入弹道，再入时间长，尽管峰值热流比传统弹头小很多，但总热量更大，烧蚀外形为非对称。烧蚀外形/气动力/气动热/温度场/烧蚀等紧密耦合，相互影响，不确定因素多，预测精度差。为了提高飞行安全性能，通常加大防热系统的安全裕度，使得防热结构重量增加。为此，诸如碳/碳等烧蚀型防热材料的烧蚀外形变形规律得到了大量研究。利用电弧加热设备、高频等离子风洞、燃气流设备等，结合对飞行回收物的分析和理论研究，对弹头烧蚀外形有了基本认识。1970年，Thyson等预测了层流和湍流中理论外形的变化[1]，并比较了地面试验结果，意识到粗糙度的重要性。随后，Baker也预测了零入射角和零迎角下的烧蚀[2]，利用地面试验设备研究了有迎角下樟脑制成的球锥试样烧蚀特性，并推断是物体表面粗糙度引起了湍流热流的增加。转捩依据粗糙度的高低先后发生，热流也随之发生变化。不同流态下的烧蚀速率相差甚大，烧蚀外形与边界层的流动状态密切相关。例如，球锥外形的三向编织碳/碳端头，依据弹头再入飞行过程中边界层流动状态的变化，一般会出现三类6种烧蚀外形[2]：层流烧蚀外形(A、B)、转捩烧蚀外形(C、D)和湍流烧蚀外形(E、F)，如图1所示[3],RD、RE、RF分别为外形D、E、F的极径。

在再入弹头热防护计算研究的早期阶段，忽视端头外形的小不对称性，例如，20世纪70年代，Crowell研究了在弹头零迎角再入过程中，碳/碳材料的端头烧蚀外形变化规律[4]。Chin研究了端头的非对称烧蚀快速计算方法，发展了适用于轴对称外形的三维烧蚀计算程序[5]，方程基于球面、可变形、移动坐标系，计算表明石墨端头在10°迎角下烧蚀外形呈现显著的非对称性。之后，Dirling[6]根据端头回收体(Nosetip Recovery Vehicle, NRV，见图2)项目的飞行试验实物[7]，针对弹头非对称外形，采用材料表面粗糙度的随机分布来解释周向边界层转捩不一致性，并引入了Monte-Carlo分析法。Hall和Nowlan对控制方程采取有限差分-积分法和统计散布分析的近似技术[8]，分析了小不对称端头的气动特性。

图1 弹头典型烧蚀外形
Fig.1 Classical ablation shapes of missile nosetip

图2 整体石墨NRV再入飞行返回端头烧蚀外形
Fig.2Ablative shape of nosetip of recovered bulk graphite NRV flight test

国内从20世纪70年代也开始了再入弹头烧蚀防热的工程计算方法研究，在庄逢甘院士、张涵信院士等推动下，姜贵庆等开展了烧蚀防热机理研究及其工程计算研究[9-11]，先后建立了零迎角烧蚀计算方法、有迎角解耦计算方法等，并基于炮风洞开展了粗糙度对转捩、热交换等影响的试验研究。20世纪80年代以来，对非对称烧蚀问题进行了持续研究，其中又以在风洞模拟环境下测量气动力、热、烧蚀的研究为主[12-15]。常用的烧蚀防热快速算法大部分为对称烧蚀[16-17]，在非对称烧蚀外形方面的研究并不多见，袁湘江等[18]根据Dirling[19]介绍的方法，采取解耦计算方法，大致复现了零迎角状态NRV的小不对称烧蚀整体外形，模拟的沟槽不够明显，这可能与计算方法和姿态角选取有一定关系。由于端头外形对升阻比和稳定性影响很大，时至今日仍有报道，最近，陈自发等针对有迎角或滚转角的情况，计算得出碳基端头迎风面烧蚀量大于背风面的结论[20]。

在所有烧蚀外形计算方法中，解耦计算方法是最常见、最简单的，即首先采用初始外形沿弹道计算所有时刻的热环境，然后根据热环境再计算烧蚀外形，不考虑烧蚀外形变化对热环境的影响，可简单表述为“初始外形→热环境→烧蚀外形”。烧蚀量小的载人飞船返回舱、战术导弹等常采取解耦计算[4-5,21]。而烧蚀防热型的战略弹头，由于烧蚀量大，外形的重要性凸现出来，不能再采取解耦计算。烧蚀对气动力影响很大，体现在3个方面：烧蚀外形变化，会明显改变飞行器的气动性能；烧蚀过程中的质量引射对边界层状态的影响也会反映到气动性能上；烧蚀形成的小不对称与质量、惯性的分布不对称性，对运动产生较大影响，反过来影响气动性能。文献[10]指出气动力、气动热、弹道运动的耦合计算内容包括：飞行器运动特性、外形计算、气动特性、烧蚀计算。弹头再入过程中，如果考虑它们之间的作用，即为所谓(紧)耦合计算。由于涉及面广，难度很大，目前还未见到这方面的研究报道。

总的来说，国内对带迎角、侧滑角的飞行状态下的烧蚀模拟技术欠缺；如何合理地模拟弹头表面粗糙度影响值得研究；为了提高模拟精准度，还需要考虑每一时刻的端头外形对气动热环境的作用。

为此，本文烧蚀外形计算中，不考虑烧蚀外形变化时气动特性的改变对运动特性(弹道参数)的影响，沿着弹道，根据上一时刻的烧蚀外形，计算当前时刻热环境参数(热流、恢复焓、温度、边界层外缘速度、压力等)，然后根据烧蚀方法，刻画出这一时刻的外形，作为下一时刻热环境的输入依据，即所谓烧蚀外形松耦合计算，可简单表示为“热环境↔烧蚀外形”。本文在零迎角、松耦合计算方法验证算例的基础上[22]，研究了由于粗糙度沿周向非均匀分布对烧蚀外形产生的影响，建立了带迎角、侧滑角的松耦合计算方法，实现了非对称烧蚀外形预测。

1 烧蚀外形模拟技术

1.1 坐标变换

气动力和气动热参数计算均在坐标原点位于弹头头部的坐标系中进行，而烧蚀导致头部物面向后退缩，坐标原点随之后移，每一时刻需要重新确定物面坐标值。此外，头部顶点还径向偏离原始物体顶点。对于此两类情况，一个较好的办法是采用瞬时坐标系，通过坐标变换，建立新旧坐标系之间的联系(图3)。设Oxyz为初始时刻的固连直角坐标系，坐标原点O位于弹头顶点。定义O*x*y*z*为当前计算时刻的瞬时直角坐标系，坐标原点O*始终位于弹头烧蚀外形的前部顶点处，各坐标轴的方向与初始坐标系一致。子午面分布不变。瞬时坐标系中的某一物面点P，在原始坐标系中为(xP，yP，zP)，子午角为φ0。瞬时坐标系O*x*y*z*与初始坐标系Oxyz之间的关系为

(1)

式中：δx、δy、δz为两个坐标系的坐标原点之间的距离Δ分别在x、y、z轴上的投影。

图3 初始坐标和瞬时坐标示意图
Fig.3Schematic of original and instantaneous coordinates

旋成体非常适合采用球坐标系(R,θ,φ)来描述，设初始时刻的球坐标原点N′位于x轴上，某一时刻移原点至N，R为原点到物面的极径；θ和φ分别为球心角和子午角。则初始直角坐标系与球坐标系之间的关系为

(2)

由烧蚀计算可确定各子午面新的R值，并由式(2)计算物面上各点的直角坐标x、y、z，然后寻找整个物面上的最小x坐标点，从而确定δx、δy、δz。最后由式(1)确定新坐标(x*，y*，z*)，开始下一时刻的流场参数计算。

1.2 烧蚀外形数值求解

烧蚀外形严格说来是一个三维问题，但对于实际飞行的弹头而言，烧蚀外形沿周向的变化率远远小于纵向变化率。因此，通常当作二维问题处理，仅在每个子午面内解二维烧蚀外形方程。球坐标系下，描述烧蚀外形变化历程的控制方程为[4，23]

(3)

(4)

记式(4)右边为函数f，它是Rθ、R、θ、φ、来流马赫数Ma∞的函数，即

则式(4)的空间半离散方程为

对应的修正方程为

(Rt)j=f(Rθj+ε,Rj,θj,Ma∞,φ)

(5)

式中：ε为一小量。式(5)展开得

(6)

2 姿态角等效处理

(7)

由此，可将有迎角和侧滑角的绕流问题转化为单纯总迎角绕流问题。

对所研究的子午面，将其变换到纵向坐标轴沿等效流动方向(与x坐标轴成迎角αEB)的气流

图4 等价变化坐标示意图
Fig.4Schematic of equivalence transformation of coordinate

坐标系轴线，可以进一步化为零迎角绕流问题。运用等效体方法，可得到相对于等效零迎角旋成体绕流的物面等价倾角，即

(8)

式中：n为物面外法线方向的单位矢量。在球坐标系中，定义物面方程F为

Fr,φ,x=r-fx,φ=0

(9)

式中：r为物面上的点到轴线(x轴)的距离。则物面外法线单位向量为

(10)

式中：ex、er、eφ为球坐标系中的3个单位向量。由式(11)定义物面角(切向和周向)，即

(11)

式中：-π/2≤θb≤π/2；-π/2≤δφ≤π/2。由式(11)可导出

(12)

将式(11)和式(12)代入式(10)，可得

n=-exsinθb+ercosθbcosδφ-eφcosθbsinδφ

可记为

n=-isinθb+jcosθbcosφ-δφ+

kcosθbsinφ-δφ

(13)

式中：i、j、k为空间直角坐标系中的3个单位向量。

速度矢量在瞬时坐标系中的表达式为

V∞=V∞icosβcosα+jcosβsinα+ksinβ

(14)

将式(13)、式(14)代入式(8)，可得

θbEQ=arcsin(cosαcosβsinθb-sinαcosβcosθb·

cosφ-δφ-sinβcosθbsinφ-δφ)

(15)

因为φ=0°为背风面，φ=180°为迎风面，所以式(15)中的α必须反号，即-α。定义等价迎角为

αEB=θb-θbEQ

(16)

至此，建立了任意子午面的等效迎角与实际飞行迎角和侧滑角的关系。通过以上变换，将有迎角和侧滑角的绕流问题化为等价零迎角绕流问题，这样就可以运用以往成熟的零迎角情况下的计算方法，求得表面压力分布、激波形状以及表面热流分布等流场参数。

3 热环境和烧蚀计算方法

3.1 边界层外缘参数计算

球锥外形弹头的外部流场考虑端头烧蚀形状的变化，表面压力分布按照端头部位和锥身部位分别计算。为了计算激波层中熵梯度效应对热流的影响，首先计算每个时刻烧蚀外形的头部激波形状，获得实际激波角和激波脱体距离，再计算激波后参数，根据这些参数对边界层外缘焓进行修正，利用均值法，获得边界层外缘平均速度。

3.2 表面热流计算及粗糙壁热增量修正

对于不可穿透物形表面的层流、湍流热流，采用有效长度法计算；考虑熵梯度效应，应用平均值法；粗糙壁热增量即烧蚀表面粗糙度对热流的影响，采用粗糙壁热流放大因子进行修正，即

qj=Fjqjh0,u0,ρ0,μ0j=l,t

(17)

式中：Fj为粗糙壁热增因子；q为热流；h为焓值；u为边界层外缘速度；ρ为密度；μ为黏性系数。上标“0”表示外部主流区；下标l、t表示层流和湍流状态。

转捩区热流用间歇因子Γ加权平均，即

qtr=(1-Γ)ql+Γqt

(18)

转捩区间歇因子Γ的值在0(层流)～1(湍流)之间变化，在不同雷诺数范围内有不同的表达式，即

(19)

式中：Re1为转捩开始雷诺数；Re2为平均雷诺数；Re3为结束雷诺数；Reθ为当地雷诺数。

3.3 烧蚀速率计算方法

端头的物面法向线烧蚀速率由材料的有效焓计算，将物面视为准稳定的加热面，对于碳/碳复合材料，有

(20)

4 粗糙度随机分布模拟

弹体表面粗糙度在烧蚀外形计算中扮演着重要的角色，它不仅影响边界层转捩，也影响表面热流，决定了烧蚀速率的快慢，使得端头外形出现深浅不一的沟槽，逐渐呈现小不对称性。通常情况下，表面粗糙度分布是随机的，另外还有一些小不对称因素也是随机分布的，因此，烧蚀外形也具有随机的小不对称性。

kj=kL+σKXj

(21)

式中：kL为测量试样得到的名义粗糙度高度；σK为粗糙度测量值的名义标准偏差；Xj为伪随机数列，由标准的伪随机数列程序产生，周期足够长，以保障其随机性。在任意的一组随机抽样值中，选出其最大值Kπ1，最小值Kπ2，作为对应子午线φ=φ*和子午线φ=π+φ*的粗糙度值，通过这两个值，确定中间子午线上的Kπ值，即

(22)

式中：φ*为最大粗糙度值所在的子午角。

5 计算结果分析

5.1 结果验证

弹头外形参考NRV的设计[6]，为球-锥旋成体构型，端头半径为31.75 mm，底部半径为151.13 mm，总长为1 168.4 mm，半锥角为6°。根据Dirling的建议[6]，粗糙度因子取为15.5 μm，均方差为8.1 μm。综合文献给出的部分弹道参数[6-7]，见表1，并将弹头姿态角沿弹道全程设为0。验证计算结果见图5(纵横比非1∶1，下面所有的外形、轮廓图均如此)，给出了端头外形。可知，模拟出端头上多道深浅不一的烧蚀沟槽，与回收物烧蚀外形类似。端头驻点后退量约为7.6 mm，落在Otey和English预测的7.4～8.9 mm范围内[6]，而对飞行试验回收的端头进行测量，获得的后退量约为8.5 mm。考虑到弹道参数不全，粗糙度影响等原因，可以认为本次验证计算结果合理，程序可靠。

表1 弹道参数Table 1 Parameters of trajectory

图5 验证计算端头烧蚀外形
Fig.5Verification simulation of ablation shape of nosetip

5.2 姿态角与端头烧蚀外形非对称性

不同迎角和侧滑角会引起迎风面和背风面的热流不同，烧蚀速率也就出现差异，导致烧蚀外形的不对称性。这里将粗糙度方差设置为0，即消除粗糙度不均带来的影响。图6给出了同一飞行条件下有/无姿态角的端头烧蚀外形，dH表示当前时刻外形与原始外形之间变化量；图7给出了其中3个状态下纵向(z=0)剖面轮廓图，横坐标为轴向坐标x，纵坐标为剖面上轮廓的周向直径r。端头顶点后退量均约为7.4～7.5 mm，轮廓线的上半部分为背风面，下半部分为迎风面，不对称性随总迎角的增大而增强。图8给出了Λ=15°的姿态角弹道、不同高度下的端头后退量。可知，有姿态角时，迎风面烧蚀量大于背风面，呈现非对称烧蚀外形；此外，在30 km以上，烧蚀量较小，25 km以下，烧蚀量增大。当迎角和侧滑角相同时，来流方向在子午角45°平面里，如图9所示(α=β=10.6°，Λ=15.0°)，这意味φ=45°的子午面轮廓不对称性最显著，φ=135°的子午面轮廓为对称性，纵向对称剖面(φ=90°)和水平剖面(φ=0°)的轮廓一致，均为程度较轻的小对称性。

图6 端头烧蚀小不对称外形
Fig.6 Small asymmetric ablation shape of nosetip

图7 端头烧蚀小不对称纵向剖面轮廓图
Fig.7Outlines of small asymmetric ablation shape of nosetip at z=0 section

5.3 粗糙度与端头烧蚀外形非对称性

在粗糙度研究中，模型姿态角取为0°。由于热流会明显影响烧蚀外形，为了便于观察，这里选用一条环境严酷的弹道进行模拟。由图10可知，烧蚀的不对称性是随机分布的，其中RL为层流粗糙度因子，RT为湍流粗糙度因子，MVL为层流粗糙度因子分布均方差，MVT为湍流粗糙度因子分布均方差。无论是层流粗糙度还是湍流粗糙度，随着粗糙度的降低，烧蚀外形的不对称性程度减轻。与湍流状态相比较，层流状态下的粗糙度因子对热流作用更显著。这可能源自层流流动对扰动更敏感，更容易向转捩发展，热流变化相对大些，烧蚀速率差异较大；湍流抗干扰能力较强，流动稳定性更好，热流变化也相对缓和，烧蚀速率较为均衡。

5.4 耦合-解耦计算结果分析

在耦合-解耦计算策略影响研究中，模型姿态角取为0°，粗糙度均方差设置为0，烧蚀外形将具有对称性。图11给出了沿同一条弹道计算， 耦合计算与非耦合计算的结果。由图可知，耦合计算结果显示为湍流烧蚀外形，顶点后退量约为28.4 mm；而解耦计算结果显示为转捩烧蚀外形， 顶点后退量约为27.4 mm。耦合情况下，流场参数根据烧蚀外形计算，开始时转捩点位于声速点附近，转捩后的热流增加，烧蚀量变大，随着外形变化，转捩点向端头顶点移动，致使端头越来越尖，趋于双锥外形。而解耦情况下，采用原始外形计算热流，对于保持不变的外形，转捩位置基本固定，转捩点后热流高，烧蚀量大；转捩点前为层流，烧蚀量小。与之前的非耦合计算结果比较[22]，通过优化解耦计算方法，消除了以前的非物理波动，模拟结果不再为凹陷外形。而耦合计算考虑了“热环境↔烧蚀外形”的双向作用，模拟更精细，结果更符合真实情况。此外，差分离散依据情况自动选择隐式或显示格式，与无波动、无自由参数、耗散的(Non-oscillatory containing No free parameters and Dissipative，NND)格式相比较[22]，对烧蚀外形仿真的模拟程度二者相当，但前者效率更高。

图8 端头不同子午剖面的轮廓图
Fig.8 Outline of nosetip at different meridian sections

图9 迎角与侧滑角相等时气流所在子午面
Fig.9Meridian section of flow direction with angle of attack being equal to sideslip angle

图10 层流和湍流粗糙度条件下的随机烧蚀外形
Fig.10 Asymmetric ablation shapes with different roughs in laminar flow and turbulent flow

图11 耦合-解耦计算烧蚀外形
Fig.11Ablation shapes with coupled/decoupled simulation

6 结 论

本文通过坐标变换、外形方程数值求解等外形重构技术，建立了三维烧蚀外形计算方法。通过定义总迎角和等效迎角的处理方法，将零迎角烧蚀方法拓展到非零迎角、侧滑角飞行条件下，应用于不同姿态角下高超声速飞行器的热环境模拟问题，适用于旋成体弹头小倾角再入烧蚀外形计算。

将Monte-Carlo随机抽样技术应用于弹头物面的层流和湍流粗糙度，结合粗糙壁边界层转捩及热流增益效应研究，模拟了弹头烧蚀外形的非对称沟槽。降低粗糙度可以减轻小不对称程度，因此，通过提高弹头表面致密性和表面光洁度将有助于烧蚀外形的对称性和连续性。

将气动热、烧蚀和三维外形变化进行双向耦合计算，采用此策略避免了解耦计算中转捩点位置固定，而是随着烧蚀外形的不断变化，逐渐向端头顶点移动，符合实际情况，因此提高了对弹头再入烧蚀过程的模拟程度。

在现有零迎角烧蚀外形计算方法的基础上，考虑飞行姿态角和表面粗糙度两个重要因素，结合松耦合计算，建立了一套新的烧蚀计算方法，成功地对烧蚀非对称外形进行了预示，这将为精细的气动特性计算提供重要的外形参数。

[1] THYSON N, NEURINGER J, PALLONE A. Nose tipe change predictions during atmospheric re-entry[C]∥AIAA 5th Thermophysics Conference. Reston, VA: AIAA, 1970.

[2] BAKER R L. Low temperature ablator nosetip shape change at angle of attack[C]∥AIAA 10th Aerospace Sciences Meeting. Reston, VA: AIAA, 1972.

[3] 李素循. 烧蚀形状的地面实验模拟[C]∥航空航天工业部科学技术司论文集, 1990：379-384.

LI S X. Wind tunnel tests simulate ablation shapes[C]. Proceedings of Department of Science and Technology，Ministry of Aerospace Industry, 1990：379-384 (in Chinese).

[4] CROWELL P G. Finite difference schemes for the solution of the nosetip shape change equations: SAMSO-TR-76-8 [R]. Philadelphia: General Electric Company, 1976.

[5] CHIN J H. Shape change and conduction for nosetips at angle of attack[C]∥7th Fluid and Plasma Dynamics Conference. Reston, VA: AIAA, 1974.

[6] DIRLING R B. Asymmetric nosetip shape change during atmospheric entry[C]∥AIAA 12th Thermophysics Conference. Reston, VA: AIAA, 1977.

[7] OTEY G R, ENGLISH E A. High-βre-entry vehicle recover[J]. Spacecraft，1977, 14(5): 290-293.

[8] HALL D W, NOWLAN D T. Aerodynamics of ballistic re-entry vehicles with asymmetric nosetips[J]. Spacecraft, 1978, 15(1): 55-61.

[9] 姜贵庆, 刘连元. 高速气流传热与烧蚀热防护[M]. 北京: 国防工业出版社, 2003.

JIANG G Q, LIU L Y. Heat transfor of hypersonic gas and ablation thermal protection[M]. Beijing: National Defense Industry Press, 2003(in Chinese).

[10] 黄志澄. 航天空气动力学[M]. 北京:中国宇航出版社, 1994.

HUAGN Z C. Aerospace aerodynamics[M]. Beijing: China Astronautic Publishing House, 1994 (in Chinese).

[11] 薛成位, 王喜军, 张文杰, 等. 战略弹头气动攻关论文集[M]. 绵阳: 中国空气动力研究与发展中心, 1990.

XUE C W, WANG X J, ZHANG W J, et al. Proceedings of strategic missile's aerodynamics tackled[M]. Mianyang: China Aerodynamics Research and Development Center, 1990 (in Chinese).

[12] 蔡金狮. 再入体滚转共振与小不对称气动力[J]. 空气动力学学报, 1985(3): 72-81.

CAI J S. Roll resonance and small asymmetrical aerodynamics of reentry body[J]. Acta Aerodynamica Sinica, 1985(3): 72-81(in Chinese).

[13] 黄兴中, 蔡慧娟. 小不对称体滚动力矩风洞实验研究[J]. 空气动力学学报, 1985(2): 69-73.

HUANG X Z, CAI H J. Roll moment research of small asymmetrical body with wind tunnel test[J]. Acta Aerodynamica Sinica, 1985(2): 69-73(in Chinese).

[14] 白葵, 冯明溪, 付光明. 小不对称再入体滚转气动力测量技术[J]. 流体力学实验与测量, 2002, 16(3): 63-67.

BAI K, FENG M X, FU G M. Experimental technique for rolling aerodynamic of slight asymmetric re-entry body[J]. Experiments and Measurements in Fluid Mechanics, 2002, 16(3): 63-67 (in Chinese).

[15] 赵俊波, 付增良, 梁斌, 等. 再入弹头小不对称俯仰气动特性测量技术研究[J]. 实验流体力学, 2015, 29(5): 55-59.

ZHAO J B, FU Z L, LIANG B, et al. Research on the wind tunnel test techniques for micro-pitching-aerodynamics of re-entry body[J]. Journal of Experiments in Fluid Mechanics, 2015, 29(5): 55-59(in Chinese).

[16] 张毅, 苗育红, 周江华. 烧/侵蚀耦合效应及其对再入体落点的影响[J]. 飞行力学, 1997, 15(2): 91-96.

ZHANG Y, MIAO Y H, ZHOU J H. The influence of the coupled ablation/erosion effect on reentry-body impact-point[J]. Flight Dynamics, 1997, 15(2): 91-96(in Chinese).

[17] 陆海波, 陈伟芳, 袁雷, 等. 再入体碳基材料烧蚀特性分析与工程计算[J]. 弹道学报, 2008, 2(3): 107-110.

LU H B, CHEN W F, YUAN L, et al. Ablation analysis and engineering calculation on carbon based material of reentry vehicles[J]. Journal of Ballistics, 2008, 2(3): 107-110(in Chinese).

[18] 袁湘江, 杨茂昭, 潘梅林. 小不对称烧蚀外形的典型化处理方法[J]. 空气动力学学报, 1997, 15(3): 386-392.

YUAN X J, YANG M Z, PAN M L. Typical processing method of slight asymmetric ablative shape[J]. Acta Aerodynamica Sinica, 1997, 15(3): 386-392(in Chinese).

[19] DIRLING R B. Performance technology program(PTP-SII) Vol. VI. A statistical model of nosetip shape change for reentry vehicles: BMO-TR-80-53 [R]. Irvine: Science Application, Inc., 1980.

[20] 陈自发, 张晓晨, 王振峰, 等. 高超声速飞行器碳基头锥烧蚀外形计算[J]. 航空学报, 2016, 37(S1): S38-S45.

CHEN Z F, ZHANG X C, WANG Z F, et al. Hypersonic aircraft’s carbon-based nose ablation shape calculation[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(S1): S38-S45(in Chinese).

[21] 黄振中. 烧蚀端头的瞬变外形及内部温度分布[J]. 空气动力学学报, 1981(1): 53-65.

HUANG Z Z. Transient changing shape and temperature distribution of ablative nosetip[J]. Acta Aerodynamica Sinica, 1981(1):53-65(in Chinese).

[22] 国义军, 童福林, 桂业伟. 烧蚀外形方程差分计算方法研究(II:耦合计算)[J].空气动力学学报, 2010, 28(4):441-445.

GUO Y J, TONG F L, GUI Y W. Finite difference schemes for solution of the nosetip shape change equation(Part II, coupling calculation)[J]. Acta Aerodynamica Sinica, 2010, 28(4):441-445 (in Chinese).

[23] 张志成. 高超声速气动热和热防护[M]. 北京: 国防工业出版社, 2003.

ZHANG Z C. Aerothermal heating and thermal protection of spacecraft at hypersonic conditions[M]. Beijing: National Defend Industry Press, 2003(in Chinese).

[24] HALL D W. The three-dimensional shock and pressure (3DSAP) approximate flow field technique: SAMSO-TR-77-C145 [R]. Philadelphia: General Electric Company, 1977.

Simulationof3Dasymmetricablationshapeofreentrymissile

ZHOUShuguang1,*,GUOYijun2,HELixin2,LIUXiao2

1.StateKeyLaboratoryofAerodynamics,ChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China2.ComputationalAerodynamicsInstitute,ChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China

Strategicmissilesareheatedbyhypersonicgaswhiletheyarereenterintotheearth’satmosphereatsmallattitudeangles.Thethermalenvironmentissoharshthatthelossofsurfacematerialfromspacecraftthroughevaporationormeltingcausedbyfrictionwiththeatmosphereisaconsiderableamountofablation,therebyformingavisibleasymmetricshape.Shapechange,matterinjecttoboundarylayerandanumberofdifferenttypesofpatternsbyablationwillaffecttheaerodynamics,aerothermodynamicsandtemperaturedistributioninthematters.Thermalenvironment,ablation,andtheshapearethecauseandeffectforeachother,andthereforemaketheimpactscattertospreadout.Accordingtoequivalenttransformanalysis,theflowfieldrelativetothevehicleatsomeanglesofattackandsideslipcanbeequaltothatatatotalangleofincidence.Theverticalaxisatameridianplaneistransformedtotheflowdirection,whereforetheflowfliedisatthezeroangleofattack.RandomsamplingofsurfaceroughsisconductedwithMonte-Carlostatisticalmeasuretosimulatetherandomdistributionofthetransitionregionsatdifferentcircumferentialpositions.Theroughsizesarecreatedbytherandomizer.Thecalculationprocedurehasbeenimprovedwiththesetechniques,andthelooselycoupledsimulationthatcanexhibitthemovementofthetransitionpointtowardsthevertexofthenosetipisalsoemployed.Thenumericalcalculationmethodisalsousedtosolethenonlinearpartialdifferentialequationfortheablativeshape.Theirregularshapeofthenosetipofballisticmissilesatsmallincidenceangleinreentryintotheearth’satmospherecanbesuccessfullysimulatedwiththemethodsappliedproposed.

reentrymissile;asymmetricbody;thermalprotection;aerothermal;ablation;carbon/carbonmaterial;looselycoupledsimulation

2017-05-08；

2017-05-31；

2017-06-20；Publishedonline2017-06-291323

URL：http://hkxb.buaa.edu.cn/CN/html/20171211.html

NationalKeyResearchandDevelopmentProgramofChina(2016YFA0401200)

.E-mailzhoushuguang111@163.com

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2017.121397

2017-05-08；退修日期2017-05-31；录用日期2017-06-20；网络出版时间2017-06-291323

http://hkxb.buaa.edu.cn/CN/html/20171211.html

国家重点研发计划(2016YFA0401200)

.E-mailzhoushuguang111@163.com

周述光，国义军，贺立新，等．再入弹头三维非对称烧蚀外形模拟J．航空学报,2017,38(12):121397.ZHOUSG,GUOYJ,HELX,etal.Simulationof3DasymmetricablationshapeofreentrymissileJ.ActaAeronauticaetAstronauticaSinica,2017,38(12):121397.

V11

A

1000-6893(2017)12-121397-12

李明敏)