0 引言

当高超声速飞行器在近空间大气层内进行远程机动飞行时, 将面临具有复杂流场、高焓、中低热流和持续长时间气动加热等特征的气动热力学环境。尽管面临的飞行环境如此恶劣, 但是高超声速飞行器的热防护设计要求却极为苛刻, 即在满足重量尽可能轻的约束条件下实现高效、可重复使用的非烧蚀热防护。这些苛刻要求对热防护系统综合性能优化提出了严峻挑战。为了实现更合理的结构选材和防热优化, 最大限度减小热结构设计冗余, 必须更准确地预测飞行器的热力环境特性。由于高超声速飞行器采用机身/发动机高度一体化设计, 其在近空间大气层内进行持续长时间、宽马赫数域的远程机动飞行过程中，外部流场的气动力/气动热和内部结构的传热/应力/变形等各物理场间的相互作用将变得极为强烈。因此，要准确地预测飞行器的热力环境特性，必须统筹考虑流场、热和结构等各物理场间的相互耦合关系，建立多场耦合一体化分析方法，从而为高超声速飞行器热防护设计提供有效的理论与技术支撑。

国外学者Thornton, Dechaumphai和Wieting等人^[1-3]最早提出了流场、热和结构之间多场耦合问题，并基于双向耦合方法发展了流场-热-结构一体化有限元方法，为分析更为真实复杂的多场耦合问题奠定了基础。其后，相继发展了迭代耦合方法^[4-5]和基于分区耦合方法的交错策略^[6]以及将非稳态问题看成由一系列“准稳态问题”的混合耦合策略^[7]和统一求解耦合策略^[8]。为了进行高时间精度的流场-热-结构耦合分析，Miller等人^[9-10]采用多循环方法发展了基于松耦合的多物理时间推进策略。国内学者黄唐和毛国良等人^[11]率先研究了基于松耦合策略的流场、热和结构一体化数值模拟方法。随着对于持续长时间气动热力环境下流场-热-结构耦合问题重要性的认识，陆续开展了流场、热和结构多场耦合计算研究和分析^[12-19]。随着高超声速飞行器研制的不断发展，多场耦合分析的思想和应用逐渐深入到工程实际，但是对于复杂热力环境作用下的多物理场耦合特征和规律的认识还不够深入，相关耦合性分析方法有待进一步发展和完善。

本文首先通过对高超声速流场-热-结构耦合问题的机理过程分析，提出了基于稳态飞行的松耦合分析策略，然后在此基础上发展了高超声速流场-热-结构耦合分析平台HyCCD (Hypersonic Coupling Computational Dynamics)，最后针对激波相互作用下的高超声速进气道前缘结构进行了流场-热-结构耦合计算分析研究。

1 耦合分析策略 1.1 流场-热-结构耦合模型

高超声速流场-热-结构耦合问题的耦合机理是流体介质内的流场气动力/气动热等物理属性和行为与固体介质内的结构传热/应力/变形等物理属性和行为通过流-固耦合界面发生相互作用和影响的物理过程。因此耦合策略的制定既需考虑流体介质与固体介质之间特征参量的面耦合，也需考虑它们各自内部特征参量的体耦合，为此本文考虑如图 1所示的高超声速流场-热-结构耦合模型。

在流体介质内部，气动力/气动热之间的耦合由统一流动控制方程描述。利用计算流体力学(CFD)方法求解流动控制方程，获得流体介质内部的气动力/气动热参量，通过流-固耦合界面向固体介质提供热载荷(壁面热流q)和力载荷(壁面压力p)。在固体介质内部，结构传热/应力/变形之间的耦合由热传导控制方程和热弹性力学控制方程描述。考虑结构温度和变形的耦合效应，利用基于有限单元法(FEM)的计算热-结构动力学(CTD/CSD，或CTSD)统一求解热传导控制方程和热弹性力学控制方程，获得固体介质内部的热/力耦合响应参量，通过流-固耦合界面向流体介质提供温度条件(壁面温度T)和结构变形条件(壁面位移u)。

