考虑力/热/结构多场耦合效应的飞行弹道预测

代光月，曾磊，刘深深，冯毅，唐伟，桂业伟

中国空气动力研究与发展中心 计算空气动力研究所，绵阳 621000

新一代临近空间高超声速飞行器被期望能够执行快速全球到达任务，成为近年来兴起的研究热点。这类飞行器需在高马赫数下长时间巡航，在设计时除需要在布局方面尽可能提高升阻比外，还要求结构能够尽可能的做到轻质，所以越来越多的轻质材料或大型薄壁设计结构被应用于高超声速巡航飞行器。这种轻质的大型薄壁结构在高马赫数长时间飞行条件下，飞行环境、结构温度、结构应力/应变、飞行弹道之间的多物理场耦合严重，强烈的气动加热会造成结构温度场改变，进而引起材料属性、几何变形、结构应力、模态和刚度变化，而防热结构性能和外形的改变又会对气动力、热环境带来较大的影响。对于无控飞行器来说，气动力的改变会直接影响飞行轨道，导致飞行与预测轨道发生偏离；对于有控飞行器来说，气动力变化会对飞行器的配平控制舵偏产生影响，而高超声速飞行器在舵面控制能力方面一般都比较有限，如多场耦合带来的气动力变化超出了舵面的控制能力，会影响飞行安全，所以在飞行器设计阶段，气动力/热/结构多物理场耦合带来的影响必须给予考虑。

在传统的飞行器研制模式中，气动力/热/结构多场耦合分析主要用于飞行器热防护系统设计，而飞行弹道和控制系统的设计与气动布局和防热结构多场耦合分析是解耦进行的，多场耦合研究往往以提供热走廊的形式作为弹道设计边界。随着临近空间高超声速飞行器对预测精度要求的不断提高，传统的“气动力—轨道参数—热环境—热响应”单向顺序解耦分析方法不能完全满足工程实际需要，需要发展考虑气动力/热/结构/弹道多物理场耦合的一体化分析方法。

在气动力/热/结构多物理场耦合研究方面，自20世纪80年代以来，国内外学者开展了大量研究工作，发展了多种多场耦合求解方法。早在1958年，Roger[1]就对典型气动力/热/结构多场耦合问题涉及到的各物理因素间的相互耦合关系进行了总结，分析了不同物理场间的强弱耦合关系；Dechaumphai等[2-4]最早开展了相关的数值方法和试验验证研究。Michopoulos和Farhat[5]于2005年提出了多场、多域和多尺度的概念，并建立了四场(气动环境场、结构温度场、结构应力场和结构应变场)两域(流体域和固体域)的耦合数学模型。Culler和McNamara[6]对高超声速气动力/热/结构耦合问题中的耦合关系及单/双向耦合方法进行了较为细致的研究，相关结果显示，随着飞行时间增长，双向耦合的影响会逐渐增大，通过采用准静态的耦合求解方式，可大幅降低计算量。在采用全数值方法开展多场耦合研究方面，Lohner等[7]通过对耦合过程涉及的物理场控制方程及其求解方式进行详细分析，提出了紧耦合和松耦合的概念，并基于其定义的松耦合算法提出了一套适用于气动力/热/结构多场耦合的松耦合计算策略，开发了相应的数值模拟软件。

国内相关工作开展相对较晚，但发展迅速。中国空气动力研究与发展中心桂业伟等针对气动力/热/结构多场耦合问题，提出了新的多场耦合关系和特有的时空耦合概念[8]，在空间上，将多场耦合分为面耦合、体耦合和面体耦合；在时间上，提出了时间全耦合、时间松耦合和时间修正耦合的不同耦合层次，并针对临近空间高超声速飞行器工程应用背景，开发了FL-CAPTER热环境/热响应耦合计算平台[9]。西北工业大学张伟伟等[10]采用松耦合的方法建立了热气动弹性仿真模型，将气动力/热/结构多场耦合问题分解成了气动热计算、结构温度场计算、热结构计算和气动弹性计算，在时间域内实现了高超声速热气动弹性的仿真。南京航空航天大学张兵和韩景龙[11]采用共享内存技术的方式开发了适用于通用CFD软件和有限元软件的多场耦合计算平台，并针对典型高超声速机翼进行了气弹分析。哈尔滨工业大学周印佳等[12]采用分区求解方法，针对超高温陶瓷材料的耦合传热问题进行了数值研究。耿湘人[13]、季卫栋[14]等在气动加热/结构传热分析一体化数值计算方法方面开展了大量工作。

