基于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數值模擬技術發展提供的建議與支持。