高超声速导弹多场耦合仿真

2019-08-13徐世南吴催生

宇航学报 2019年7期

徐世南，吴催生

(中国空空导弹研究院，洛阳 471000)

0 引言

新一代高超声速导弹采用超燃冲压发动机等新型动力，从而实现短时间内打击长距离军事目标[1]。为了保证燃料供给充分，又不增大飞行阻力，弹体结构需要在减小弹径的同时，增加发动机舱长度，因此大长细比成为高超声速导弹结构形式的一个发展方向。此结构形式保证导弹实现高速度、超机动性和广泛作战空间的同时，也增大了结构的柔性，导致结构在气动载荷作用下会产生更大的变形，气动与结构间的耦合效应显著增强[2-5]。

由于飞行速度快，高超声速导弹面临更加严峻的热环境，结构刚度的改变造成结构的承载能力降低，飞行器内部温度的增高致使元器件等部件损害率加大，这也是目前高超声速飞行器面临的共性问题。面对热防护设计和结构设计带来的挑战，众多试验项目以此为背景相继进行：北京航天长征飞行器研究所采用主动式热疏导技术对金属尖化前缘模型进行试验研究[6];吴大方等[7]对某飞行器翼面进行了热/振联合试验;成都飞机设计研究所对飞行器热环境下的结构固有振动特性进行了试验[8];中国空气动力研究所与发展中心分别进行了MF-1航天模型飞行试验以及激波风洞试验[9-10]。但是试验要求技术难度大，经费高，采用数值仿真对飞行器气动力、结构温度等进行有效的预测，可为高超声速导弹热防护设计与结构设计提供指导。

由于耦合效应不可忽视，对高超声速导弹流场与结构温度仿真研究需采用多场耦合方法。目前关于耦合数值计算方法，可以分为两大类：一体化求解方法和分区求解方法。一体化求解方法也被称为直接耦合求解方法，是将不同子系统的控制方程放在一个系统方程组中求解[11]。黄唐等[12]进行了二维流固热一体化数值模拟，以二维圆柱绕流问题为例，流场和热响应部分分别采用总变差消格式和迦辽金方法进行方程离散。季卫栋等[13]将流场与结构控制方程写成统一形式，对整个物理场进行统一的迎风格式有限体积方法离散，实现了二维圆管模型的一体化数值模拟。虽然一体化求解中流体与结构的解变量完全耦合，求解精度更高，但由于流体和结构的控制方程具有完全不同的物理特性和数学性质，目前仅在简单的二维问题展开研究。分区求解方法是将系统分为流体和固体两个独立的模块[14]。刘健等[15]、周印佳等[16]与黄杰等[17]均采用分区耦合方法，对高超声速二维圆管进行数值模拟，计算结果与试验值相吻合。夏刚等[18]将新型上风格式与结构传热的Galerkin方法相结合，实现钝体高超声速飞行器的分区耦合计算。张胜涛等[19]采用自适应耦合时间步长和混合差值策略实现了高超声速飞行器进气道前缘结构的分区耦合仿真分析。肖军等[20]通过在子代步上气动与结构方程交替求解实现了一种全隐式分区耦合算法。李凯伦等[21]基于有限元数值分析方法建立了一种多场分区耦合模型，对高超声速环境中功能梯度薄板热气动弹性进行了研究。梁杰等[22]基于一种分区求解策略，采用直接模拟蒙特卡洛方法，对某飞行器气动力、气动热进行了数值模拟。

目前关于高超声速多场耦合的研究主要在耦合方法的实现上，但是顺利完成高超声速导弹的设计，不仅需要有效的多场耦合仿真技术去实现气动热力学环境预示，还需对各场之间耦合关系定量化分析，否则会造成设计不当而引发导弹高速飞行时结构破坏等不利后果。通过研究高超声速导弹温度场、压力场与结构变形之间的相互关系，分析耦合效应对其影响，对导弹热防护以及结构设计具有指导意义。

本文基于分区耦合方法，以大长细比导弹为研究对象，实现其流/固/热数值模拟。对导弹不同攻角下的结构变形和温度场、压力场进行计算分析，并对导弹气动力、气动热与结构变形间的耦合效应进行定量化研究。

1 耦合分析策略和控制方程

1.1 耦合方法

