1 数值算法

1.1 气动热数值算法

不考虑体积力和内热源情况下，在直角坐标系下的流体动力学N-S方程的积分形式为

式中：W为守恒变量；F_c为无粘通量；F_v为粘性通量。针对以上流动控制方程，采用有限体积法进行空间离散，对流通量采用AUSM+^[10]空间离散格式，该格式的数值耗散较小，激波分辨率较高，且鲁棒性较强。AUSM+格式将无粘通量分为对流项F_e和压力项P，在单元界面处其表达式为

AUSM+格式通过单元界面处的马赫数进行流场上、下游的判断：

式中：Ma_i+1/2,j,k为单元界面处的马赫数。而单元界面处的压力项则可分解为

为获得单调解，采用完全迎风的二阶MUSCL格式^[11]离散分裂后的无粘通量,并采用Minmod限制器使格式达到空间二阶精度。流体控制方程中的粘性项采用中心格式进行空间离散，湍流模型采用两方程Menter’s SST k-ω模型^[12]。针对非定常问题，采用三阶精度的显示多步Runge-Kutta方法^[13]求解，其在n~n+1时间步内,有

采用隐式残差光顺技术^[14]加速收敛，并且针对高超声速气动热问题，流体导热系数和粘性系数通常采用Sutherland公式或分子动力学计算，其对物面热流密度的计算精度有重要影响。由于本文的算例不考虑化学反应、热力学状态（内能松弛）等，以上控制方程和数值算法不涉及高超声速条件下的真实气体效应。

1.2 结构热分析数值算法

在无体积热源的假设下结构瞬态热传导的控制方程为

式中：ρ₀为结构材料密度；c为材料比热容；k_x,k_y和k_z分别为材料3个方向的导热系数。其中比热容和导热系数一般为温度的函数。针对本文热防护系统的热分析问题，其外表面边界条件为壁面热流密度Q_aero和壁面热辐射量Q_rad，其表达式分别为

式中：∂T/∂n为壁面法向温度梯度；ε为壁面热辐射率；σ为玻尔兹曼常数；T_wall为壁面温度；T_at为大气环境温度。对式(7)进行有限元离散可得到总体合成矩阵求解方程

式中：C为热容矩阵；K为热传导矩阵；P为温度载荷列阵；T和T˙分别为节点温度列阵和节点温度对时间的导数列阵。针对第n~n+1时间步，用Gaierkin格式离散得

求解式(11)即可得到结构各个时刻节点的瞬时温度，且本文结构传热采用ABAQUS有限元软件进行分析。

2 耦合分析模型及策略

锥体高速飞行时产生的巨大气动热会造成结构温度升高，而结构温度的升高导致激波层和边界层内的气体与壁面之间的温度梯度减小，这将导致壁面热流密度也将减小即气动热与结构传热之间存在着强烈的耦合效应，如图1所示。

图1 气动热与结构传热之间的耦合关系

Fig.1 Coupled relationship between aerodynamic heating and structural heat transfer

本文采用分区迭代推进的方法进行气动热与结构传热的分析，其推进求解策略如图2所示，流场与结构分为两个不同的区域，分开建模与求解，其中流场采用有限体积法（Finite volume method，FVM）求解，而结构热传导采用有限元法（Finite element method，FEM）求解，其特点为：

图2 分区迭代推进策略

Fig.2 Division iterative marching strategy

(1) 在任意n~n+1时间步内FVM（或FEM）求解过程中壁面温度（或壁面热流密度）不变，即各子学科求解时边界条件冻结。

(2) 各子学科在相同的时间节点进行数据交换，以保证耦合分析的协调性以及耦合时间精度。

本文气动热与结构传热的分区迭代推进方法的主要分析流程为：

(1) 首先建立数值分析模型，定义来流条件和壁面初始温度T₀，进行定常流场的计算，计算结果作为瞬态耦合分析的流场初始条件。

(2) 进行t_i时刻的非定常流场分析，将计算获得的热流密度Q_i通过插值方法传递给结构传热有限元模型外表面。

(3) 进行t_i时刻的结构有限元瞬态传热分析。

(4) 判断t_i是否达到瞬态分析的总时间t_total，若t_i达到t_total，结束耦合计算；若t_i未达到t_total，将计算获得的壁面温度通过插值方法传递给流体模型。

(5) 返回步骤(2)，进行下一个时间步t_i=t_i+Δt的求解，直至t_i=t_total，结束耦合分析。