总体来说，国内外针对气动力/热/结构多场耦合问题，提出了多种研究方法，分析的思路也相对比较清晰，难点主要是工程应用时如何兼顾耦合求解的精度与效率问题。现有气动力/热/结构多场耦合问题的研究主要集中在两个方面，一是从耦合模型方面如何提升多场耦合求解分析的精度和效率，二是相关高超声速飞行器热气动弹性问题的研究，包括静热气动弹性问题和动热气动弹性问题。目前的耦合模型忽略了结构变形和振动对热环境的影响及结构热生成与弹性变形之间的热力学耦合[15]，考虑气动力/热/结构多场耦合效应带来的飞行弹道变化则相对较少。

因此，本文针对临近空间高超声速飞行器气动力/热/结构多场耦合问题，基于FL-CAPTER软件平台现有多场耦合分析能力，将其拓展至飞行力学领域，构建了一套考虑气动力/热/结构多场耦合效应的弹道模拟新方法，并针对给定舵偏角下自主配平控制的助推-压缩楔组合体外形，开展了不同耦合时间尺度下的飞行弹道特性分析，初步探讨研究多场耦合效应对飞行弹道的影响。

1 多场耦合效应下的弹道预测方法

1.1 FL-CAPTER软件多场耦合分析方法

FL-CAPTER软件是中国空气动力研究与发展中心自主研发的热环境/热响应耦合计算分析平台[9]。该软件采用基于优选锚点方案的逐段热环境加沿弹道结构温度场和典型时刻应力应变场的计算模式，支持试验状态或给定弹道条件下的气动力/热/结构多场耦合计算；软件经过了大量标模算例与型号外形验证，具备较好的计算精度，关于该软件的详细介绍可参见文献[16-20]。本文新构建的考虑多场耦合效应的弹道仿真方法基于该软件内设的多场耦合分析方法发展而来，此处对耦合方法进行简单说明。

图1展示了FL-CAPTER软件采用的多场耦合分析方法。计算时，流场通过准定常求解三维Navier-Stokes方程获得，温度场通过数值求解三维热传导方程获得，应力应变场通过数值求解弹性力学微分方程组获得。其中，ti是第i阶段的起始时刻，ti+1是第i阶段的结束时刻；ΔtG=ti+1-ti为第i阶段考虑变形的间隔时间，可随i值变化而变化；ΔtT为温度场计算的时间步长；qw(t)为t时刻的表面热流；Tw,i为第i阶段起始时刻的壁面温度；ωi为第i阶段起始时刻的结构变形量。在每一阶段，均同时计算ti和ti+1时刻的流场；当温度场以ΔtT的时间步长向前推进时，对应时间曲线上的环境边界通过两端流场解插值获取。结构应力应变分析视为准稳态，针对每一阶段的结束时刻ti+1开展应力应变场求解并据此更新计算外形。同时，计算过程中充分考虑交界面上的数据传递，流场计算时均假定存在ti时刻的壁温分布，并引入对应温度条件下的热壁修正公式[21]以加快收敛速度：

(1)

式中：Ti为ti时刻的温度分布；Ti+N·Δi为从ti时刻推进N个时间迭代步时的结构温度分布；q为热流；h为焓值；hre为恢复焓。

图1 FL-CAPTER软件采用的多场耦合分析方法Fig.1 Coupling analysis method used in FL-CAPTER software

1.2 弹道仿真方法及验证

