高超声速弹丸气烧动蚀数值模拟研究

陈 蕊，董 刚

(南京理工大学 瞬态物理国家重点实验室， 南京 210094)

【装备理论与装备技术】

高超声速弹丸气烧动蚀数值模拟研究

陈 蕊，董 刚

(南京理工大学 瞬态物理国家重点实验室， 南京 210094)

弹丸飞行过程中会受到来流的气动加热，对于高超声速弹丸，气动加热会使弹丸产生烧蚀，改变弹丸形状。针对该问题，数值模拟了高超声速弹丸流场与结构传热和弹丸的烧蚀变形过程，对弹丸烧蚀过程特性进行研究。模拟结果表明，弹丸头部驻点位置烧蚀情况最为严重，弹丸烧蚀程度随飞行马赫数增加而增加；弹丸头部半径对烧蚀影响不明显，减小头部半锥角会减少烧蚀；钨的抗烧蚀特性优于钢材。研究结果对高超声速弹丸热防护设计有指导意义。

高超声速；气动烧蚀；弹丸；数值模拟

超高速发射系统武器是新一代的动能武器，这些新式动能武器以高超声速飞行，依靠自身的高动能产生破坏力，与常规的化学能武器相比，其穿甲、防空、反导等能力有巨大改观，对未来战争将产生重大的影响[1-2]。

新概念火炮发射弹丸炮口初速通常都要大于5马赫数，达到了高超声速的水平，而且这类弹丸多为水平发射，在飞行过程中高度的变化不明显，整个高超声速飞行过程是在大气最为稠密的底层进行。与常规炮弹相比，由于弹丸速度显著增加，使得弹丸与其周围气体质点的摩擦加剧产生很大的热量，在飞行初期受到剧烈的气动加热[3-4]。气动热可引起弹丸头部烧蚀，弹丸外形的一系列变化还将增加弹丸的飞行阻力，产生非平衡力矩，导致弹丸飞行失稳和外弹道轨迹变化。因此，研究弹丸表面气动热和烧蚀所引起的外形变化对高超声速弹丸的设计具有重要参考意义。

T Chiang等[5]采用熔点烧蚀模型，通过有限元方法数值模拟了高超声速弹丸头部的烧蚀情况，研究指出熔点温度的选择不仅要根据弹丸材料，还要考虑弹丸的飞行环境对于材料的氧化程度。国内对于常规兵器气动加热引起的烧蚀现象研究还很少，史金光[6]对高超声速弹翼前缘表面气动加热及烧蚀现象进行了研究，模拟了多种条件下弹翼前缘的瞬时烧蚀情况，但没有考虑弹丸的烧蚀情况。杨晨[7]采用软件Fluent对弹丸气动烧蚀及内部传热进行了研究，仅对不同马赫数条件下弹丸烧蚀情况进行分析，且其中弹丸烧蚀速度采用经验常数。因此，对超高速飞行无尾翼回旋体弹丸烧蚀过程的特性需要进一步研究。本文数值研究了高超声速飞行条件下弹丸的烧蚀特性，分析了弹丸外形、飞行马赫数以及弹丸材料等多种因素对烧蚀特性的影响规律。

1 数理方程与计算方法

1.1 基本假设与控制方程

一般情况下，弹丸的质量分布和结构轴对称，且材料的热物性参数不变。弹丸飞行时，在边界层流动中，当化学反应速度比流体速度相对快得多时，化学反应能够将流体粒子很快地调整到局部热力学平衡值。基于这些情况，为简化高超声速弹丸飞行时气动加热和烧蚀的数学模型，本文采用了如下假设：(1)弹丸烧蚀过程中的氧化反应很快达到平衡状态；(2)弹丸飞行过程中攻角为零；(3)弹丸结构和质量分布轴对称；(4)材料热物性参数不变，为各向同性，无内热源。根据以上假设，本文采用的弹丸导热控制方程为

(1)