3 数据传递方法

耦合分析中涉及到热流和壁面温度的数据传递，而实际分析中FEM网格尺度通常远大于FVM网格，即FVM和FEM耦合面节点不一致，如图3所示，故需要在耦合面进行插值以实现数据传递。本文采用如下的基于虚拟空间的插值方法进行热流密度的数据传递：

图3 FVM节点与FEM节点对应关系

Fig.3 Nodes relationship between FVM and FEM

(1) 将各学科耦合面上节点从物理空间(x, y, z)通过坐标变换映射到二维虚拟空间(u, v)，x=x(u, v)，y=y(u, v)和z=z(u, v)，即将三维曲面上的空间节点转换到二维虚拟空间平面上。

(2) 在物理空间中搜索任意FEM网格节点ζ_i(x, y, z)附近的FVM网格节点η_i(x, y, z)，并将其转换到虚拟空间得到ζ_i(u, v)和η_i(u, v)。

(3) 将FVM网格节点坐标η_i(u, v)和相应的热流Q_i(u, v)代入到二次插值函数Q(u, v)，采用最小二乘法求解插值函数的系数，即

(4) 将所有FEM网格节点ζ_i(u, v)代入到已知系数的插值函数Q(u, v)中，即可求得FEM网格节点插值热流密度q_i(u, v)。

4 验证算例

采用NASA进行的高超声速圆管风洞试验模型^[15]验证了本文的分区迭代推进分析方法的正确性。圆管内径R₁=25.4 mm，外径R₂=38.1 mm，圆管材料为不锈钢，其密度ρ=8 030 kg/m³，导热系数k=16.72 W/(m·K)，比热容c=502.48 J/(kg·K)，此外圆管初始壁面温度T₀=294.4 K。来流马赫数Ma=6.47，来流温度T=241.5 K，来流压强p=648.1 Pa，迎角α=0°。二维分析模型如图4所示，其中流场网格量为60×100，壁面第1层网格高度Δh=1×10^-5 m，并且在激波处进行了局部加密，而流场分析对流项采用有限体积法AUSM+空间离散格式，湍流模型选择Menter’s SST k-ω两方程模型。圆管结构有限元网格量为40×25，采用四节点平面单元模拟。采用基于虚拟空间的插值方法进行耦合数据的传递，分析类型为瞬态分析，耦合分析时间步长Δt=1×10^-4 s，分析总时间t_total=2 s，并且将定常流场作为耦合分析的初始条件。分析获得了圆管驻点温度随时间的变化曲线，如图5所示，从图5中可观察到曲线斜率在初始阶段很大，随时间推移逐渐减小，说明结构开始阶段升温很快，随后温度上升速度逐渐减缓，最后趋于一个稳定的温度值（稳态解）。图6为2 s时刻流场和圆管结构的温度云图，从图6可观察到流场数值解捕捉到了光滑的弓形激波，且圆管最高温度出现在驻点处，沿壁面下游温度逐渐降低。

图4 流体和结构模型网格

Fig.4 Meshes for fluid and structural models

图5 驻点温度随时间变化情况

Fig.5 Temperature of stagnation point varying with time

图6 2 s时刻流体和结构温度场

Fig.6 Temperature fields for fluid and structural models at 2 s

图7为2 s时刻圆管外壁面的相对热流密度分布和相对温度分析情况。从图7中可知壁面耦合分析结果分布曲线与试验结果分布曲线吻合良好，其中驻点处热流密度耦合分析结果Q_stag为661 kW/m²，而试验结果为670 kW/m²，相对误差为1.34%；而驻点处温度耦合分析结果T_stag为442 K，而试验结果为465 K，相对误差为4.95%。

图7 2 s时刻壁面热流密度和温度分布

Fig.7 Distribution for wall heat flux and temperature at 2 s

由以上讨论可知采用本文耦合推进分析方法得到的圆管壁面热流密度及温度分布均与试验结果吻合良好，从而验证了本文提出的分区迭代推进分析方法的正确性与分析精度，可用于分析锥体高速飞行时的热环境问题。

5 锥体热环境分析

5.1 初始气动热分析