高超声速飞行器飞行的高度和速度都比较高且航程较长，在轨道仿真过程中不能按照不动的圆球形大地的情况进行分析，必须要考虑地球扁率及地球自转的影响，地球扁率影响使得重力加速度的方向并没有始终指向地心，分析飞行器运动变得更加复杂。而地球自转影响则给飞行器带来了附加的牵连加速度和科氏加速度，进而影响飞行器的飞行轨迹。

考虑地球的自转及地球扁率时，动力学方程必须补充惯性牵连力和哥式惯性力，此时质心运动的动力学方程为[22]

(2)

式中：F、G及P分别为作用在飞行器质心上的气动力、重力及发动机推力；m为飞行器的质量；V为速度；t为时间；-mωe为惯性牵连力；-mωk为科氏惯性力。

对于惯性牵连力，有

(3)

式中：

(4)

其中：Ωe为地球旋转角速度；φe为飞行器对应的地球纬度；μe为发射方位角；Re为飞行器对应点下的地球半径。

对于科氏惯性力，有

(5)

式中：

(6)

考虑到地球的扁率及飞行器飞行高度变化，飞行器的重量在不同地理位置和高度也各不相同，即

(7)

式中：r为飞行器的地心向径；R0x、R0y、R0z为坐标原点地心向径的分量。

描述飞行器运动的微分方程组为一阶常微分方程组，一般很难得到解析解，常用的求解方法是在给定初始条件的情况下，通过数值的方法进行求解，本文选用常用的四阶Runge-Kutta方法。

1.3 考虑多场耦合影响的弹道预测

本节在FL-CAPTER软件平台现有已知弹道多场耦合分析能力的基础上，将其拓展至飞行力学领域，建立一套能够考虑气动力/热/结构多场耦合效应影响的弹道预测新方法。由于问题的复杂性，研究中暂不考虑舵面控制的影响，即所有弹道均假设为固定舵偏条件下的自主配平飞行。

基于FL-CAPTER软件多场耦合平台，具体通过如下步骤实现：

步骤1依据飞行器初始外形在初始时刻t0所在位置对应的来流条件计算初始时刻的气动环境。

步骤2假设每隔dt时间考虑一次多场耦合效应，近似认为在dt时间间隔内外形不发生变化，基于飞行力学的飞行运动方程，计算出从初始时刻t0～t1=t0+dt时刻的飞行弹道。

步骤3基于算得的t1时刻的位置及来流条件，计算t1时刻的气动环境。

步骤4基于FL-CAPTER软件，开展基于t0～t1时刻弹道条件下的多场耦合计算，更新飞行器外形。

步骤5计算新外形在t1时刻的气动环境。

步骤6重复步骤2～步骤5直至完成整个飞行过程。

需注意的是，dt为考虑多场耦合变形影响的计算时间步长，即：气动/弹道耦合时间步长，并不是弹道计算时间步长。显然，当dt非常小时，即与弹道计算时间步长量级相当时，可近似认为该飞行弹道实时考虑了多场耦合效应带来的变形影响，计算结果将非常接近物理事实，更为真实可靠。然而，弹道计算时间步长通常为10-3s量级，甚至更小，在当前计算条件下，这样的计算量是难以接受的。随着计算机技术的不断发展、飞行器气动力/热/结构多场耦合计算方法的不断进步，未来dt的取值足够小是完全可能的。对比分析考虑多场耦合变形前后的飞行弹道，可获得多场耦合效应对飞行弹道的影响。

图2 考虑力/热/结构多场耦合影响的弹道仿真计算流程Fig.2 Process of trajectory simulation considering fluid-thermal-structural coupling effects

2 计算外形与初始弹道