耦合模型如图1所示。在流体介质内部，求解统一的流体控制方程得到热流qf和压力pf，并通过流-固耦合界面向固体介质提供热载荷qf和力载荷pf。在固体介质内部，求解热传导控制方程和热弹性力学控制方程得到结构壁面温度Ts和位移us，并通过流-固耦合界面向流体介质提供温度条件Ts和结构变形条件us。

基于多场耦合模型，采用分区求解方法建立如图2所示的紧耦合分析策略，流体与固体区域求解器均为瞬态求解，求解器所需数据通过耦合界面数据不断交换提供，耦合分析策略实施流程描述如下：

1)首先根据结构的初始值，即结构壁面初始温度Ts与位移us，作为瞬态耦合分析的初始条件，通过耦合界面传递给流体域。

2)从流体域内进行瞬态流场计算，得到热流qf和压力(流体应力)pf通过耦合界面传递给固体域。

3)固体域内进行瞬态求解得到结构壁面温度Ts和位移us通过耦合界面传递给流体域。

4)检查流体应力和结构位移是否满足收敛标准，如果不满足，则返回第2步继续进行迭代求解，直到满足收敛标准或者达到设置的迭代次数；如果满足收敛标准，则t0时间步的计算结束，输出流体和结构的结果，进行t1时间步的耦合计算。

5)该过程不断重复直到计算时间结束。

1.2 流体控制方程

流体控制方程采用积分形式的N-S方程，其形式如下[23]：

(1)

式中：Q为守恒变量，F为无黏通量，FV为黏性通量，t为物理时间，∂Ω为某一固定区域Ω的边界，dS为面积微元，n为控制边界法向单位矢量。

1.3 热-结构控制方程

结构传热方程计算采用Fourier定律，不考虑内热源，并假设无辐射，基于Fourier定律的三维瞬态热传导控制微分方程为Fourier定律的三维瞬态热传导控制微分方程[23]：

(2)

式中：ρ为密度，c为比热容，T为温度，t为时间，λ为导热系数。

结构方程为热弹性力学控制方程，采用刚度矩阵求解法，假设结构产生的静变形是一个缓慢的过程，不考虑惯性力，其方程为[23]：

(3)

K是热结构刚度矩阵，热-结构控制方程通过K实现耦合，us是结构位移向量，Fs是作用在结构上的力矩阵。温度的改变会造成材料特性改变，变化后的结构刚度矩阵为KT；温度还会使结构内部产生热应力，结构间的相互约束使其受热后膨胀受到限制引发预拉压应力，这种热应力产生的附加刚度为KΔT。

2 仿真校验

经典圆管绕流试验作为算例[24]。试验所用圆管为不锈钢，圆管内、外半径分别为25.4 mm和 38.1 mm，热力学参数：导热系数16.27 W/(m·℃)，比热容502.48 J/(kg·℃)，密度8030 kg/m3，弹性模量1.2×1013Pa，热膨胀系数1.68×10-5(1/℃)，泊松比0.3。圆管内壁设为等温壁，温度值与初始环境温度一致，为21.4 ℃。来流参数：温度-31.5 ℃，压力648 Pa，速度Ma6.47。图3为二维计算网格，流场与圆管网格交界面为耦合界面，流固通过耦合界面进行数据传输实现多场耦合。图4和图5分别为计算初始时圆管表面压力和热流分布与文献[24]的试验结果对比，图中p0，q0分别为驻点压强与热流，θ为物面到圆心的连线与X轴的夹角。可以看出压p强和热流q分布与文献[24]符合得很好。图6为驻点2 s时刻的温度值与文献的对比，对比结果较好。由此，校验了此耦合方法的有效性。

3 计算结果与分析

3.1 计算模型

导弹几何模型如图7所示，仅考虑弹体结构，简化弹翼、舵面和内部元器件等。导弹结构材料采用钛合金，材料参数如表1所示。在弹身前段迎风和背风位置分别布置监测点p1，p2。导弹头部顶端为坐标原点。

