基于TLBM-FVM耦合的飞行器舱内热环境跨尺度预测方法

肖光明，张超，桂业伟，杜雁霞,*，刘磊，魏东

1. 中国空气动力研究与发展中心 空气动力学国家重点实验室，绵阳 621000

2. 西安交通大学 航天航空学院 机械结构强度与振动国家重点实验室，西安 710049

飞行器在大气层中飞行时将面临以空气为流体介质的复杂流动和传热过程，不仅有机体表面高速可压缩流动条件下的强制对流换热现象，同时存在低速不可压缩流动条件下的舱内对流与导热、辐射复合的复杂换热现象。近年来，随着国内外新型飞行器的发展，特别是新一代临近空间飞行器的出现，高空高速、长航时、可重复天地往返等苛刻性能要求不断提出，给飞行器舱内热安全带来了严峻挑战[1]。由于舱内高能密度仪器设备的广泛使用，飞行器内部的热载荷以每5、6年增加一倍的速度增长[2]。不适宜的工作温度会缩短其使用寿命，从而降低系统可靠性。与此同时，鉴于一些飞行器的特殊使命，其机载设备对舱内温度环境也提出了越来越严苛的要求。因此，舱内热环境的准确预测已经成为优化飞行器热控与防热设计、减小系统冗余、提高飞行器有效载荷并保障飞行器热安全的重要基础[3]。

在采用密封式设计的情况下，飞行器舱内流动与传热过程主要表现为复杂空间内多体、分散热源作用下的空气自然对流换热、表面辐射热传递及其与机身结构导热、外部气动加热的多物理机制耦合传热现象[4]。此外，在真实飞行条件下，舱内热环境将受到飞行弹道、设备布局、防热及热控设计等多种因素的综合影响。由于内在传热机制以及外部影响因素等耦合问题的复杂性，飞行器舱内复合传热特性的研究几乎无法做到全区域、全时段的全耦合分析。同时，对舱内热环境的地面试验模拟来说，不仅试验成本高，而且存在热载荷偏差大、重力环境不同、测试设备干扰及测量数据不足等诸多难题。

尽管伴随着计算流体力学(CFD)与计算传热学(NHT)的飞速发展，针对各类复杂流动与传热过程的数值模拟技术均得到了较好的发展与应用[5]，但由于飞行器舱内耦合传热过程反映了不同传热机理及模型的相互作用，单一的数值计算方法在相关问题的研究中存在较大限制，很难同时具备部件级瞬时传热特性与整机级全时温度响应的评估能力。目前，对于飞行器舱内热环境的研究，国内外研究人员通常采用成熟的商业软件(如FLUENT、SINDA/FLUINT、I-DEAS等)进行局部的流动与传热过程模拟[6]，或者基于集总参数化方法的热网络法开展整机热分析[7]，较少关注局部特征与整体性能的相互关联与影响。

针对飞行舱内热环境的精细化预测需求，在国家数值风洞(NNW)工程支持下，本文以FL-CAPTER软件平台[8]为基础，开发了基于热格子Boltzmann方法(TLBM)的舱内对流换热计算模块，建立了可用于导热-对流耦合传热数值仿真的TLBM-有限体积法(FVM)分区求解方法，并以典型飞行器仪器舱为研究对象，开展了同时考虑结构热响应与设备散热的舱内热环境综合分析，为实现飞行器热防护/热管理一体化设计提供一定的技术支撑。

1 基于TLBM的对流换热计算方法

对于低速不可压缩对流换热问题，TLBM可以直接通过流体密度和声速的状态方程进行压力计算，而基于有限体积方法的SIMPLE(Semi-Implicit Method for Pressure Linked Equations)算法则需要额外进行压力修正方程的求解。此外，TLBM对于极端复杂几何形状的边界条件处理更加简洁有效，根据近壁格子点的几何权重前处理结果，采用适当的粒子反射模型可以准确地刻画流体通过壁面的质量、动量及能量流量等边界信息[9-10]。

考虑实际流体中质量、动量和能量传递的不同步性，采用基于多松弛时间模型的双分布热格子Boltzmann方法[11]，密度与温度分布函数分别采用D3Q19和D3Q7离散速度模型[12]，设r为空间位置向量，t为计算时刻，对应(r,t)时空演化方程的离散形式为

fi(r+eiδt,t+δt)-fi(r,t)=Ωi+δtFi

i=0,1,…,18

(1)

gj(r+ejδt,t+δt)-gj(r,t)=Λj

j=0,1,…,6

(2)

式中：ei和ej为离散速度矢量；δt为演化时间步长；Fi为Bolzmann方程的作用力项；Ωi和Λj分别为密度分布函数fi(r,t)和温度分布函数gj(r,t)的碰撞算子，具体表达式为