本文针对助推-压缩楔组合体外形开展气动力/热/结构多场耦合效应对飞行弹道影响的研究，计算外形如图3所示。其中，前体部分为一个两级压缩楔，总长约530 mm；在此基础上，为满足气动稳定性要求，增加了一级助推器，助推器与压缩楔间由过渡段连接。助推级直径约200 mm，长度约1 480 mm，连接段长度为200 mm。助推器一级安装有4片平板舵，用于稳定控制，平板舵呈十字布局，为满足配平需要，初始状态时一级助推上的气动舵面预置了8°舵偏。需要说明的是，此外形主要用于基本方法和影响规律的研究分析，与真实飞行器有一定区别。

在弹道仿真分析时，假设弹道起始高度为30 km， 模型重量为300 kg，俯仰方向转动惯量为1 400 kg·m2，以马赫数Ma=3.0作为弹道终点参数。考虑助推结束后的无动力滑翔下降过程，组合体外形在无动力、固定舵偏角条件下飞行，且只考虑俯仰方向运动。仿真得到的初始弹道曲线如图4所示。可以看到，由于气动阻力作用，飞行马赫数逐渐降低，从马赫数6.0～马赫数3.0的过程中，共飞行了约200 s，飞行高度H降低了约9 km，飞行距离L约为250 km。对于此构型，随着马赫数的降低，配平迎角α逐渐减小，起始状态马赫数为6.0时，配平迎角α约为17.5°，马赫数为3.0时，配平迎角约为17.2°，迎角随时间在平衡迎角附近振荡收敛。

图3 助推-压缩楔组合体外形Fig.3 Configuration of boost-compression-wedge model

图4 分析模型的初始弹道Fig.4 Initial trajectory of analyzed model

在多场耦合效应分析时，考虑到外形数值计算周期较长且变形主要发生在前体，仅对前体压缩楔外形开展气动力/热/结构多场耦合计算与分析，一级助推器主要起稳定配平的作用，假定在多场耦合计算时助推级外形不发生变化。

3 计算结果与分析

综合计算时间和效率等因素，本文共考虑3种气动力/热/结构多场耦合效应与飞行弹道的耦合方案，即分别取dt=5，20，40 s考虑一次多场耦合效应影响，在所取的时间间隔内，假定结构外形保持不变。

3.1 前体压缩楔多场耦合效应

首先，针对两级压缩楔外表面不同位置选取了温度关注点，分析了关注点处结构温度随时间的变化历程。关注点位置如图5所示，图中模型为半模外形，其中，点A位于模型侧面与前缘交界位置，点B位于前缘驻点中心线位置，点C、D、E、F位于压缩拐角区域，点G、H位于背风面区域。图6给出了不同耦合方案下关注点的温升历程。可以看出，前缘温度在短时间内迅速上升至一较高水平，之后则随着飞行马赫数和飞行高度的不断降低，温度逐渐回落，下降十分明显；压缩拐角区域也存在温度先逐渐上升再逐渐回落的趋势，但上升速度较前缘位置更慢，下降程度也更低；背风区随着时间的增长，温度逐渐升高，但始终保持在较低的温度水平。对不同耦合方案，随着耦合时间步长的增大，计算的温度值逐渐偏高。

图5 关注点位置Fig.5 Sketch map of concerns

图7给出了不同耦合方案下结构最大位移随时间的变化曲线，最大位移位置位于前缘处。可以看出，压缩拐角区域与背风区巨大的温差诱导拐角区域出现大的结构热应力应变，最终引起压缩楔外形尤其是前缘位置出现较大位移，温差越大位移量越大。对当前飞行条件，结构最大位移的变化趋势与迎风面温度变化趋势保持一致，均为先快速增大后缓慢减小直至结束。位移量的最大值出现在约50 s时刻，此时前缘上翘约18.5 mm。对比不同时间步长的耦合方案，可以看出，结构位移量对耦合时间步长并不敏感，3种耦合方案下最大位移量量值差异不超过1 mm。

图6 不同耦合方案下的关注点温升历程Fig.6 Temperature histories of concerns under different coupling schemes

图7 不同耦合方案下结构最大位移量时间历程Fig.7 Time history of maximum structural displacement under different coupling schemes

3.2 不同耦合方案对飞行弹道的影响

