舵翼结构模型气动热力耦合关系仿真研究

近空间飞行器拥有快速响应、突防能力强等特点，能够在短时间内精准打击长距离军事目标，在国防科研、空间技术研究中具有极其重要的战略意义[1-2]。飞行器在高马赫数飞行过程中，周围空气受到激波的强烈压缩，对飞行器表面产生剧烈的气动热力载荷。舵翼结构作为飞行器关键部件，为了确保燃料供给能够满足高速、长距离飞行，要求尽量减小结构尺寸，结构刚度的降低导致更大的变形量，影响其气动外形，进一步增强气动热、气动力与结构响应之间的耦合效应[3-5]，直接影响飞行稳定性，严重时造成偏离飞行轨道甚至掉落，影响精确打击效果，这也是目前近空间飞行器热防护设计与结构设计过程中必然要面对的问题。

为此，多位研究者开展了相关的试验项目。夏吝时等[6]对主动式热疏导简化前缘金属开展地面热环境模拟试验。吴大方等[7-8]对近空间飞行器复合材料先后开展了1 100 ℃和1 500 ℃的热振联合实验。蔡礼港、欧朝、李强等[9-11]在电弧风洞、激波风洞中对飞行器结构与热防护材料开展了实验研究。然而，开展试验周期长、难度大、费用高，使用数值仿真方法对近空间飞行器结构温度、变形进行预测，能够更为高效地为热防护设计与结构设计提供指导。李佳伟等[12]采用一体化求解方法，将流场与结构温度场作为一个物理场，采用统一的控制方程进行求解，避开了传统耦合方法在迭代计算时场之间的大量数据交换。季卫栋等[13]采用统一的积分方程组作为气动加热和结构热物理过程的控制方程，对物理场进行统一的迎风格式有限体积方法离散，对二维圆管模型的气动加热和结构传热问题进行了模拟。陈梅洁等[14]通过顺序耦合方法，基于Transition k-k1-ω湍流模型计算得到了温度场和应力场的分布。梁杰等[15]采用蒙特卡洛方法，通过分区求解策略对大型航天器离轨再入陨落过程中的气动力和气动热特性进行了数值模拟。

当前研究者更多关注于多物理场耦合方法的实现，而在实现仿真的基础上，开展各物理场之间的耦合关系分析，对于热防护材料设计及其结构设计有着更为现实的意义[16]。为此，本文通过分区耦合方法，开展舵翼结构流-热-固耦合仿真，对不同飞行迎角状态下结构变形场、温度场进行计算，实现气动力、气动热及结构耦合响应的量化分析。

1 多场耦合理论

1.1 数值计算理论

流-热-固多场耦合根据以下6个方程进行数值计算。

1) 质量守恒方程。

式(1)中：▽代表拉普拉斯算子；V为速度矢量(u、v、w)，指流体微元在x、y、z等3个方向的速度分量(m/s)； ρf为流体密度(kg/m3)； t为时间(s)。

2) 动量守恒方程。

式(2)中：p指作用在所划分的流体微元上的压强(Pa)； fx、 fy、 fz分别指的是流体微元体积力3个方向的分量(N)；

是粘性力在3个方向的分量。

3) 能量守恒方程。

式(3)中：E为总能量(J)； ρ为材料密度(kg/m3)； fw= fxu+ fyv + fzw为单位质量的流体体积力所做的功(J)。

4) 气体状态方程。

式(4)中：M表示气体摩尔质量(kg/mol)； R*为摩尔气体常数，R*=8.314 5 J/(mol·K)； T为温度场(K)； Z为压缩因子。

式(5)中： λ为材料的导热系数(W/(m·K))； C为比热容(J/kg·K)； Q为物体内部的热源强度(W/kg)。

6) 结构响应方程。

式(6)中，[σ]为应力矩阵；[D]为结构材料弹性矩阵；[ε]为结构的总应变矩阵。

1.2 多物理场耦合方法

舵翼结构在近空间环境中飞行时，气动加热过程是一个持续的、非稳态的过程，并且是流场、温度场和结构场相互耦合作用的过程，因此采用流-热-固耦合方法进行多物理场耦合问题的求解。热流固耦合模型如图1所示，在流体域内使用基于有限体积法的Fluent求解器，对流场的热流和压力进行求解，使用System Coupling模块，通过流固耦合界面将热载与气动力载荷传递至固体域。在固体域内使用Transient Structural求解器，对固体域的温度与变形进行求解，使用System Coupling模块通过流固耦合界面将温度场与变形场传递至流体域，继续进行流体域的下一步计算，直到达到设置的耦合计算时间。

2 仿真算法校验

