吸气式高超声速飞行器的研究已成为当今世界航空航天发展的一个重要方向，但同时，与之伴随的气动热防护问题也成为一项亟待解决的重要技术难题[1]。如图 1所示，以超燃冲压发动机为动力的高速飞行器在大气层中以高马赫数飞行时，将面临十分严峻的气动热力学环境。特别是，其前体产生的激波系与进气道唇缘弓形激波发生干扰，形成激波-边界层、激波-剪切层、激波-激波相互作用等非线性物理现象[2]。激波干扰产生的复杂流场波系结构，在高速飞行器表面产生局部高温高压载荷，给飞行器的性能和结构带来极大安全隐患，同时也给高超声速飞行器热防护问题提出了很大挑战。
|
| 图 1 前体-进气道激波干扰示意图 Fig. 1 Schematic diagram of forebody-inlet shock wave interference |
Edney[3]根据激波干扰点的位置，将激波干扰分为6种类型(Ⅰ~Ⅵ型)，其中Ⅳ型干扰形成超声速射流冲击壁面，产生极大的压强和热流[4-5]，对飞行器安全影响最大，因此备受关注，也是本文主要研究内容。另外，高速飞行器外部流场的气动加热与热防护结构内的热传导相互作用不可分割，因此要准确预测“Ⅳ”型激波干扰下飞行器的气动热环境，必须考虑气动加热与结构传热的相互耦合作用[6]。
目前，国内外学者针对高超声速流-固-热多场耦合问题已经开展了大量研究工作。例如，Thornton和Dechaumphai[7]基于有限元方法提出了一种双向耦合计算方法来求解不锈钢平板在高超声速流动下的气动加热传热与变形问题，计算结果证明了多场耦合计算的重要性。Kazemi-Kamyab等[8]在分区双向耦合的基础上发展了一种强耦合迭代计算方法，采用高阶隐式时间迭代方法来求解流场与结构非稳态共轭传热计算问题，从而提高时域上的计算精度。Chen等[9]发展了一种高超声速流场-热-结构耦合分析平台(HyCCD)，分析了高超声速前缘在高温高压条件下的热力学行为，同时引入自适应时间步长计算方法，提高计算效率。Qin等[10]利用MPCCI多平台耦合计算软件将FLUENT与ABAQUS两种商用软件进行联合使用，计算模拟了高超声速条件钝体外形在不同支杆外形下流-热-固多场耦合特性，分析了不同外形支杆对高超声速钝头体的减阻防热效果。Zhao等[11]依据结构热响应特征时间较小的特点，发展了一种新的多物理场松耦合计算方法，数值模拟了马赫数为6.47下三维非定常圆管流-热-固耦合算例，得到的热壁热流、冷壁热流以及圆管温度与试验值吻合。Zhang等[12]基于自动控制理论发展了一种基于PID(Proportional-Integral-Differential)控制的自适应时间步长的松耦合高效计算方法，用于计算分析高超声速流中的共轭传热问题。
然而需要指出的是，上述计算方法都归属于多场分区耦合迭代计算方法，即将流场和固体结构独立分开建模，基于“准稳态耦合”策略，在耦合交界面进行数据传递与交换，同时在时间域上进行耦合交替迭代。后来，NASA兰利研究中心的Wieting等[13]认为多物理场分区耦合计算方法将连续的物理场进行分割建模，耦合界面需要额外的数据插值方法和交换策略，这种计算方法虽然较容易实现，但无法准确模拟流-热-固多物理场耦合特性。李鹏飞和吴颂平[14]在对高速飞机机身结构进行气动加热与结构传热的耦合数值计算分析后，认为为了准确分析热壁对传热计算的影响，很有必要开展气动加热与结构传热一体化数值计算方法研究。姜贵庆等[15]为了满足多场耦合计算问题中连续变化的物理条件，提出了基于有限元方法的流-热-固一体化计算方法，流场与结构温度传热计算采用统一网格和计算程序，计算结果表明一体化计算可以保证热结构的计算精度。
本文基于多物理场连续变化的物理条件，发展了一种高超声速流-热-固一体化计算方法。该方法基于有限体积法，将高超声速流场与结构温度场统一到同一控制方程中，并进行统一离散与统一求解，避开了传统分区耦合方法所需的数据传递策略。同时，提出一种新的双温阻模型以保证交界面上物理场的连续性。选取典型高超声速绕流二维钝体算例验证本文方法的可靠性与正确性，最后采用本文计算方法对高超声速“Ⅳ型”激波干扰开展流-热-固一体化计算与分析研究。
1 一体化数值计算方法 1.1 一体化控制方程及求解方法首先推导一体化控制方程。固体结构传热控制方程积分形式(不含热源，暂忽略热辐射换热)可写为
| $ \frac{\partial }{{\partial t}}\int_\mathit{\Omega } {\rho {C_{\rm{s}}}T{\rm{d}}\mathit{\Omega }} - \oint_{\partial \mathit{\Omega }} {k\nabla T \cdot \mathit{\boldsymbol{\hat n}}{\rm{d}}S} = 0 $ | （1） |
式中: ρ、Cs、T、k分别为密度、固体材料比热容、温度和热传导系数；Ω为控制体体积；dS为控制体表面面积微元；
对于结构外部流场计算，采用积分形式的可压缩Navier-Stokes方程作为控制方程，方程形式与式(1)相似，因此，流场计算方程与固体结构传热方程可写成如下统一控制方程：
| $ \frac{\partial }{{\partial t}}\int_\mathit{\Omega } \mathit{\boldsymbol{W}} {\rm{d}}\mathit{\Omega } + \oint_{\partial \mathit{\Omega }} {\left( {{\mathit{\boldsymbol{F}}_{\rm{c}}} - {\mathit{\boldsymbol{F}}_{\rm{v}}}} \right){\rm{d}}S} = {\bf{0}} $ | （2） |
式中：W为守恒量；Fc为对流通量；Fv为黏性通量。W、Fc、Fv的表达式为
| $ \mathit{\boldsymbol{W}} = \left[ {\begin{array}{*{20}{c}} \rho \\ {\rho u}\\ {\rho v}\\ {\rho w}\\ {\rho E} \end{array}} \right],{\mathit{\boldsymbol{F}}_{\rm{c}}} = \left[ {\begin{array}{*{20}{c}} {\rho V}\\ {\rho uV + {n_x}p}\\ {\rho uV + {n_y}p}\\ {\rho wV + {n_z}p}\\ {\rho HV} \end{array}} \right] $ |
| $ {\mathit{\boldsymbol{F}}_{\rm{v}}} = \left[ {\begin{array}{*{20}{c}} 0\\ {{n_x}{\tau _{xx}} + {n_y}{\tau _{xy}} + {n_z}{\tau _{xz}}}\\ {{n_x}{\tau _{yx}} + {n_y}{\tau _{yy}} + {n_z}{\tau _{yz}}}\\ {{n_x}{\tau _{zx}} + {n_y}{\tau _{zy}} + {n_z}{\tau _{zz}}}\\ {{n_x}{\mathit{\Theta }_x} + {n_y}{\mathit{\Theta }_y} + {n_z}{\mathit{\Theta }_z}} \end{array}} \right] $ |
| $ \left\{ \begin{array}{l} {\mathit{\Theta }_x} = u{\tau _{xx}} + v{\tau _{xy}} + w{\tau _{xz}} + k\frac{{\partial T}}{{\partial x}}\\ {\mathit{\Theta }_y} = u{\tau _{yx}} + v{\tau _{yy}} + w{\tau _{yz}} + k\frac{{\partial T}}{{\partial y}}\\ {\mathit{\Theta }_z} = u{\tau _{zx}} + v{\tau _{zy}} + w{\tau _{zz}} + k\frac{{\partial T}}{{\partial z}} \end{array} \right. $ |
上述控制方程用于流场计算时，u、v、w分别为控制体3个方向的速度；V=nxu+nyv+nzw为面法向速度，其中nx、ny、nz为面法矢分量；τij为应力张量分量；H为总焓，H=E+p/ρ，其中E为流体控制体总能, p为压力。另外控制方程满足理想气体状态方程p=ρRT，R为理想气体常数。对于固体结构传热计算时，结构无变形满足u=v=w=0，式(2)中对流通量项Fc =0，黏性通量项中，Θx=k∂T/∂x，Θy=k∂T/∂y，Θz=k∂T/∂z; E=CsT为固体结构单元控制体内能。
采用有限体积法对上述一体化控制方程进行离散求解。其中采用AUSM+[16]迎风格式求解对流通量项，同时基于Venkatakrishnan[17]限制器采用MUSCL[18]插值方法进行线性重构，以实现二阶精度计算。采用二阶中心差分格式计算黏性通量项。时间推进采用LU-SGS(Lower-Upper Symmetric Gauss-Seidel)隐式时间迭代计算方法，其中流场当地时间步长计算公式为
| $ \Delta t_{\rm{f}}^I = {N_{{\rm{cfl}}}}\frac{{{\mathit{\Omega }_I}}}{{{{\left( {{\mathit{\Lambda }_{\rm{c}}} + C{\mathit{\Lambda }_{\rm{v}}}} \right)}_I}}} $ | （3） |
式中: Ncfl为CFL数；Λc和Λv分别为控制体I单元所有边界的对流谱半径和黏性谱半径的累加；ΩI为第I个单元的体积；C为常数。
针对结构传热方程的计算时间步长，需要注意的是，结构传热过程中产生的热脉冲需要限制在当地固体结构单元中，才能保证计算的稳定性。因此，结构传热方程当地时间步长计算公式为
| $ \Delta t_{\rm{s}}^I = \frac{{{{\left( {{\rm{ \mathsf{ δ} }}x} \right)}^2}}}{{2\alpha }},\alpha = \frac{k}{{\rho {C_{\rm{s}}}}} $ | （4） |
式中：δx为单元最小边长; α为热扩散率。
另外，为保证壁面热流计算的准确性，需要选取合适的湍流模型。对多种湍流模型进行对比分析后，确定选取Menter[19]提出的SST(Shear Stress Transport)k-ω两方程湍流模型，该模型不需要阻尼函数，能较好计算边界层近壁面热通量。
1.2 流-固耦合界面参数定义由于外部流场与固体结构热力学特性差异较大，本文一体化计算方法为了保证流-固耦合界面热通量与温度的连续性，需要对交界面的热力学参数进行特别计算。
图 2给出了流-固交界面左右单元意图，Ωf和Ωs分别表示流体单元(左边)和固体结构单元(右边)。耦合界面的温度计算式为
| $ T = \frac{{{k_1}{T_1} + {k_{\rm{r}}}{T_{\rm{r}}}}}{{{k_1} + {k_{\rm{r}}}}} $ | （5） |
|
| 图 2 流-固交界面左右单元意图 Fig. 2 Schematic diagram of left and right cells at fluid-solid interface |
式中：下标l、r代表耦合界面左、右单元。
交界面上温度梯度▽T的准确计算对热流的准确预测起着至关重要的作用，本文温度梯度▽T的计算方法采用高斯格林[20]方法，并采用如下方法进行修正[21]：
| $ \left\{ {\begin{array}{*{20}{l}} {\nabla {T^*} = \left( {\nabla {T_1} + \nabla {T_{\rm{r}}}} \right)/2}\\ {\nabla T = \nabla {T^*} - \left( {\nabla {T^*} \cdot {\mathit{\boldsymbol{r}}_{{\rm{lr}}}} - \frac{{{T_{\rm{r}}} - {T_1}}}{{{L_{{\rm{lr}}}}}}} \right){\mathit{\boldsymbol{r}}_{{\rm{lr}}}}} \end{array}} \right. $ | （6） |
式中: Llr为左右单元中心之间的距离; rlr为单位向量。
值得注意的是，流-固交界面的热传导系数的计算对一体化传热计算起到关键作用。固体结构内部导热系数较大，交界面的热传导系数采用简单加权平均会对计算产生很大误差，为了提高计算的准确性与合理性，本文提出了一种双-热阻抗模型用于计算交界面热传导系数。其中热阻抗表征的物理意义为热量在介质中传导时所遇到的阻力，其定义为
| $ {R_{\rm{t}}} = \frac{L}{{kA}} $ | （7） |
式中:L为热通量传导的距离; A为热量传导的截面面积。因此，可以得到交界面的热阻关系式为
| $ {R_{{\rm{t}},{\rm{bnd}}}} = {R_{{\rm{t}},1}} + {R_{{\rm{t}},{\rm{r}}}} $ | （8） |
其中：Rt,l和Rt,r分别为交界面左右单元的热阻。将式(7)代入式(8)中，得到
| $ \frac{L}{{kA}} = \frac{{{L_1}}}{{{k_1}A}} + \frac{{{L_{\rm{r}}}}}{{{k_{\rm{r}}}A}} $ | （9） |
因此，可得交界面当量导热系数k为
| $ k = \frac{{{k_1}{k_{\rm{r}}}}}{{{k_1}{L_{\rm{r}}} + {k_{\rm{r}}}{L_1}}} \cdot \left( {{L_{\rm{r}}} + {L_1}} \right) $ | （10） |
为了提高非定常计算效率，本文非定常计算中，真实物理时间步长采用于经典PID控制器[22]调节的自适应控制算法进行实时计算，以取代传统固定时间步长，伪时间步长选取当地时间步长。
图 3为自适应时间步长PID控制器示意图，整个自适应控制器由处理器和控制器2个部分组成。自适应控制基本原理为：处理器将时间步长Δt作为输入控制参数，进行数值迭代输出误差值e；控制器根据输出误差值e和给定输入的容错值tol进行反馈控制，计算合适的新时间步长Δt。
|
| 图 3 时间步长PID控制器示意图 Fig. 3 Schematic diagram of PID control system in time step-size |
通过PID控制器得到的真实物理时间步长Δt计算为
| $ \Delta t_{n + 1}^{\rm{r}} = {\left( {\frac{{{e_{n - 1}}}}{{{e_n}}}} \right)^{{k_{\rm{P}}}}}{\left( {\frac{1}{{{e_n}}}} \right)^{{k_{\rm{I}}}}}{\left( {\frac{{e_{n - 1}^2}}{{{e_n}{e_{n - 2}}}}} \right)^{{k_{\rm{D}}}}}\Delta {t_n} $ | （11） |
式中:n为计算迭代时间步；kP、kI、kD为PID控制器参数, 分别取值为kP =0.075，kI =0.175, kD =0.015；en为温度误差值，其计算式为
| $ {e_n} = \frac{{e_n^*}}{{{\rm{tol}}}},e_n^* = \frac{{\left\| {{T^n} - {T^{n - 1}}} \right\|}}{{\left\| {{T^n}} \right\|}} $ | （12） |
其中:tol为给定值，本文给定tol=1.0。
为了保证PID控制器的稳定性和计算迭代的收敛性，需要对时间步长的大小和增长率进行一定约束：
| $ \Delta {t_{\min }} \le \Delta {t_n} \le \Delta {t_{\max }} $ | （13） |
| $ \Delta t_{n + 1}^\alpha = \frac{{{\alpha _{{\rm{ref}}}}}}{\alpha }\Delta {t_n} $ | （14） |
| $ \Delta {t_{n + 1}} = \min \left( {\Delta t_{n + 1}^{\rm{r}},\Delta t_{n + 1}^\alpha } \right) $ | （15） |
式中：Δtmin和Δtmax分别为设定的最小和最大的非定常计算时间步; αref为时间步长控制率，一般0.2＜αref＜0.4，理想取值αref≈0.2;α的表达为
| $ \alpha = \mathop {\max }\limits_n \frac{{\left\| {{T^n} - {T^{n - 1}}} \right\|}}{{\left\| {{T^{n - 1}} - {T^{n - 2}}} \right\|}} $ | （16） |
选取经典不锈钢圆管前缘气动加热与传热试验[23]来验证一体化计算方法的准确性，该算例已多次被用于验证高超声速流-热-固多场耦合计算。圆管外径尺寸为Rout =38.1 mm, 内径尺寸为Rin =25.4 mm, 结构材料为不锈钢，表 1给出了不锈钢材料热性能参数[24]。
图 4给出了验证算例的计算网格和边界条件。流场与结构传热计算采用统一网格，图中绿色部分为流场计算网格，红色部分为圆管结构传热计算网格。为了更好捕捉流场激波结构，流场网格在激波位置附近进行加密。圆管的热边界条件为：两端为绝热壁，内壁为等温壁，初始温度为294.4 K。表 2给出了初始来流计算条件，Ma∞、T∞、p∞、Re∞分别为来流马赫数、温度、压力和单位雷诺数。
|
| 图 4 计算网格与边界条件 Fig. 4 Computational grids and boundary conditions |
| Parameter | Ma∞ | T∞/K | p∞ /Pa | Re∞ /m-1 |
| Value | 6.47 | 241.5 | 648.1 | 1.31×106 |
由于近壁面温度梯度的准确计算直接依赖于近壁面网格单元的尺寸，因此为了更加准确计算壁面热流，采用当地网格雷诺数Relocal=ρUΔh/μ[25]来确定物面法向第1层网格的合适高度，其中ρ、U、μ为流体密度、速度和动力黏性系数，Δh为网格高度。文献[25]指出：Relocal应不大于3才能保证壁面热流计算结果的准确性。因此，流场网格在壁面法向第1层高度取1×10-6 m，当地网格雷诺数约为1.21，满足网格计算要求。同时，为保证热流计算的连续性，圆管网格在流固交界面处第1层高度也取1×10-6 m。
2.2 计算结果与分析为了与参考文献和试验结果更好进行对比，采用提出的一体化计算方法对不锈钢材质圆管开展高超声速气动加热与结构传热非定常数值模拟，初始真实时间步长取Δt0 =1×10-4 s。
图 5直观地给出了初始时刻计算得到的圆管外部流场的密度云图与风洞试验纹影图[26]的对比结果，从图中可以看出计算预测的激波位置和流场结构与试验结果吻合较好，说明本算例暂不考虑高温空气化学非平衡效应是合理的。另外，图 6为归一化后的圆管表面压力(p/p0，p0为驻点压力)分布与试验值的对比，从图中可以看出压力沿圆管的分布和试验值吻合较好。
|
| 图 6 圆管表面压力分布对比 Fig. 6 Comparison of surface pressure distributions on cylinder |
图 7给出了初始时刻(t=0 s)圆管外表面归一化热流(q/q0，q0为驻点热流)分布与风洞试验值的对比曲线，结果显示计算预测的热流分布与试验结果保持一致。另外，t=0 s时刻，圆管表面驻点处热流值最大，本文方法计算得到的驻点处热流值为49.21×104 W/m2，与采用Fay-Riddell[27]工程经验公式计算得到的驻点热流值48.27×104 W/m2以及采用黏性激波层[28]计算方法得到的47.02×104 W/m2基本一致。但是与风洞试验值67.0×104 W/m2差别较大，这是由于数值计算中暂未考虑风洞试验中来流湍流度对气动加热热流的影响所导致。
|
| 图 7 t=0 s圆管表面热流分布对比 Fig. 7 Comparison of surface heat flux distributions on cylinder at t=0 s |
图 8展示了不同时刻下流场与圆管结构沿中心线(y=0 m)温度变化曲线。从曲线图可以看出，高速自由来流在圆管前端产生脱体弓形激波，温度从初始的241.5 K急剧跃升到2 162 K，与文献[12]计算结果2 163.7 K基本吻合。另外，本文计算得到的弓形激波的位置约在-54.5 mm处，与激波经验公式[29]计算得到的-54.4 mm基本吻合，说明了本文计算方法可以较准确捕捉高速流场激波位置。初始时刻，在圆管驻点处，流场温度急剧降低至294.4 K，在驻点区域会产生一个明显的薄热边界层，其厚度约为激波层厚度的2.9%。从曲线图可以看出，在薄热边界层中流场温度变化非常明显，产生很大温度梯度，因此在该区域会产生很高的壁面热流。另外，图 8也给出了圆管温度沿中心线随时间变化的放大曲线图，可以观察到随着时间的推移，热量逐渐传入不锈钢圆管内部，圆管温度逐渐升高且趋势明显。
|
| 图 8 不同时刻流场与圆管沿中心线温度分布 Fig. 8 Temperature distribution along centerline within fluid and cylinder domains at different time |
图 9直观地给出了不锈钢圆管内部温度分布随时间变化的云图(t=0.1, 0.5, 1.0, 2.0 s)。可以观察到，随着时间的变化，外部高速气流产生的热量不断传入圆管结构内部，圆管内部温度逐渐升高，高温区域也不断扩大，驻点区域温度始终最高且不断升高。t=2.0 s时刻，本文计算得到驻点温度为389.2 K，与Dechaumphai等[26]的计算结果388.8 K吻合很好，误差在0.1%左右。
|
| 图 9 不同时刻圆管内部温度云图 Fig. 9 Temperature contours within cylinder domain at different time |
为了更加准确对比验证本文计算结果的可靠性，图 10给出了圆管驻点温度(T0)随时间变化曲线与文献[30-31]的对比结果，可以看出本文驻点温度变化趋势与参考文献结果相同，t=2.0 s时，最大温度差值约为3.9 K，计算相对误差在0.9%以内。另外，图中还给出了驻点温度变化率随时间变化曲线，可以观察到初始时刻驻点温度升高速度最快，温升速度逐渐减小趋于平稳。
|
| 图 10 圆管驻点温度随时间的变化 Fig. 10 Temporal variation of stagnation point temperature on cylinder |
图 11为圆管驻点热流与热流变化率随时间的变化曲线。随着时间的推进，圆管壁面温度逐渐升高，驻点区域热边界层内的温度梯度逐渐减小，导致驻点热流逐渐下降。t=2.0 s时，热流下降约6.5%，与Dechaumphai等[26]计算的2 s内驻点热流下降8%接近。同时从虚线可以看出驻点热流下降速度随时间变化也逐渐减小趋于平稳。
|
| 图 11 圆管驻点热流随时间的变化 Fig. 11 Temporal variation of stagnation point heat flux on cylinder |
为了更清楚地验证本文方法计算结果，表 3给出了t=2 s时刻驻点温度和初始时刻驻点热流值与不同参考文献的对比。从表中可以看出本文计算结果与各参考文献值都吻合较好，计算误差在允许范围内。综上，本算例验证了本文流-热-固一体化计算方法的可靠性和可行性，并可用于下文的计算与分析。
另外，为了对比本文一体化计算方法与传统分区耦合计算方法的计算效率，表 4列出了不同计算方法下，迭代步数与计算时间统计结果。本算例在图形工作站(主频2.3 Hz，内存16 GB)上进行串行计算。
| Method | Iterative step | Calculation time/h |
| Partition-coupling calculation | 20 000 | About 20 |
| Integrated calculation (fixed time step) |
20 000 | About 12 |
| Integrated calculation (adaptive time step) |
800 | About 0.8 |
从表 4中可以看出，相比传统分区传递计算方法，一体化计算方法无数据交换计算过程，可以节省计算时间。采用自适应时间步长算法，可以进一步大幅度提高计算效率。
3 “Ⅳ型”激波干扰算例 3.1 计算模型及网格选取与2.1节中尺寸与材质完全相同的不锈钢圆管作为计算模型，开展高超声速“Ⅳ型”激波干扰流-热-固一体化计算分析研究。为了获得高超声速斜激波与正激波干扰产生的“Ⅳ型”激波干扰，本算例计算模型还包含了斜激波发生器，图 12给出了计算模型尺寸和激波发生器的相对位置。表 5给出了初始流场来流条件，U∞为来流速度。
|
| 图 12 计算模型 Fig. 12 Computational model |
| Parameter | Ma∞ | U∞ | T∞/K | p∞/Pa |
| Value | 9.5 | 2 871.28 | 226.5 | 1 197 |
图 13给出了“Ⅳ型”激波干扰算例的计算网格和边界条件。同样，流场与结构传热计算采用统一四边形网格，图中绿色部分为流场计算网格，红色部分为圆管结构传热计算网格。流场网格单元数约176 013，圆管内部结构网格单元约为49 900。流场与结构耦合交界面法向第1层网格高度都取值为1×10-6 m，当地网格雷诺数约为1.98，圆管的热边界条件与2.1节中相同。
|
| 图 13 “Ⅳ型”激波干扰算例的计算网格与边界条件 Fig. 13 Computational grids and boundary conditions of "Type Ⅳ" shock wave interference case |
对不锈钢材质圆管分别开展高超声速“Ⅳ型”激波干扰气动加热与结构传热一体化的定常与非定常数值计算，给出计算结果并进行分析。
3.2.1 定常计算结果在进行气动加热与结构传热一体化分析之前，首先要了解高超声速“Ⅳ型”激波干扰的流动特征。图 14显示了圆管前缘的“Ⅳ型”激波干扰流动的稳态压力云图与流线结构。可以观察到，流场具有明显的弓形激波和透射激波的特征。计算发现当入射斜激波与圆管前方的弓形激波接近正激波的位置相交时，发生“Ⅳ型”激波干扰现象。激波发生器产生的斜激波与圆管前端产生的弓形激波干扰产生一道明显的透射激波，将弓形激波分为上下两个部分，且激波干扰产生的高温高压使得弓形脱体激波的脱体距离明显增大。同时，上、下两道弓形激波波后的亚声速区中产生上、下两道剪切层，在这两层剪切层之间产生一股超声速“喷流”，超声速“喷流”冲击圆管壁面产生一道较短的弓形滞止激波，形成一个小的滞止加热区域，引起壁面热流、压力急剧上升。由此看出，“Ⅳ型”激波干扰流场结构是非常复杂的。
|
| 图 14 “Ⅳ型”激波干扰稳态流场结构 Fig. 14 Steady flow field structure of "Type Ⅳ" shock wave interference |
图 15展示了“Ⅳ型”激波干扰流场与圆管结构内部的稳态温度分布。自由来流经过弓形激波第1次加热后，超声速“喷流”产生的二次弓形滞止激波对其进一步加热，流场温度从初始的226.5 K急剧增加，最高温度达到3 722.1 K。不锈钢圆管由于外部流场的持续气动加热，内部结构温度整体升高。其中，二次弓形激波在圆管壁面形成了一个局部的高温高压滞止加热区域，由此导致圆管结构在该区域形成一个高温区，最高温度可达3 406.2 K。
|
| 图 15 流场与圆管稳态温度云图 Fig. 15 Steady temperature contours within fluid and cylinder domains |
为了定量分析“Ⅳ型”激波干扰产生的高温高压效应，图 16给出了圆管表面稳态热流与压力系数(Cp)分布曲线。可以明显看出，“Ⅳ型”激波干扰产生的超声速“喷流”冲击圆管壁面导致表面热流与压力系数出现显著峰值。超声速“喷流”冲击壁面的作用位置大致位于120°的地方，在此处，表面热流与压力系数分别高达4 348.1 kW/m2和15.8。
|
| 图 16 圆管表面稳态热流与压力系数分布曲线 Fig. 16 Curves of steady surface heat flux and pressure coefficient distributions on cylinder |
值得一提的是，本定常算例也展示了本文一体化计算方法的优势：可以高效求解复杂高超声速流场下的流-热-固耦合定常计算问题。传统分区耦合计算方法在求解定常问题时，需要在时间域上进行大量的耦合交替迭代直至流场与固体结构温度场定常收敛，因此必将导致计算效率的大幅度降低，相比之下，本文一体化方法可以很好解决这一问题。
3.2.2 非定常计算结果图 17给出了非定常计算中不锈钢圆管内部温度分布随时间变化的云图(t=0.1, 0.5, 1.0, 2.0 s)。从图中可以直观地观察内部结构的温度分布随时间的演变历程。由于“Ⅳ型”激波干扰的作用，产生的超声速“喷流”不断冲击壁面，产生的高热量首先在冲击区域快速积累，导致该区域的圆管内部结构温度急剧上升，形成一个显著的高温区。随着时间的推移，“喷流”冲击产生的高热流持续传入圆管结构内部且不断延伸，导致结构内部整体温度不断升高，同时冲击区域的高温区域也不断扩大。由此可见，“Ⅳ型”激波干扰的热冲击作用导致结构内部产生集中的高温区，在此高温区域必然会出现较大的热应力集中现象，这也将是后续研究的重要内容。
|
| 图 17 圆管内部温度云图随时间变化历程 Fig. 17 Temporal variation of temperature contours within cylinder |
图 18和图 19分别给出了圆管表面的压力系数分布和热流分布随时间的变化历程。同时，为了更好地对比分析“Ⅳ型”激波干扰产生的高温高压效应，图中也给出了无斜激波干扰时，初始时刻圆管表面的压力与热流分布，并用此时的驻点压力系数(Cp0)和热流值对预测的“Ⅳ型”激波干扰壁面压力系数和热通量进行统一归一化处理。
|
| 图 18 圆管表面压力系数随时间变化 Fig. 18 Pressure coefficient variation with time along cylinder surface |
|
| 图 19 圆管表面热流随时间变化 Fig. 19 Heat flux variation with time along cylinder surface |
从图 18和图 19可明显观察到，“Ⅳ型”激波干扰的确导致圆管表面压力和热流急剧增大，且峰值都出现在超声速“喷流”冲击壁面区域。相比于无斜激波干扰的流动，“Ⅳ型”激波干扰的射流冲击作用使得壁面最大压力系数增大约9倍，壁面最大热流增大约4.7倍，计算结果与文献[32-33]吻合。由此可见，此类“Ⅳ型”激波干扰将极大放大外部高速流场对飞行器结构的气动热与气动力载荷作用，给高速飞行器飞行安全造成极大安全隐患。另外，从图中可以看出，随着时间的推移，圆管表面压力分布几乎不随时间变化，这说明圆管壁面压力受圆管整体结构温度变化的影响很小。然而，壁面热通量的峰值随着计算时间的推进逐渐减小。这是由于初始时刻壁面热流峰值最大，圆管结构吸收气动加热热量后，壁面温度快速上升。壁温的升高导致温度边界层内的温度梯度减小，壁面热流峰值减小，削弱了外部气动加热效应从而导致圆管结构温升速率的降低。同时，结构温升速率的降低反过来也会导致壁面热通量的变化率逐渐减小。这些现象正是反映了流-热-固相互作用的耦合特性。
综上分析，高速飞行器在进行长航时飞行时，“Ⅳ型”激波干扰作用产生的高温高压冲击载荷，极有可能造成飞行器热防护结构的热力学破坏从而形成严重的飞行安全隐患。因此，“Ⅳ型”激波干扰作用将给高速飞行器的热防护设计与选材带来严峻挑战，需引起重要关注。
4 结论本文针对高超声速进气道前缘“Ⅳ型”激波干扰产生的气动加热与结构传热耦合计算问题，提出一种基于有限体积法的流-热-固一体化计算方法，并采用经典高超声速二维圆管绕流算例验证了一体化方法的正确性与稳定性。采用所提出的一体化计算方法对高超声速前缘“Ⅳ型”激波干扰中的流-热-固耦合问题进行一体化定常/非定常计算分析，得到以下结论：
1) 采用经典高超声速二维圆管绕流非定常算例对本文发展的一体化计算方法进行验证分析，计算结果与参考文献和风洞试验结果吻合较好，证明了本文高超声速流-热-固一体化计算方法的可行性与可靠性，同时也证明了自适应时间步长控制方法的正确性。
2) 本文所提出的高超声速流-热-固一体化求解方法可以高效解决流场与结构传热稳态求解问题，较快计算出稳态结构与流场的温度分布。计算方法进行全物理场统一迭代，可以很好地改善传统分区耦合算法多次迭代计算效率低的不足。
3) 对高超声速前缘“Ⅳ型”激波干扰流-热-固耦合进行定常/非定常一体化计算分析，计算发现，“Ⅳ型”激波干扰作用产生的超声速“喷流”不断冲击壁面，使得壁面最大压力系数增大约9倍，壁面最大热流增大约4.7倍，给高速飞行器的热防护设计与选材带来严峻挑战。
为进一步探索高超声速流-热-固一体化求解方法的应用范围，后续将研究“Ⅳ型”激波干扰作用产生的高温高压对结构热力学变形产生的影响。
| [1] |
王江峰, 伍贻兆, 季卫栋, 等. 高超声速复杂气动问题数值方法研究进展[J]. 航空学报, 2015, 36(1): 159-175. WANG J F, WU Y Z, JI W D, et al. Progress in numerical simulation techniques of hypersonic aerodynamic problems[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(1): 159-175. (in Chinese) |
| Cited By in Cnki (7) | Click to display the text | |
| [2] | ALBERTSON C, VENKAT V. Shock interaction control for scramjet cowl leading edges:AIAA-2005-3289[J]. Reston, VA:AIAA, 2005. |
| [3] | EDNEY B. Anomalous heat transfer and pressure distributions on blunt bodies at hypersonic speeds in the presence of an impinging shock: Technical Report 115[R]. Stockholm: The Aeronautical Research Institute of Sweden, 1968. |
| [4] | LU H B, YUE L J, XIAO Y B, et al. Interaction of isentropic compression waves with a bow shock[J]. AIAA Journal, 2013, 51(10): 2474-2484. |
| Click to display the text | |
| [5] |
肖丰收.若干典型高超声速激波干扰流动特性研究[D].合肥: 中国科学技术大学, 2016: 21-32. XIAO F S. Research on flow characteristics of some typical hypersonic shock wave interactions[D]. Hefei: University of Science and Technology of China, 2016: 21-32(in Chinese). |
| Cited By in Cnki | Click to display the text | |
| [6] |
桂业伟, 刘磊, 代光月, 等. 高超声速飞行器流-热-固耦合研究现状与软件开发[J]. 航空学报, 2017, 38(7): 020844. GUI Y W, LIU L, DAI G Y, et al. Research status of hypersonic vehicle fluid-thermal-solid coupling and software development[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(7): 020844. (in Chinese) |
| Cited By in Cnki (3) | Click to display the text | |
| [7] | THORNTON E A, DECHAUMPHAI P. Coupled flow, thermal, and structural analysis of aerodynamically heated panels[J]. Journal of Aircraft, 1988, 25(11): 1052-1059. |
| Click to display the text | |
| [8] | KAZEMI-KAMYAB V, VAN ZUIJLEN A H, BIJL H. Analysis and application of high order implicit Runge-Kutta schemes for unsteady conjugate heat transfer:A strongly-coupled approach[J]. Journal of Computational Physics, 2014, 272: 471-486. |
| Click to display the text | |
| [9] | CHEN F, LIU H, ZHANG S T. Coupled heat transfer and thermo-mechanical behavior of hypersonic cylindrical leading edges[J]. International Journal of Heat and Mass Transfer, 2018, 122: 846-862. |
| Click to display the text | |
| [10] | QIN Q H, XU J L, GUO S. Fluid-thermal analysis of aerodynamic heating over spiked blunt body configurations[J]. Acta Astronautica, 2017, 132: 230-242. |
| Click to display the text | |
| [11] | ZHAO X L, SUN Z X, TANG L S, et al. Coupled flow-thermal-structural analysis of hypersonic aerodynamically heated cylindrical leading edge[J]. Engineering Applications of Computational Fluid Mechanics, 2011, 5(2): 170-179. |
| Click to display the text | |
| [12] | ZHANG S T, CHEN F, LIU H. Time-adaptive, loosely coupled strategy for conjugate heat transfer problems in hypersonic flows[J]. Journal of Thermophysics and Heat Transfer, 2014, 28(4): 1-12. |
| Click to display the text | |
| [13] | WIETING A R, DECHAUMPHAI P, BEY K S, et al. Application of integrated fluid-thermal-structural analysis methods[J]. Thin-Walled Structures, 1991, 11(1-2): 1-23. |
| Click to display the text | |
| [14] |
李鹏飞, 吴颂平. 类航天飞机前身结构与高超声速流场的耦合传热模拟分析[J]. 航空动力学报, 2010, 25(8): 1705-1710. LI P F, WU S P. Numerical simulation of fluid-solid-thermal interaction in hypersonic flows[J]. Journal of Aerospace Power, 2010, 25(8): 1705-1710. (in Chinese) |
| Cited By in Cnki | Click to display the text | |
| [15] |
姜贵庆, 童秉纲, 曹树声. 以有限元方法为主体的计算气动热力学[J]. 力学与实践, 1992, 14(3): 1-8. JIANG G Q, TONG B G, CAO S S. Computational aerothermodynamics based on finite element method[J]. Mechanics and Engineering, 1992, 14(3): 1-8. (in Chinese) |
| Cited By in Cnki (5) | Click to display the text | |
| [16] | LIOU M S. A sequel to AUSM:AUSM+[J]. Journal of Computational Physics, 1996, 129(2): 364-382. |
| Click to display the text | |
| [17] | VENKATAKRISHNAN V. On the accuracy of limiters and convergence to steady state solutions: AIAA-1993-0880[R]. Reston, VA: AIAA, 1993. |
| [18] | VAN LEER B. Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[J]. Journal of Computational Physics, 1979, 32(1): 101-136. |
| Click to display the text | |
| [19] | MENTER F R. Two-equation eddy-viscosity turbulence models for engineering applications[J]. AIAA Journal, 1994, 32(8): 1598-1605. |
| Click to display the text | |
| [20] | HASELBACHER A, BLAZEK J. Accurate and efficient discretization of Navier-Stokes equations on mixed grids[J]. AIAA Journal, 2000, 38(11): 2094-2102. |
| Click to display the text | |
| [21] |
李佳伟, 王江峰, 杨天鹏, 等. 钝体外形气动加热与结构传热一体化数值模拟[J]. 推进技术, 2019, 40(1): 33-43. LI J W, WANG J F, YANG T P, et al. Numerical simulation of integrated aeroheating-structural heat transfer study for blunt body[J]. Journal of Propulsion Technology, 2019, 40(1): 33-43. (in Chinese) |
| Cited By in Cnki | Click to display the text | |
| [22] | VALLI A M P, CAREY G F, COUTINHO A L G A. Control strategies for timestep selection in finite element simulation of incompressible flows and coupled reaction-convection-diffusion processes[J]. International Journal for Numerical Methods in Fluids, 2005, 47(3): 201-231. |
| Click to display the text | |
| [23] | WIETING A R. Experimental study of shock wave interference heating on a cylindrical edge: NASA-TM-100484[R]. Washington, D.C.: NASA, 1987. |
| [24] |
李佳伟, 王江峰, 杨天鹏, 等. 高超声速飞行器前缘流-热-固一体化计算[J]. 国防科技大学学报, 2018, 40(6): 9-16. LI J W, WANG J F, YANG T P, et al. Fluid-thermal-structural study of integrated algorithm for aerodynamically hypersonic heated leading edges[J]. Journal of National University of Defense Technology, 2018, 40(6): 9-16. (in Chinese) |
| Cited By in Cnki | Click to display the text | |
| [25] | PAPADOPOULOS P, VENKATAPATHY E, PRABHU D, et al. Current grid-generation strategies and future requirements in hypersonic vehicle design, analysis and testing[J]. Applied Mathematical Modelling, 1999, 23(9): 705-735. |
| Click to display the text | |
| [26] | DECHAUMPHAI P, THORNTON E A, WIETING A R. Flow-thermal-structural study of aerodynamically heated leading edges[J]. Journal of Spacecraft and Rockets, 1989, 26(4): 201-209. |
| Click to display the text | |
| [27] | FAY J A, RIDDELL F R. Theory of stagnation point heat transfer in dissociated air[J]. Journal of the Aeronautical Sciences, 1958, 25(2): 73-85. |
| Click to display the text | |
| [28] | HOLCOMB J E, CURTIS J T, SHOPE F L. A new version of the CVEQ hemisphere viscous shock layer program for equilibrium air: NASA TN AEDC-TMR-85-V7[R]. Washington, D.C.: NASA, 1985. |
| [29] | BILLIG F S. Shock-wave shapes around spherical and cylindrical-nosed bodies[J]. Journal of Spacecraft and Rockets, 1967, 4(6): 822-823. |
| Click to display the text | |
| [30] |
耿湘人, 张涵信, 沈清. 高超飞行器流场和固体结构温度场一体化计算新方法的初步研究[J]. 空气动力学学报, 2002, 20(4): 422-427. GENG X R, ZHANG H X, SHEN Q. Study on an integrated algorithm for the flow fields of high-speed vehicles and the heat transfer in solid structures[J]. Acta Aerodynamica Sinica, 2002, 20(4): 422-427. (in Chinese) |
| Cited By in Cnki (46) | Click to display the text | |
| [31] |
黄杰.高超声速飞行器流热固多物理场耦合计算研究[D].哈尔滨: 哈尔滨工业大学, 2013: 33-36. HUANG J. Study on hypersonic vehicle fluid-thermal-structure muti-physics coupling calculation[D]. Harbin: Harbin Institute of Technology, 2013: 33-36(in Chinese). |
| Cited By in Cnki (13) | Click to display the text | |
| [32] | WIETING A R, HOLDEN M S. Experimental study of shock wave interference heating on a cylindrical leading edge at Mach 6 and 8: AIAA-1987-1511[R]. Reston, VA: AIAA, 1987. |
| [33] |
张胜涛, 陈方, 刘洪. 高超声速进气道前缘流场-热-结构耦合分析[J]. 空气动力学学报, 2017, 35(3): 436-443. ZHANG S T, CHEN F, LIU H. Fluid-thermal-structural coupling analysis on leading edge of hypersonic inlets[J]. Acta Aerodynamica Sinica, 2017, 35(3): 436-443. (in Chinese) |
| Cited By in Cnki | Click to display the text |