在获得结构位移量的基础上，开展了多场耦合效应对飞行弹道的影响研究。首先基于位移量，计算得到了变形后新外形的静态气动力数据。图8给出了不同耦合方案下对应时间的配平迎角。可以看出，初始配平迎角约为17.5°，受多场耦合效应影响，压缩楔前缘不断上翘，前体升力增大，在给定舵偏角条件下，抬头力矩作用使得配平迎角进一步增大，特别是由于飞行器前缘位置力臂较长，变形带来的力矩增量相对较大，导致配平迎角随着变形量的增大急剧增大，在50 s附近，变形后外形对应的配平迎角约为34°。之后，随着变形量的逐渐回落，配平迎角呈逐渐减小趋势。

图8 不同耦合方案配平迎角变化Fig.8 Variation of trimming angels of attack under different coupling schemes

图9给出了不同耦合方案下气动力特性随时间的变化曲线。可以看出，初始弹道下飞行器的升阻力变化较小，当考虑多场耦合效应带来的结构变形影响之后，前缘上翘导致升力、阻力同时增大，但升阻比降低，且随着气动/弹道耦合时间间隔的增大，气动特性振荡加剧。

图9 不同耦合方案下气动力特性变化Fig.9 Variation of aerodynamic characteristics under different coupling schemes

图10给出了考虑变形前后固定舵偏角时飞行弹道的变化比较。可以看出，考虑多场耦合效应后，前体上翘导致配平迎角增大，并带来了一定的升力增量，使得飞行高度增加；同时，飞行过程中的阻力也随之增加，对于无动力情形，飞行中消耗的能量增加，飞行器减速效应更明显；由于变形后升阻比有所降低，航程变短。此外，从飞行迎角随时间的变化曲线可以看出，即使本文采用固定舵偏的自主配平飞行策略，变形后的飞行器会自动寻找新的配平位置，并在其附近呈振荡收敛趋势，表明考虑热变形后的外形仍然是稳定的，如果外形的静稳定裕度不够或者考虑热变形后外形处于中立稳定或者静不稳定状态，飞行器将可能出现迎角发散等破坏性结果。

对比不同气动/弹道耦合时间步长下的气动力特性变化曲线，可以看出，当耦合时间步长过大时，突然施加的结构变形量较大，气动力的阶跃类似于冲击载荷，会对弹道带来大的影响，迎角振荡的幅度很大，收敛到新的配平迎角的时间增长，在还没有达到新的迎角收敛状态时开始下一轮的耦合计算，会导致振荡的叠加。实际上，飞行器真实飞行过程，多场耦合效应影响带来的结构变形应是一个渐变过程，不会是相隔较长时间的一个剧烈冲击，这也就是说，当选取的时间步长过大时，会带来非物理的振荡，导致计算结果失真。由此可见，在气动力/热/结构多场耦合效应与飞行弹道的耦合过程中，耦合时间步长的选取非常关键，须兼顾计算效率与精度。

图10 考虑变形前后弹道变化Fig.10 Variation of trajectory before and after deformation

3.3 修正耦合方案下的结果对比

通过3.2节的分析可以看出，在弹道仿真时过大的气动/弹道耦合时间步长会带来非物理的振荡，导致计算结果失真；但气动/弹道耦合时间步长的减小会使得计算量显著增大，计算周期明显延长。为解决这一问题，在研究过程中，注意到不同时间步长耦合方案下的结构变形量差异较小，对本文计算，3种耦合方案下最大位移量量值差异不超过1 mm。为此，针对气动力/热/结构多场耦合效应和飞行弹道耦合时时间步长的选取问题进一步展开研究，提出了基于变形量回溯插值技术的双时间步修正耦合方法，具体实施步骤如下:

步骤1兼顾计算效率和精度，获得较大时间步长下的结构变形量。

步骤2针对时间步长过大、外形突变带来的非物理振荡现象，假定结构变形是按照某种规律渐变(如：线性变化)，以曲线拟合的方式(如线性插值)将变形量插值到较小的时间步长内。