分析模型来自于NASA技术报告 NASA TP 2334^[16]，其几何形状为直二次圆锥体，如图8所示。来流马赫数Ma=9.86，来流静压p=59.92 Pa，来流静温T=48.88 K，迎角α=0°，壁面温度T_wall=300 K。划分了CFD结构网格，如图9所示，六面体网格总量约50万，壁面第1层网格高度小于1×10^-5 m，以保证量纲一参数y⁺≤1。进行了初始定常流场求解，空间离散格式采用二阶精度的AUSM+格式，湍流模型采用Menter’s SST k-ω模型。

图8 分析模型几何外形及尺寸

Fig.8 Geometry and sizes of analysis model

图9 CFD模型的网格及边界条件

Fig.9 Mesh and boudary conditions of CFD model

获得的流场压力云图和温度云图如图10所示。数值解捕捉到了光滑的离体激波，锥体驻点处的压强超过了7 000 Pa，而正激波后气体温度则超过了950 K，与300 K的初始壁面温度形成强烈的温度梯度，产生很高的热流密度值。锥体表面压力分布情况如图11(a)所示，而锥体表面热流密度计算值和试验值分布情况如图11(b)所示，驻点处热流密度试验值为443.2 kW/m²，而本文的数值计算值为429.4 kW/m²，与试验值的相对误差为3.1%，除了驻点以外的其他区域热流密度计算值与试验值均吻合得很好，这验证了本文分析模型的正确性。

图10 头部压力云图和温度云图

Fig.10 Pressure and temperature fields of head area

图11 物面压强和热流分布

Fig.11 Pressure and heat flux distributions of wall

5.2 分区迭代分析

结构瞬态传热分析的有限元模型如图12所示。结构材料为不锈钢，密度ρ=8 030 kg/m³，导热系数k=16.27 W/(m·K)，比热容c=502.48 J/(kg·K)，结构初始温度T₀=300 K。采用以上CFD分析获得的初始流场进行气动热和结构传热的瞬态耦合分析，耦合方法为串行耦合迭代方法，耦合面数据传递采用虚拟空间插值法，耦合时间步Δt=0.000 1 s，瞬态分析总时间t_total=100 s。

图12 结构传热分析有限元网格

Fig.12 Finite element mesh for structural heat transfer analysis

分析获得了锥体结构前缘点（流场驻点）处的温度随时间变化的情况，如图13所示。从图13中可以观察到随着时间的推移，结构前缘点温度逐渐上升，这是由于外部连续的气动热作用于结构物面，导致结构温度持续上升。此外，随着时间的推移前，缘点温度上升的趋势逐渐减小，即曲线斜率逐渐减小，时间趋于无穷大时斜率将趋于0，结构温度场达到稳定（稳态解）。

图13 驻点处温度随时间的变化情况

Fig.13 History of heat flux at stagnation point

图14为25,50,75和100 s时刻结构剖面的温度云图，从图14中可观察到各时刻结构最高温度分别为498.0,546.0,576.1和597.9 K，结构最高温度均在锥体前缘驻点处，且随着时间的推移高温区域越来越大。分析获得了不同时刻锥体表面压强和热流密度的分布情况，如图15所示，从图15中可观察到时间的变化对锥体表面的压强几乎没有影响，而表面热流密度却随时间的增加而降低。这主要是因为随着时间的推移锥体表面温度逐渐升高，激波层及边界层与壁面之间的温度梯度将逐渐减小，这就造成了锥体表面热流密度的降低；而壁面压强主要受控于流体力学而非热力学，故壁面压强几乎不改变。

图14 不同时刻结构剖面的温度云图

Fig.14 Temperature field of structural profile at different anlysis time

图15 不同时刻壁面压力和热流分布

Fig.15 Pressure and heat flux distributions at different time

6 结 论

(1) 本文提出了气动热与结构传热的分区迭代推进分析方法，其考虑了气动热和结构传热之间的耦合效应，适用于本文所研究的高超声速钝头体热环境的精确分析。而传统的非耦合分析方法未考虑结构温度升高对气动热的反馈效应。

(2) 进行了圆管验证算例分析，2 s时刻圆管驻点处的热流密度和温度的计算值与试验值的相对误差分别为1.34%和4.95%，从而验证了本文分区迭代推进分析方法的合理性与分析精度。

(3) 最后研究了直二次圆锥体的热环境，锥体驻点温度随时间的推移而升高，且上升趋势逐渐变缓，最终趋于稳态值。此外时间的变化对锥体表面的压强几乎没有影响，而表面热流密度却随时间的增加而降低。

参考文献