计算模型与网格示意图分别如图8、图9所示。仿真计算模型流场区域分为外流场和内流场，外流场为导弹外部气流，内流场为导弹结构内部气体。外流场对称面采用对称边界条件，其它面采用外部流边界条件，在此边界条件上设置压力、速度和温度载荷，载荷值分别为5529 Pa，Ma5，-56 ℃。固体结构场为弹体结构，在弹体尾部采用固定边界条件。分别在外流场与弹体外壁面交界处，内流场与弹体内壁面交界设置外场与内场耦合边界条件，通过耦合边界实现分区耦合的数据传输。计算攻角分别为0°，5°，10°，15°，20°，仿真时间40 s，初始内、外流场和结构温度假设为20 ℃，压力值5529 Pa。

弹性模量/GPa热膨胀系数/(℃-1·10-6)热传导系数/(W·m-1·℃-1)比热容/(J·kg-1·℃-1)75～1099.20～11.066.8～11.8610～702

为研究导弹气动力、气动热与结构变形之间的相互关系，将计算工况分为耦合计算(考虑耦合效应影响)、仅考虑热因素(结构仅有热变形)的计算、仅考虑气动力因素(结构仅有气动变形)的计算和无结构变形因素计算四种情况。

3.2 结构变形分析

导弹结构在X，Y，Z三个方向发生变形，以攻角20°的位移云图为例，如图10所示。由图10可知，弹体结构发生拉伸、膨胀和弯曲变形。导弹为轴对称结构且仿真计算未考虑侧偏角，弹体结构在Y方向相对轴线仅发生膨胀而无弯曲变形，所以主要研究弹体结构轴向(X方向)和横向(Z方向)变形。

弹体结构沿导弹轴向的变形量为Δx=(Δx1+Δx2)/(2L)×100%，其中L为导弹全长，Δx1为结构迎风侧沿X方向的位移改变量，Δx2为结构背风侧沿X方向的位移改变量。弹体结构轴向变形曲线如图11所示，横坐标为弹体结构所处位置，x=0为导弹头部顶端位置。

由图11可知，导弹高速飞行时产生气动加热，使结构产生轴向的拉伸，头部气动加热最剧烈，拉伸量大。随着攻角增加，轴向变形量减小，因结构沿轴向拉伸时，挠度随攻角的增加而增加，相应沿轴向变形量减小。

弹体结构沿导弹横向的变形量为Δz=(Δz1+Δz2)/(2R)×100%，其中R为导弹半径，Δz1为结构迎风侧沿Z方向的位移改变量，Δz2为结构背风侧沿Z方向的位移改变量。其结构横向变形曲线如图12所示。

由图12可知，0°攻角时温度和气动力均布，弹体结构未发生轴向变形，随着攻角增大，结构发生弯曲变形，且头部变形量大。气动力/热共同作用使结构产生弯曲变形：弹体末端采用固定边界条件，类似于悬臂梁模型，梁承受Z方向的气动力使结构发生横向弯曲；同时气动加热作用使结构发生轴向拉伸，结构迎风侧温度高，拉伸量大，结构背风侧温度低，拉伸量小，造成弹体结构沿背风侧位置发生弯曲；气动力和气动加热在头部作用明显，弯曲变形量大。

研究气动加热与气动力对导弹结构变形的影响，以弹身前段位置为例，变形量由p1和p2测点数据计算而得，其中对结构轴向变形的影响见图13，对结构横向变形的影响见图14。

由图13可知，弹体结构的轴向变形主要由气动加热引起；气动力使弹体结构产生轴向压缩但影响很小，可忽略不计。由图14可知，弹体结构在横向的弯曲变形由气动热和气动力共同引起，且导弹采用大长细比结构时，结构柔度增加，气动力造成结构轴向的弯曲变形影响更大。

3.3 温度分析

导弹温度云图如图15所示。在弹体轴向位置，来流在头部位置发生摩擦和压缩最剧烈，动能转化为热能最多，气动加热剧烈，且随着攻角增大，温度最大值点从头部最前端往头部迎风侧后侧偏移。攻角为0°时温度均布，随着攻角增大，迎风侧与背风侧位置温度差也随之增大。

研究耦合效应对温度仿真结果影响，以弹身前段迎风侧p1监测点为例，如图16所示。当攻角较小时，耦合效应对对温度仿真结果影响不大；随着攻角增大，耦合效应对仿真结果的影响增大，因为导弹大攻角飞行时，弹体结构因气动热和气动力发生弯曲变形，攻角由α增加为α+Δα，且此Δα较大造成弹体结构的温度分布发生改变，加剧了弹体结构的气动加热。