步骤3基于插值后的变形量获取较小时间步长条件下的新外形，计算变形后的气动力数据，更新弹道数据。

为检验上述方法的可行性，假定结构变形量随时间呈线性变化，利用线性插值的方法，将耦合时间步长20 s方案得到的变形量和耦合时间步长40 s方案得到的变形量进行插值处理，得到修正近似耦合时间步长5 s方案下的变形量，并与真实计算的5 s时间步长方案下的气动特性和弹道特性进行了对比分析。

图11给出了真实计算的5 s时间步长耦合方案和修正近似耦合方案下的气动力特性对比曲线。可以看出，采用本文修正耦合方法计算得到的升阻力特性与5 s步长耦合方案计算结果规律基本保持一致，由于间隔步长过大引起的非物理振荡不再明显。图12给出真实计算的5 s时间步长耦合方案和修正近似耦合方案下的弹道特性对比曲线。可以看出，基于变形量插值的修正耦合方法由于变形带来的非物理冲击振荡明显减弱，整体弹道与5 s间隔耦合计算方案基本一致，验证了本文提出的基于变形量回溯插值的修正耦合方法的可行性，表明该方法在精度允许范围内，可大幅提高计算效率。

图11 5 s时间步长耦合方案和修正近似耦合 方案下气动力特性比较Fig.11 Comparison of aerodynamic characteristics under 5 s time step coupling scheme and modified approximation scheme

图12 5 s时间步长耦合方案和修正近似耦合方案下弹道特性对比Fig.12 Comparison of flight trajectory characteristics under 5 s time step coupling scheme and modified approximation scheme

4 结 论

本文初步探讨与研究了多场耦合效应对飞行弹道的影响。将通常的气动力/热、结构传热、结构应力/应变的耦合进一步拓展至与飞行弹道的耦合，构建了一套考虑气动力/热/结构多场耦合效应的弹道预测新方法，实现了考虑气动力/热/结构多场耦合效应的飞行弹道预测功能。

1) 建立了考虑多场耦合效应的飞行弹道预测方法。通过对助推-压缩楔外形开展多场耦合效应对飞行剖面的影响研究，表明本文建立的耦合求解方法可用于飞行器弹道仿真。

2) 对于助推-压缩楔组合体外形，考虑多场耦合效应后，变形将带来配平迎角增大，飞行器升力、阻力同时增大，升阻比降低，弹道飞行高度增加，飞行马赫数降低，航程变短等一系列影响，这也验证了更精准的弹道预测需综合考虑多场耦合效应带来的影响。

3) 气动/弹道耦合计算时间步长选取对弹道仿真结果存在较大影响，当步长选取过大时，可能会带来非物理振荡，导致计算结果失真。

4) 兼顾计算效率与精度，提出了基于变形量回溯插值技术的修正耦合方法。仿真结果表明，修正后的弹道仿真精度得到明显提升，削弱了由于时间步长选取过大造成的非物理振荡问题。

无论是旨在飞行器防热结构设计的气动力/热/结构多场耦合研究还是飞行器的飞行弹道设计，两者都是难度很深、工作量很大的复杂工作，本文仅仅是在二者强耦合需求背景下考虑二者耦合的一些初步尝试。需注意的是，本文研究的飞行器虽能够反映高超声速流动特征，但与真实飞行器在尺寸、模型质量和惯量等参数方面都存在一定区别，结果表明多场耦合效应对飞行器配平与弹道存在一定影响，但定量的分析需基于本文方法针对具体的真实飞行器才具有实际意义。同时，本文目前所开展的研究中飞行器假定为给定舵偏角下的自主配平控制，没有考虑控制系统带来的影响。下一步，不仅要对所发现的气动/弹道耦合计算时间步长展开进一步的影响研究，还将继续深化飞行弹道的耦合研究，分析不同弹道参数与气动力/热/结构多场耦合特性的影响规律，逐步拓展考虑控制影响的耦合规律。