1.2 基于稳态飞行的松耦合分析策略

基于上述流场-热-结构耦合模型，本文建立了如图 2所示的基于稳态飞行的松耦合分析策略。图中，Δt_F为流场计算时间步长；Δt_TS为固体热-结构体耦合计算时间步长；Δt_C为流-固面耦合计算时间步长，可取为若干倍的固体热-结构体耦合计算时间步长，即Δt_C=nΔt_TS，(n=1, 2, 3, …)。

根据耦合物理过程的时间演化自适应地选取流-固面耦合计算时间步长Δt_C，对于持续长时间耦合计算分析效率具有重要意义。本文将耦合计算时间步长Δt_C的选取看作为一个自动控制问题，采用经典的PID (proportional-integral-derivative)控制器对耦合计算时间步长进行自适应控制^[20]，则有：

(1)

式中：下标和上标“n”表示耦合时间步；TOL为人为给定的输入容差；K_P、K_I和K_D为PID控制器的经验参数。本文采用文献[21]建议的数值，即K_P=0.075，K_I=0.175，K_D=0.01；r_n为耦合界面温度在t_n时刻的变化度量：

(2)

为了防止耦合计算时间步长的过大或过小增长，引入如下步长限制为：

(3)

(4)

式中：Δt_min和Δt_max分别为人为给定的最小和最大耦合计算时间步长；m和M分别为人为给定的最小和最大步长增长率。

耦合策略的具体实施流程可以概述为：

1) 在初始t₀时刻，首先给定固体域某一初始恒温或温度场分布以及初始受力和位移约束，然后通过界面信息传递方法将固体域的壁面温度和位移量传递给流体域，其中壁面温度作为流场计算的边界条件，而位移量用来更新流场网格；

2) 基于施加的温度边界条件和更新的流场网格，在流体域内进行稳态流场计算，最终得到稳态流场的壁面热流密度和壁面压力；

3) 通过界面信息传递方法将稳态流场的壁面热流密度和壁面压力传递给固体域，分别作为固体热-结构动力学计算的热载荷和力载荷；

4) 基于施加的热载荷和力载荷，在固体域内进行瞬态热-结构动力学计算，时间推进从t₀时刻至t₀+Δt_C时刻，最终得到t₀+Δt_C时刻的固体结构热/力耦合响应参量；

5) 通过界面信息传递方法将t₀+Δt_C时刻固体域的壁面温度和位移量传递给流体域，其中壁面温度作为流场计算的边界条件，而位移量用来更新流场网格；

6) 至此完成一个耦合计算时间步长Δt_C的计算，然后进入下一个耦合计算时间步长继续进行计算，以此类推直至整个耦合时间结束。

2 控制方程及数值方法 2.1 流动控制方程求解

考虑包含NS种组分的化学反应完全气体混合物，则可压缩粘性化学非平衡流动的Navier-Stokes控制方程组在笛卡尔直角坐标系下的守恒积分形式写为^[22-23]：

(5)

其中，

(6)

分别为流场守恒变量、对流通量、粘性通量和化学反应源项。式中：ρ为混合气体的密度，ρ=∑ρ_s，其中ρ_s为组分s的密度；Y_s为组分s的质量分数，Y_s=ρ_s/ρ；u_i为混合气体在坐标x_i方向上的速度分量；p为混合气体的热力学压力；E和H分别为混合气体的单位质量总能和单位质量总焓；J_si为组分s的质量扩散通量；q_i为混合气体的热流密度矢量；σ_ij为混合气体的粘性应力张量；为组分s的单位体积质量净生成率。

基于有限体积法对上述流动控制方程进行离散求解。其中，对流通量采用M-AUSMPW+混合迎风格式，并用多维限制器MLP提高至三阶精度；粘性通量采用二阶中心差分格式进行离散；时间推进采用LU-SGS隐式方法。为了加速收敛过程，采用当地时间步长和隐式残值平均等加速收敛措施。