其中：t为时间(s)；λ为材料导热系数(W/(m·K)；ρ为材料密度(kg/m3)；cp为材料比热容(J/(kg·K)；λ/ρcp=a为热扩散率(m2/s)。

随着飞行时间增加，弹丸发生烧蚀，弹丸的形状不断变化。对于弹丸边界不断变化的问题，使用在运动表面上一个坐标保持不变的坐标系是很方便的。本文采用修正的Landau变换的坐标变换方法[6]，将运动边界问题变换成为固定边界问题，其坐标变换关系式为

(2)

式中：x0(y)是弹丸y方向坐标处的初始弹丸表面x坐标，s(y,t)是t时刻弹丸y高度处表面从初始位置沿x轴后退的距离，Ls为弹丸的长度。因此，φ(y,x,t)=0总是对应于弹丸瞬时烧蚀表面，φ(y,x,t)=1总是对应所研究弹丸底部表面。变换坐标系示意图如图1。

图1 变换坐标系示意图

采用修正的Landau变换对弹丸导热微分方程进行坐标变换后，式(1)可转换为如下的二维轴对称控制方程：

弹丸轴对称截面(y≠0)：

(3a)

弹丸轴线处(y=0)：

(3b)

对于弹丸烧蚀表面后退速率的求解，可以采用如下公式[6]

(4)

式中：β是弹丸表面上某点处的切线与x轴的夹角；F是单位质量材料的熔化热(J/kg)；q为表面热流密度(W/m2)。对式(4)进行积分可以得到弹丸表面的后退深度。

对于控制方程式(3)和式(4)，本文采用显式有限差分法进行离散，其中，等式左边的时间导数项采用1阶向前差分格式，温度对空间坐标的1阶偏导数采用1阶迎风差分格式，温度对空间坐标的2阶偏导数采用2阶中心差分格式，温度对空间坐标的混合偏导数采用1阶迎风差分格式。表面后退速率的偏导数采用1阶向后差分格式。

1.2 边界条件与初始条件

本文选取典型的球头-圆锥-圆柱体结构的弹丸飞行过程中的烧蚀现象进行数值模拟。由于烧蚀主要发生在头部区域，因此计算区域取为球锥型弹头及部分圆柱体，其中计算区域示意图如图2。

图2 计算区域示意图

图中Rn为头部半径，R为弹丸半径，θc为头部半锥角，Ln为弹丸头部长度。由于采用的是二维轴对称的控制方程，因此计算区域为图2中的阴影部分。

热流密度ql的求解基于化学反应很快达到平衡状态的假设条件，利用不可压缩的Blasius关系式，并通过Eckert参考物性法考虑可压缩效应进行层流热交换计算[8-9]，得到弹丸表面热流密度计算公式为

(5a)

驻点处的热流密度公式如下

(5b)

其中上标“*”代表参考值，下标“e”表示边界层外缘处的值。Reθ为边界层外缘雷诺数，μ为动力黏性系数，ue为边界层外缘速度，h为焓，下标aw和w分别表示绝热面和壁面。Pr为普朗特数，due/dx为驻点速度梯度。

对于本文所研究的气动烧蚀问题，其边界条件主要包括弹丸前迎风面气动加热热流密度条件，并且假设弹丸后表面绝热。初始条件假设弹丸内部处处均匀等温。

对于弹丸瞬时烧蚀表面，热流密度对温度的影响存在关系式[6]

(6)

根据式(6)，便可以获得温度的边界条件。初始条件假设弹丸内部处处均匀等温。

1.3 弹丸的烧蚀模型

在高温气流作用下，表面材料消失变形的现象(包括材料表面的熔化、燃烧、升华、氧化等)统称为烧蚀。在本文中烧蚀模型采用熔化模型[10-11]，即将弹丸温度达到材料熔点的部分去除模拟弹丸飞行过程中的烧蚀过程。考虑到当弹丸在飞行过程中，高温下弹丸材料会和空气中氧气等发生化学反应生成金属氧化物，而金属氧化物熔点一般较纯金属低，则需采用较低的金属氧化物熔点温度作为烧蚀判据。在每个时间步中通过对式(4)的积分得到弹丸表面烧蚀后退距离，弹丸表面后退并重新进行坐标变换，使得计算坐标φ(y,x,t)=0总是对应于弹丸瞬时烧蚀表面，φ(y,x,t)=1总是对应所研究弹丸底部表面。

2 计算结果与分析

2.1 计算算例及条件

本文考察了不同弹丸头部外形(头部半径Rn和半锥角θc)，弹丸材料和弹丸飞行马赫数对其烧蚀行为的影响，计算的算例及其参数见表1所示。

表1 计算参数

在所有计算中，飞行条件为海平面标准气象条件，弹丸半径R为10 mm，模拟飞行时间为2 s。算例1为参考算例，针对该算例，本文开展了网格无关性测试。计算结果表明，当网格数量为100×200时(x×y)，增加网格数量对计算的烧蚀外形没有显著影响，因此，所有算例均选用该网格数量进行计算。

图3给出了算例1(参考算例)在不同时刻的瞬时烧蚀外形及其温度分布。对于算例1的钨材料弹丸，本文的计算中假设钨的氧化反应过程很快完成(参见1.1节的基本假设)，因此以氧化物WO3的熔点(1746 K)作为烧蚀的标准。从图3可以看出，弹丸飞行初期(t=0.5 s)，头部迅速升温达到熔点发生烧蚀；随后头部继续烧蚀后退，高温区域逐渐沿锥面蔓延(t=1.0 s)；随着飞行时间的增加，弹丸内部温度迅速增加(t=1.5 s)；到t=2.0 s时，弹丸头部驻点烧蚀后退深度约为3.6 mm，约占弹丸头部初始长度的12%，说明弹丸头部已发生了明显的烧蚀现象。从弹丸头部烧蚀形状来看，烧蚀初期弹丸头部主要呈现平头形状(t=0.5 s和1.0 s)，后期随着锥面温度的增加，弹头形状又逐渐变成烧蚀后的球面形状。

图3 算例1的弹丸烧蚀外形和内部温度分布变化

2.2 弹丸外形影响

图4(a)和图4(b)分别给出了不同的弹丸头部外形，弹丸头部驻点烧蚀后退深度及其后退速度随时间的变化关系。从图4(a)可以看出，在t=1 s前，不同外形的弹丸头部驻点烧蚀后退深度相差不大；t=1 s之后，初始头部半径较大(算例3)的弹丸头部驻点烧蚀后退深度逐渐变小，略小于头部半径为2.5 mm(算例1，2)的弹丸，而头部锥角增大的弹丸(算例2)其驻点烧蚀后退深度增加。图4(b)的弹丸头部驻点烧蚀后退速度曲线表明，在烧蚀发生初期，其烧蚀后退速度急剧增加，随后保持相对稳定，其中算例3的烧蚀后退速度最小。

图4 不同弹丸形状驻点烧蚀随时间变化

文献[12]给出了高超声速弹丸的表面热流密度的经验关系式，如下

(7)

式中：u∞为弹丸飞行速度；ρ∞为表面气体的密度；N和M为系数，对于不同情况，M取值在3～3.7之间，N在0.5～0.8之间；C与弹丸其他条件相关，对于球头外形弹丸驻点，C=1.83×10-8R-1/2(1-hw/ho)，其中R为弹丸半径，hw和ho分别代表壁焓和总焓。在弹丸飞行速度和弹丸表面气体密度不变的条件下，式(7)可表达为

(8)

(9)

由此可知，当弹丸头部半径Rn减小或半锥角θc增大，都可以使热流密度ql增加，因此使得烧蚀深度和速度增加，这与图4数值计算的结果定性相符。

图5给出了不同弹丸外形(算例1，2和3)条件下，t=2 s时的弹丸烧蚀外形变化，总体来讲3种初始外形的弹丸，烧蚀后其外形形状变化相差不大。在该时刻，头部半径为3 mm(算例3)的弹丸头部驻点烧蚀深度相对初始头部长度(Ln)百分比为11%，头部半径为2.5 mm(算例1)的弹丸为12%，而半锥角增加为17°的弹丸(算例2)这一百分比则达到约15.4%。上述结果表明，头部半径增加会使烧蚀深度略微减少，防热效果更好；而弹丸头部半锥角增大，则会使烧蚀略有加剧。但是，弹丸头部半径的增大也会进一步增加头部的钝化程度，使得弹丸压差阻力增加，因此，需要设计者根据具体要求在气动热防护和低阻外形中寻找平衡。在弹丸其他外形参数确定的情况下，为减少弹丸烧蚀程度应当选择较小的头部半锥角。

图5 不同弹丸头部形状的烧蚀外形(t=2 s)

2.3 弹丸飞行速度影响

图6(a)、图6(b)分别给出了不同弹丸飞行速度条件下(算例1和算例4)，弹丸头部驻点烧蚀后退深度和烧蚀速率随时间变化。可以看出，马赫数为7的弹丸(算例4)在接近1s时才开始烧蚀，且最终驻点烧蚀后退深度远小于马赫数为8的弹丸(算例1)。在烧蚀初期都出现驻点烧蚀后退速度急剧上升的阶段，但马赫数为7的弹丸驻点烧蚀速度始终小于马赫数为8的弹丸。从式(7)也可以看出，弹丸表面热流密度与飞行速度的3次方成正比，因此，弹丸的飞行速度对弹丸所受气动加热有重要的影响，对弹丸的烧蚀起主要作用，这一点与本文计算结果亦相符。

图7为不同弹丸飞行速度(算例1，4)条件下，t=2 s时烧蚀后退外形变化。从图7可以看出，经过相同的飞行时间，马赫数为7的弹丸烧蚀后退深度较小，其头部驻点烧蚀后退深度约为0.98 mm，约占弹丸初始头部长度的3.2%，烧蚀现象不明显；而马赫数为8的弹丸出现明显的烧蚀现象，其头部驻点烧蚀后退深度是马赫数为7情况下的3.7倍。小马赫数时烧蚀较缓慢，在t=2 s时头部烧蚀呈平头状态，锥面变化不大。然而，尽管小马赫数下弹丸烧蚀不明显，但平头形状的弹丸也会影响弹丸的飞行气动参数，改变其外弹道飞行轨迹。

图6 不同弹丸飞行速度驻点烧蚀随时间变化

图7 不同弹丸飞行速度的烧蚀外形(t=2 s)

2.4 弹丸材料影响

以上计算分析的弹丸均是以金属钨作为弹丸材料，本节针对钢材料的弹丸进行数值模拟，并将结果与以钨为材料的弹丸结果进行对比。图8(a)、图8(b)分别给出了钨和钢两种材料弹丸的头部驻点烧蚀后退深度和烧蚀速率随时间变化，取钢的氧化物(Fe3O4)的熔点 1 867 K作为烧蚀温度。从其中可以看出，钢材料弹丸(算例5)比钨材料弹丸(算例1)更早开始发生烧蚀，且最终驻点烧蚀后退距离远超过钨材料弹丸。钢材料弹丸驻点烧蚀速度同样在开始烧蚀时急剧增加，且烧蚀速度高于钨材料弹丸。

图9为不同弹丸材料(算例1、5)条件下t=2 s时烧蚀后退外形变化，可以发现不同材料弹丸最终烧蚀外形不同：钨弹丸头部在t=2 s时呈现球形烧蚀外形，而钢弹丸头部烧蚀后则呈平头形状。图10给出了钢弹丸(算例5)在几个典型时刻下的烧蚀外形及其温度分布变化。与钨弹丸不同的是，钢弹丸高温区域始终主要分布于头部，锥面及弹丸内部升温缓慢，因此导致了钢材料弹丸始终保持平头形状后退，且烧蚀后退速度高于钨材料，这主要是由于材料本身导热和氧化物熔点等的不同所引起的。相比之下，钨材料的抗烧蚀性能更好，而钢材料弹丸的气动参数可能对弹丸的烧蚀形状产生不利影响。

图8 不同材料弹丸驻点烧蚀随时间变化

图9 不同弹丸材料的烧蚀外形(t=2 s)

图10 算例5的弹丸内部温度分布

3 结论

本文通过气动加热、温度分布和烧蚀外形变化的数值计算，模拟了弹丸飞行过程中的瞬时烧蚀情况，并计算分析了不同条件下弹丸的烧蚀情况。根据计算结果可以发现，弹丸烧蚀主要发生在头部，且随弹丸飞行速度的增加而愈发显著；当弹丸飞行速度一定时，弹丸烧蚀与弹丸形状有一定关系，其中弹丸头部半径对烧蚀影响不明显，而增大头部半锥角会加重烧蚀；钨材料弹丸的抗烧蚀特性优于钢材料弹丸。本文的研究对高超声速弹丸材料选择和外形设计具有一定的参考意义。由于弹丸烧蚀外形的改变，势必对其飞行过程中的气动参数和稳定性产生影响，进而影响其外弹道飞行轨迹，因此，下一步工作将在本文的基础上计算气动烧蚀对弹丸外弹道飞行特性的影响。

[1] 宋博,沈娟.美国的X-51A高超声速发展计划[J].飞航导弹,2009(5):36-40.

[2] 张攀峰,詹世革.从国家自然科学基金资助看高超声速流动研究的发展现状[J].航空学报,2015,36(1):1-6.

[3] 吕红庆,王振清,王永军,等.高超声速钝头体气动热分析[J].导弹与航天运载技术,2008(3):41-45.

[4] 李海燕,唐志共,杨彦广,等.高超声速飞行器高温流场数值模拟面临的问题[J].航空学报,2015,36(1):176-191.

[5] CHIANG T,HOFFMANN K,SIDDIQUI M.Analysis of nosetip recession of the tungsten projectiles in hypersonic flight[C]//31st Aerospace Sciences Meeting.Reno,NV,USA,1993.

[6] 史金光.高超声速弹翼表面气动烧蚀现象的研究及其仿真数值计算[D].南京:南京理工大学,2002.

[7] 杨晨.高超声速飞行器的气动烧蚀及内部传热研究[D].南京:南京理工大学,2008.

[8] 蒋友娣,董威,陈勇.高超声速钝头体变熵流表面热流计算[J].航空动力学报,2008,23(9):1591-1594.

[9] 张志豪,孙得川.飞行器气动加热烧蚀工程计算[J].兵工学报,2015,36(10):1949-1954.

[10] EWING M E,LAKER T S,WALKER D T.Numerical Modeling of Ablation Heat Transfer[J].Journal of Thermophysics and Heat Transfer,2013,27:615-632.

[11] TSE-FOU ZIEN.Recent Developments in Modeling Hypersonic Melting Ablation[C]//6th AIAA Theoretical Fluid Mechanics Conference,Honolulu,Hawaii,2011.

[12] 中国人民解放军总装备部军事训练教材编辑工作委员会.高超声速气动热和热防护[M].北京：国防工业出版社,2003.

NumericalSimulationsofAerodynamicAblationfortheHypersonicProjectile

CHEN Rui, DONG Gang

(National Key Laboratory of transient physics, Nanjing University of Science & Technology, Nanjing 210094, China)

The projectile can be subjected to the aerodynamic heating during the flight. For hypersonic projectiles, the aerodynamic heating causes the ablation of the projectiles and there leads to the change of the projectile shape. In present study, the heat transfer between flow and structure and the ablation process of a projectile were numerically simulated. Numerical results show that the ablation of stagnation point on the projectile is worst, and the degree of projectile ablation increases with the Mach number. The head radius of the projectile has no obvious influence on ablation, while reducing the semi-cone angle of projectile head can decrease the ablation. The ablation resistance of tungsten is superior to that of steel. The results in present study are practically significant for the thermal protection design of hypersonic projectile.

hypersonic; aerodynamic ablation; projectile; numerical simulation

2017-06-10；

2017-06-27

总装重点实验室开放基金项目(9140C300502130C30105)

陈蕊(1993—)，女，硕士研究生，主要从事燃烧爆炸与爆轰研究。

董刚(1970—)，男，研究员，主要从事新概念发射与推进理论和技术研究。

10.11809/scbgxb2017.10.007

本文引用格式：陈蕊,董刚.高超声速弹丸气烧动蚀数值模拟研究[J].兵器装备工程学报，2017(10):30-35.

formatCHEN Rui, DONG Gang.Numerical Simulations of Aerodynamic Ablation for the Hypersonic Projectile[J].Journal of Ordnance Equipment Engineering,2017(10):30-35.

TJ413

A

2096-2304(2017)10-0030-06

(责任编辑周江川)