(3)

(4)

Bolzmann方程的作用力项Fi的具体表达式为

(5)

F=-gβ(T-Tref)

(6)

式中：g为重力加速度矢量；β为流体体积膨胀系数；T和Tref分别为流体实际温度及参考温度。

TLBM在应用于真实飞行器舱内对流换热模拟时还存在复杂流/固界面的高效识别与高保真处理难题，为保证复杂曲面边界的求解精度，进一步基于非平衡态外推数据重构，并结合流/固界面的自适应分辨率识别方法(详见2.2节)发展了无滑移黏性壁面、给定速度、压力及温度等边界类型的通用处理格式。

2 TLBM-FVM耦合方法

2.1 时空耦合策略

针对飞行器结构传热与舱内自然对流换热的同时模拟需求，基于“分区求解、界面耦合”的思想[13]，发展了流体与固体分区求解的TLBM-FVM耦合策略，如图1所示,其中，流体区域采用TLBM方法进行空气对流换热的计算，g(t)为随时间变化的重力加速度，qf(tf)为流体区域计算时刻tf对应的壁面热流；固体区域采用FVM方法进行防热结构与仪器设备内部温度响应的模拟，Ts(ts)为固体区域计算时刻ts对应的壁面温度，关于FVM数值方法的详细介绍可参考文献[14-15]。

图1 TLBM-FVM的分区耦合策略Fig.1 Zonal coupling strategy of TLBM-FVM

当热量在流体与固体之间进行跨介质传递时，由于介质间热响应时间的显著差异，流/固耦合传热过程将呈现时间多尺度效应。特别对于飞行器封闭舱内以重力和温差为主要驱动力的自然对流换热问题，结构及设备表面温度变化将对流场结构产生直接影响，导热与对流两种传热机制之间存在强耦合关联。因此，舱内热环境的预测将无法借鉴可用于飞行器外部气动热与结构传热耦合的准平衡求解策略，即对流换热与结构导热分别采用稳态与非稳态计算模式。

采用TLBM进行自然对流数值模拟时，为保证计算的稳定性，一般要求时间步长Δtf满足

(7)

式中：Mamax为最大模拟马赫数，对于不可压缩流体的TLBM计算，通常要求Mamax≤0.3；Pr和Ra分别为无量纲的普朗特数和瑞利数；Lref和ΔlTLBM分别为计算参考长度与单个LBM格子对应的真实物理长度；v为流体黏性系数。

在采用Runge-Kutta方法进行固体导热的非稳态求解时，时间步长Δts需要满足von Neumann稳定条件：

(8)

式中：ΔlFVM为有限体积网格大小；αs为固体的热扩散系数。

根据计算时间步长Δtf和Δts，同样可以初步判断热扰动在不同区域传播的特征时间大小。为兼顾耦合计算的精度及效率，进一步考虑流/固耦合边界热扰动特征时间Δtc(如热载荷及温度变化率等)，发展了异步迭代、协同推进的TLBM-FVM时间耦合方法，如图2所示,图中Step 1)′表示下一轮迭代的Step 1)。具体步骤如下：

图2 TLBM-FVM的时间耦合策略Fig.2 Time coupling strategy of TLBM-FVM

1) 根据当前时刻及前一时刻的TLBM计算结果给定FVM求解所需的边界条件，流体域的界面信息向固体域传递。

2) FVM时间推进，推进步数M的具体计算表达式为

M=max(1, Δtc/Δts, Δtf/Δts)

(9)

3) 根据FVM计算结果，给定TLBM求解所需的边界条件，固体域的界面信息向流体域传递。

4) TLBM时间演化，演化步数N为

N=max(1, Δtc/Δtf, Δts/Δtf)

(10)

5) 重复步骤1)～4)，直到计算结束时刻tend。

2.2 界面处理方法

在建立了TLBM-FVM时空耦合策略的基础上，为提高TLBM对于复杂曲面边界的适应性，进一步发展了自适应分辨率的TLBM-FVM耦合界面识别算法，并基于FL-CAPTER软件具有的高鲁棒性数据插值方法实现了不同求解域之间的高效高精度D-N(Dirichlet-Neumann)信息交互。

图3给出了热格子Boltzmann计算网格与有限体积计算网格示意图，BEl,n为第n块FVM网格对应的第l个边界单元，NDijk,m为第m块TLBM网格对应的(i,j,k)格子点(实心代表流体域，空心代表固体域)，ISijkm,ln为边界单元BEl,n与邻近格子点NDijk,m方向射线的交点。可以看出，界面识别作为TLBM-FVM耦合计算前处理的关键问题，其核心在于如何快速有效地求解FVM网格边界单元与TLBM网格的交点ISijkm,ln，判断格子点NDijk,m的流/固区域属性，建立边界单元BEl,n与对应格子点NDijk,m的相互映射关系，从而为非平衡态外推数据重构提供几何权重，并指定界面数据插值的交互网格点信息。