为验证本文选用多物理场耦合算法的准确性，利用经典的圆管绕流试验作为验证算例。该试验是1987年美国科学家 Allan R.Wieting在NASA LANGLEY 8-foot高温风洞中进行的，试验如图2所示。

风洞试验中，不锈钢圆管是完全对称的圆柱体，因此本文中采用1/2模型进行计算，以提高计算效率。利用SolidWorks软件建立流场及圆柱模型，采用结构化网格，在边界层处和激波位置处进行网格加密处理。流体域和固体域在耦合面上网格节点是连续的，保证流固耦合面上的几何相容性，模型网格如图3所示。流体域网格数量为22 866，固体域网格数量为3 249，网格质量均大于0.985。流体域入口、出口设置为远场边界条件，流场参数设置如表1所示，耦合面按流固耦合边界条件设置。圆管内径为25.4 mm，外径为38.1 mm，材料为1Cr18Ni9Ti不锈钢，材料参数如表2所示，圆管初始温度为294.4 K。

仿真耦合求解时间设置为2s，与文献[17]相同。图4为仿真流场温度分布云图与文献试验激波纹影，二者形状相吻合。图5表示了仿真结果与文献中的流场温度突变位置、流场温度，二者分别相差0.071%与2.369%。图6表示了仿真结果与文献中的圆管结构温度分布、应力分布，二者分别相差0.345%和0.567%，具体数据如表3所示。由上述分析可知，仿真结果与文献所述基本吻合，说明本文所使用的多场耦合方法是准确有效的。

3 舵翼结构模型耦合计算结果与分析

3.1 有限元建模

舵翼结构尺寸如图7(a)所示，材料为D6AC高温合金钢，材料参数如表4所示，计算模型流体域、固体域网格分别如图7(b)、图7(c)所示。模型计算域的边界条件如图7(d)所示，壁面边界选择无滑移边界条件，将初始壁温设置为216.65 K，流体域外边界设置为压力远场边界条件，模拟飞行高度为20 km，来流速度5Ma，压强5 529.3 Pa，温度216.65 K，密度0.088 91 kg/m3。固体结构场设置为弹性结构，尾部施加固定约束，对称面添加对称约束，将舵翼结构与外流场交界面设置为流-固耦合面，在界面上进行分区耦合之间的数据传输。根据舵翼结构飞行特点，设置飞行迎角为-10°、-5°、0°、5°、10°进行计算，耦合时间为40 s。

为了分析舵翼结构在近空间飞行状态下的气动力、气动热与结构响应之间的关系，将仿真工况设置为耦合计算(考虑气动力、气动力耦合效应)、仅气动热计算(仅考虑热变形)、仅气动力计算(仅考虑气动力变形)和无结构变形因素计算4种工况。

3.2 结构变形分析

舵翼结构在X、Y、Z 3个方向发生变形，在不同迎角状态下发生拉压、弯曲和热膨胀变形，结构在Z方向只发生热膨胀变形而无拉压和弯曲变形，因此只对舵翼结构X方向(纵向)和Y方向(横向)的变形进行分析。

对舵翼结构施加耦合工况，图8为不同迎角状态下的舵翼结构纵向变形云图，图9为舵翼结构不同位置的纵向变形曲线，其中横坐标为舵翼结构所处位置，X=0处为前缘顶端位置，如图7(a)所示。

舵翼结构在近空间飞行状态下发生气动加热现象，使舵翼结构产生纵向变形。由于舵翼结构尾部采用固定约束，类似于悬臂梁受载方式，当迎角为-10°时，舵翼结构尖角的平分线与来流方向夹角最小，此时舵翼结构弯曲挠度最小，而纵向变形量最大。随着迎角增大挠度随之增加，纵向变形量逐渐减小。

图10为舵翼结构不同位置的横向变形曲线，与纵向变形规律不同，随着迎角增大，横向变形量逐渐增大。舵翼结构尾部采用固定约束，受到横向气动力作用，发生弯曲变形。另一方面，由于气动加热使结构发生热膨胀，迎风侧温度较高，膨胀量大，背风侧温度低，膨胀量小，造成舵翼结构沿背风侧发生弯曲。在气动力和气动热的共同作用下，舵翼结构发生弯曲变形，迎角越大，上述现象越明显，横向变形量也就越大。

对比图9、图10数据可知，当舵翼结构迎角较小时，以纵向变形为主，横向变形较小，随着迎角增大，横向变形明显增大并超过纵向变形。为了进一步分析上述变形规律产生的原因，对气动热、气动力以及耦合作用下的横向、纵向变形进行分析。