2.2 热弹性力学控制方程求解

考虑结构温度和变形的相互耦合效应，利用有限单元法(FEM)进行统一耦合求解。将变分原理应用于热传导控制方程和热弹性力学控制方程以及它们相应的边界条件，最终可得到如下形式的有限单元矩阵方程为：

(7)

式中：和分别为节点位移列阵、节点位移对时间的一阶导数列阵和二阶导数列阵；{T}、和分别为节点温度列阵、节点温度对时间的一阶导数列阵和二阶导数列阵；{F}为力载荷列阵；{Q}为热载荷列阵；[M]为结构质量矩阵；[C]为结构阻尼矩阵；[K]为结构刚度矩阵；[C_T]为热容矩阵；[K_T]为热传导矩阵；[C_Tu]和[K_uT]分别为热弹性阻尼矩阵和热弹性刚度矩阵，它们反映了温度与变形之间的耦合效应。

利用ANSYS Mechanical APDL完成热-结构耦合计算，采用的单元类型为SOLID226。该单元具有20个节点，每个节点具有五个自由度。对于热-结构耦合分析，节点自由度主要是三个方向的位移分量和温度。

2.3 界面信息传递方法

综合考虑反距离加权插值法^[24]和径向基函数插值法^[25]的特点和耦合变量(位移、温度、压力和热流)的物理特性，提出如图 3所示的混合插值策略。壁面温度和位移量的传递需要满足其分布特性，采用反距离加权插值法；壁面热流和壁面压力的传递需要保证其插值前后的守恒特性，采用径向基函数插值法。同时考虑到在插值过程中线性方程组可能会出现求解失败，则采用反距离加权插值法以保证耦合计算的顺利进行。

2.4 网格变形方法

考虑固体结构的小变形问题，基于反距离加权插值^{[24, 26]}实现复杂外形网格变形方法。在已知网格边界节点位移u_i^b后，则网格内部节点的位移u_jⁱⁿ可由反距离加权插值法插值得到：

(8)

式中：NB为网格边界节点的个数；w_i为距离权函数，其表达式为：

(9)

式中：d_ji为网格内部节点r_j与网格边界节点r_i之间的欧氏距离，d_ji=‖r_j－r_i‖，r为位置矢量；u为距离的幂次，取u=2。当所有网格节点的位移u_j已知后，则可以更新整个流场网格：

(10)

式中：x_j^old为原网格节点的坐标；x_j^new为更新后的网格节点坐标。

3 计算结果及分析 3.1 计算模型

高超声速飞行器前体预压缩斜激波与进气道前缘激波可能相交，产生复杂的激波相互作用，导致进气道前缘的气动加热问题变得更为严重。如图 4所示，Case 1为超临界状态，斜激波入射在进气道唇口以内，进气道前缘处于远场自由来流条件；Case 2为临界状态，斜激波正好入射在进气道唇口上，与进气道前缘激波发生相互作用；Case 3为亚临界状态，斜激波入射在进气道唇口以外，进气道前缘处于斜激波波后的自由来流条件。本文以进气道前缘为计算模型，考虑Case 2工况，它包含Case 1和Case 3两种流动条件(见表 1)，持续耦合计算时间为11s。

外部流场计算采用化学非平衡气体模型，壁面为非催化条件；固体结构初始温度为300K，初始应力为0，热应力参考温度为300K。如图 5所示，考虑外部流场的气动力/气动热载荷作用，同时考虑壁面的辐射效应，壁面发射率为ε=1.0；内壁面施加均匀分布的压力载荷p_∞=1197 Pa；模型两端位置约束条件采用固支约束。进气道前缘的几何模型与文献[19]中的圆柱前缘模型相一致，CFD网格采用结构化划分，大小为201×111×11，在激波和近壁面附近进行加密处理；FEM网格采用非结构化划分，最小单元边长为0.5mm。前缘材料采用钛合金，其物理属性见表 2。