图3 TLBM-FVM的耦合界面示意图Fig.3 Schematic diagram of TLBM-FVM coupling interface

图4为TLBM-FVM耦合界面前处理的具体操作流程。首先，利用传统的结构网格或非结构网格生成技术完成固体区域的FVM网格生成，并提取对应的边界单元信息。其次，根据边界单元所属FVM网格的分块信息及几何大小自适应设置TLBM网格的区域划分及空间分辨率，完成格子点生成。在计算格子点不同方向射线与边界单元交点之前，先进行边界单元相对格子空间的区域自适应，通过格子射线与边界单元交点范围的合理判断，可有效提升求解效率。然后，在获得交点信息的基础上，完成流体域、固体域及边界等不同类型格子点的标识，找到边界格子点至边界单元的单向映射关系。最后，判断边界格子点至边界单元的单向映射是否为满射，即判断是否所有的边界单元都具有对应的边界格子点，当边界单元较小时，有可能与格子各方向射线均不存在交点，此时需要进一步搜索邻近的流体格子点与之对应，进而完成TLBM格子点与FVM网格边界单元之间双向映射关系的构建。

图4 TLBM-FVM耦合界面前处理流程图Fig.4 Pre-processing flow chart of TLBM-FVM coupling interface

3 耦合算例验证

为验证TLBM-FVM耦合方法的有效性，对经典的有限厚度平板表面对流换热与结构传热耦合问题进行了求解[16-17]，如图5所示，图中δT(x)和δ(x)分别表示温度边界层厚度及速度边界层厚度，平板厚度δs=10 mm，参考长度L=200 mm，底部恒温边界设置Tb=600 K，平板热导率λs=0.288 4 W/(m·K)。入口不可压缩空气的来流速度U∞=12 m/s，温度T∞=1 000 K，出口压力P∞=1.03×105Pa。分别选取空气密度ρ=0.352 5 kg/m3、普朗特数Pr=0.662 9、动力黏性系数μ=3.95×10-5kg/(m·s)，则平板表面流动的参考雷诺数ReL为

(11)

图5 有限厚度平板的耦合传热Fig.5 Conjugate heat transfer of a flat plate with finite thickness

图6给出了基于TLBM-FVM耦合方法计算得到的平板表面对流换热系数沿水平方向的分布曲线及中间位置(x=100 mm)沿厚度方向温度分布曲线。可以看出，计算结果与Luikov[16]基于微分传热方程(DHT)及边界层方程(BL)得到的理论值均吻合较好。

图6 TLBM-FVM耦合验证算例计算结果Fig.6 Calculation results of TLBM-FVM coupled verification example

4 舱内热环境耦合分析算例

4.1 几何模型

如图7所示，针对典型飞行器仪器舱开展舱内热环境耦合分析建模研究，其几何模型主要包括防热结构及各类仪器设备。防热结构由CMCs(Ceramic Matrix Composites)一体化尖锐前缘、高温合金防热层及隔热内衬3部分组成，舱内则装填有10个仪器设备，其中设备4和设备8为供电装置，其他大部分为测控设备。模拟舱段的总长为1.5 m，总宽为0.5 m，前缘半径为3 mm，楔角为10°，舱段防/隔热层总厚度为15 mm。

图7 典型飞行器仪器舱结构模型、设备布局及基本尺寸Fig.7 Structural model, equipment layout and basic size of a typical aircraft instrument cabin

重点研究结构/设备热传导与舱内空气对流换热的耦合传热特性，暂时不考虑不同设备之间的直接接触，各设备主要通过空气对流进行间接换热，同时忽略设备表面的辐射换热影响。

4.2 计算参数

飞行器仪器舱的外部热载荷数据主要通过FL-CAPTER软件平台的数值气动热计算模块分析获得。当飞行器以马赫数10、攻角0°在高度30 km处进行巡航飞行时，300 K冷壁面条件下的表面热流Qc及绝热壁面条件的表面温度Tinf计算结果如图8所示，图中黑色线为飞行器气动外形的轮廓线，包括迎风面(Windward)和背风面(Leeward)。

图8 典型飞行器仪器舱气动热载荷条件Fig.8 Aerodynamic heat load conditions of a typical aircraft instrument cabin

在不进行气动热耦合求解的前提下，为进一步考虑结构表面温度升高后壁温对气动加热热流及表面辐射散热的影响，防热结构外部净热载荷Qin的计算条件采用热壁修正形式：