舵翼结构纵向变形如图11所示，单纯气动热作用下的变形量略大于耦合作用，二者变化规律一致，相差很小；而气动力对舵翼结构造成压缩变形，变形量随迎角增大而降低，由于变形量非常小，对整体纵向变形的影响可以忽略。因此，舵翼结构纵向变形主要由气动热的作用引起，气动力作用极小，气动热与气动力的耦合作用不明显。

舵翼结构横向变形如图12所示，在气动热因素作用下变形为0.400 42～0.459 53 mm；在气动力作用下横向变形量较纵向变形量明显增大，但仍然很小，仅为0.000 01～0.000 08 mm。当气动热与气动力耦合作用时，变形量明显增大，达到为0.484 30～0.633 66 mm，可见耦合效应在变形量上体现得比较明显。在单纯气动力作用下，舵翼构件并没有发生明显的弯曲变形。在气动热耦合作用下，舵翼结构温度急剧升高，材料物性参数发生明显变化，弹性模量大幅降低，导致发生了更大的变形。弯曲变形改变了结构的气动外形，使温度场进一步升高，一方面加剧了因舵翼结构热膨胀和热弯曲导致的变形，另一方面继续降低了材料的弹性模量，使结构产生了更大的变形。在二者的耦合作用，横向变形量较单独气动热和单独气动力作用下有着明显的升高，迎角越大，气动力与气动热作用的耦合效果越明显。

3.3 温度分析

对舵翼结构施加耦合工况，图13为不同迎角状态下的温度分布云图，舵翼结构前缘处气动热效果最明显，温度最高，沿尾部方向逐渐降低。由于舵翼结构前缘尖角相对气流方向不是对称的，上下面迎风角度不同，导致温度分布也不相同。下表面与来流形成的角度更大，气动热更剧烈，温度更高。随着迎角增大，下表面气动热效果进一步加剧，温度进一步升高。

对舵翼结构温度场受气动热力耦合效应的影响规律进行研究，不同耦合因素下最高温度随迎角变化曲线如图14所示。在3种耦合状态下，温度都随迎角增大而升高。在只考虑气动热因素时，整体温度最低；在无结构变形耦合因素中，只考虑了舵翼结构的热膨胀而不考虑气动力导致的变形，这时温度较只考虑气动热时略高，可见单纯热膨胀变形能够小幅度增加气动热效应；在完全耦合因素中，温度明显大于上述2个工况，且迎角越大差值越高。由于在大迎角姿态飞行时，舵翼结构受气动热和气动力作用发生变形，使迎角进一步增大，导致气动加热作用更加剧烈，进而升高了舵翼结构的温度。可见，迎角越大，气动热与气动力耦合效应对结构温度场的影响就越明显。

4 结论

1) 在近空间飞行状态下，舵翼结构头部尖角处的温度、变形最大，随着飞行迎角增大，温度升高，纵向变形减小，横向变形增大。

2) 在近空间飞行状态下，舵翼结构会发生纵向拉伸和横向弯曲变形，纵向拉伸变形主要由气动加热引起，横向变形主要由气动力作用引起。

3) 在近空间飞行状态下，舵翼结构气动热、气动力的耦合效应使变形与温度增大；飞行迎角越大，耦合效应越明显。因此，当舵翼结构以大迎角状态飞行时，气动热力环境会更加恶劣，在热防护设计和结构设计时应充分考虑耦合效应。

[1] 陈雄昕，刘卫华，罗智胜，等.高超音速飞行器气动热研究进展[J].航空兵器，2014(06)：8-13.

[2] 彭治雨，石义雷，龚红明，等.高超声速气动热预测技术及发展趋势[J].航空学报，2015，36(01)：325-345.

[3] 张灿，胡冬冬，叶蕾，等.2017 年国外高超声速飞行器技术发展综述[J].战术导弹技术，2018(01)：47-50.

[4] 刘敏华，俞启东，陈升泽，等.关于未来导弹战形态及创新设计的研究[J].导弹与航天运载技术，2018(01)：1-6.

[5] 聂雪媛，黄程德，杨国伟.基于 CFD/CSD 耦合的结构几何非线性静气动弹性数值方法研究[J].振动与冲击，2016，35(08)：48-53.

[6] 夏吝时，杨凯威，孔维宣，等.金属尖化前缘模型主动式热疏导性能试验研究[J].装备环境工程，2017，14(10)：82-86.

[7] 吴大方，王岳武，蒲颖，等.高超声速飞行器复合材料翼面结构 1 100 ℃ 高温环境下的热模态试验[J].复合材料学报，2015，32(02)：323-331.

[8] 吴大方，林鹭劲，吴文军，孙陈诚.1 500 ℃极端高温环境下高超声速飞行器轻质隔热材料热/振联合试验研究[J].航空学报，2020，41(07)：223612.