3.2 初始流场

图 6所示为该进气道前缘模型在初始t=0s时刻的绕流流场结构云图。可以看出，其流场结构大致可以划分为远场自由来流区及其产生的弓形激波区、入射激波波后自由来流区及其产生的弓形激波区。在这两个弓形激波波后区之间产生由两剪切层形成的“流管”，导致壁面压力和热流密度出现很高的极值。如图 7所示，相比于未受激波干扰的流动，激波相互作用使壁面压力系数增大约8倍，而使壁面热流密度增大约5倍。由此可见，这种激波相互作用现象的存在将极大地放大外部流场的气动力/气动热载荷作用，从而给飞行器热防护设计带来更为严峻的挑战。

3.3 流场-热-结构耦合特性分析

图 8所示为进气道前缘模型的壁面温度分布和变形量随时间的演化历程。由图 8(a)可以看出，由于激波相互作用而导致的高热流犹如一道高强度能量束撞击在壁面上，引起壁面被撞击处的温度随时间迅速升高。在短短11s时间内，考虑辐射效应的壁面温度最高点已达至初始温度的近8倍，远高于壁面驻点处的温度及其温升速度。同时也可以看到，在t=11s时刻的壁面被撞击处的温升速度已经明显减弱，随着时间的再持续，该位置将率先达到热平衡状态，而其他部位的壁面温度会继续升高，直至达到热平衡状态。从图 8(b)中可以看出，在激波相互作用导致的高热流作用区域，壁面变形随时间变化较快，变形量也最大；随着时间的持续，壁面变形量最大点逐渐由激波相互作用的部位向驻点偏移。

图 9所示为进气道前缘模型的壁面热流密度分布和压力系数分布随时间的演化历程。壁面热流密度随壁面温度的逐渐升高而减弱，体现了流-固之间传热的耦合性；随着时间的不断持续，最终将会达到热平衡状态；而壁面压力受壁面温度变化的影响不大，随着时间的持续，壁面变形量逐渐增大，在壁面变形量最大区域，压力开始有所增大。

图 10所示为进气道前缘内部结构温度分布和应力分布随时间的演化历程。从图 10中可以直观地看出，在高热流撞击作用下，热量首先在撞击区域附近的结构内部迅速累积，致使该处的结构温度迅速升高。随着时间的持续，热量不断累积，逐渐向结构内部纵深区域传递；结构温度也在不断升高，逐渐向结构内部纵深区域不断延展。而且，激波相互作用撞击区域附近的结构内部出现应力集中，随着时间的持续也逐渐向结构内部纵深区域延展。

综上分析，由于激波相互作用，进气道前缘将遭受强烈的气动力/气动热载荷冲击，在实际飞行过程中即使作用时间很短，也足以造成热防护结构的热学和力学破坏，使其热防护选材和结构设计面临严峻挑战。而且，高超声速飞行器绕流流场普遍存在激波相互作用现象，这些部位将会更长时间地遭受激波相互作用的冲击，因而其热防护也同样会变得更为严峻。

4 结论

在高超声速流场-热-结构耦合分析平台HyCCD基础上，针对高超声速飞行器进气道前缘进行了初步分析。从中可以发现：激波相互作用使壁面压力系数增大约8倍，而使壁面热流密度增大约5倍。随着时间的持续，流-固之间传热的耦合性使壁面最终达到热平衡状态，而热量逐渐向结构内部纵深区域传递；壁面变形量逐渐增大，在壁面变形量最大区域，压力开始有所增大。在高热流撞击作用下，撞击区域附近的结构温度迅速升高，集中应力逐渐向结构内部纵深区域不断延展。这些规律的认识和分析表明本文发展的HyCCD平台有效地预测了持续长时间飞行条件下的热/力耦合响应时空分布特征和规律，可以用来分析更为真实复杂结构的高超声速流场-热-结构耦合问题，从而为热防护结构的热学/力学综合性能评估以及进一步完善防热结构设计提供理论依据和分析工具。