(12)

式中：Tw和Tamb分别为结构壁温及环境温度；ε为表面辐射散热系数；σ为Stephan-Boltzmann常数。对于非气动加热面，如飞行器仪器舱侧面及尾部端面，均采用绝热边界进行处理。

表1给出了各舱内仪器设备的发热功率，考虑设备4和设备8发热量较大，温升较快，通常需要对其进行额外的温度控制，模拟时可以考虑直接设置为恒温边界。研究中设备4和设备8的表面温度均取为320 K，同时整个计算过程中设备都处于开启状态。

表1 舱内设备发热功率Table 1 Heating power of cabin equipment

除了内/外热载荷外，压力及重力条件也会对飞行器舱内热环境分析计算产生重要影响。分别选取舱内压力为1个大气压，重力加速度大小为9.81 m/s2。经初步估算，当最大温差为50 K时，以飞行器结构舱段总长为参考长度，表征自然对流强弱的瑞利数将高达109量级，向基于TLBM的对流换热模拟能力提出了较大挑战。

4.3 计算结果与分析

针对4.1节中典型飞行器舱内热环境耦合分析模型进一步开展了TLBM-FVM耦合数值模拟能力测试，FVM的计算网格规模为15.7万，耦合边界单元数约2万。为保证高Ra数下TLBM的计算精度及稳定性，其模拟的格子总数达10.6亿。

经测试，采用本文发展的自适应分辨率界面识别方法，完成FVM的耦合边界单元与TLBM网格的交点计算、区域判断及双向映射建立的总时长仅为100 s左右，计算机的CPU型号为Intel©Xeon©Gold 6244，大幅提高了TLBM-FVM方法对复杂工程问题的模拟效率。

同时，选取耦合特征时间Δtc=1 s，计算结束时刻tend=70 s。TLBM求解器采用8张32 G显存的Tesla V100 GPU卡进行加速运算，完成TLBM-FVM耦合计算的CPU总时长约为480 h。

图9和图10展示了基于TLBM-FVM分区耦合方法的飞行器结构传热响应与舱内自然对流温度场及流场的计算结果，同时给出了耦合界面的热流及温度交互情况。由图9可以看出，受气动热载荷的直接影响，防热结构的高温区主要集中在尖前缘与第2级压缩面上，且70 s后压缩面上的热量已渗透至仪器舱内部，对应内表面的温度明显升高。同时，设备4与设备8作为仪器舱内的主要高温面，在浮升力驱动下将产生上升热羽流，从而使设备上部的结构温度升高。由图10可以看出，尽管在内/外载荷的共同作用下，舱内流场结构比较复杂，但由于内部空气温度整体较低，自然对流换热较为强烈的区域仍然集中在高温设备附近。

图9 FVM计算结果及耦合界面传递数据(t=70 s)Fig.9 Calculation results and coupling interface transfer data of FVM (t=70 s)

图10 TLBM计算结果及耦合界面传递数据(t=70 s)Fig.10 Calculation results and coupling interface transfer data of TLBM (t=70 s)

研究结果表明，采用TLBM-FVM耦合方法进行飞行器舱内热环境数值模拟，不仅能同时获得舱内热量整体传输性能及局部分布特征，而且可以综合考虑防热结构设计与仪器设备布局对舱内热环境的影响，TLBM-FVM耦合方法可为实现飞行器舱内/外一体化热管理设计提供分析工具。

5 结 论

本文发展了基于多松弛时间模型的热格子Boltzmann方法，提出了分区求解、协同推进的时空耦合策略及自适应分辨率的流/固界面识别算法，建立了基于TLBM-FVM耦合的飞行器舱内热环境跨尺度预测方法，并以典型飞行器仪器舱为对象开展了耦合分析模型研究及数值模拟，验证了耦合计算方法的有效性。得到的主要结论如下：

1) 采用自适应分辨率的流/固界面识别算法可以实现十亿量级以上大规模TLBM计算网格的高效预处理，包括交点求解、区域判断及边界映射等，为进一步发展TLBM-FVM耦合界面的数据重构及交互技术奠定基础。

2) 在内/外热载荷的共同作用下，飞行器舱内热环境的形成及演化机制通常较为复杂，为深入探索舱内热量及耗散机理，实现给定约束条件下舱内流动及传热过程的优化设计，有待于借助TLBM-FVM耦合分析手段系统开展防热结构、设备布局、环境参数等不同因素对舱内热环境的影响规律研究。

致 谢

感谢西安交通大学王娴教授研究团队及李明佳教授研究团队对热格子Boltzmann数值模拟技术发展提供的建议与支持。