熱防護系統(tǒng)分區(qū)協(xié)調(diào)耦合推進方法

黃 杰，姚衛(wèi)星，陳 炎,3，孔 斌,3

(1. 南京航空航天大學飛行器先進設計技術(shù)國防重點學科實驗室，南京 210016；2. 南京航空航天大學機械結(jié)構(gòu)力學及控制國家重點實驗室，南京 210016；3. 中國航空工業(yè)集團公司成都飛機設計研究所，成都 610091)

0 引 言

空天飛行器再入階段受到巨大的氣動加熱作用，為保證內(nèi)部機體結(jié)構(gòu)在可承受的溫度范圍內(nèi)，需要在機體表面附加熱防護系統(tǒng)(Thermal protection system,TPS)[1-3]，熱防護系統(tǒng)的設計為高超聲速飛行器關(guān)鍵技術(shù)之一，而陶瓷防熱瓦是應用最廣泛的隔熱結(jié)構(gòu)之一，它通過應變隔離墊粘接在機體表面，其熱控性能直接影響到飛行器的安全，必須進行熱防護系統(tǒng)的熱分析。

傳統(tǒng)的熱防護系統(tǒng)熱控分析方法將氣動熱與結(jié)構(gòu)傳熱人為地分割開來。文獻[4-7]在已知熱流密度的情況下進行了飛行器頭部防熱瓦、蓋板式熱防護系統(tǒng)以及多層隔熱結(jié)構(gòu)的熱分析，獲得了熱防護系統(tǒng)和機體結(jié)構(gòu)的溫度-時間歷程。其求解過程為先進行氣動加熱的分析，再將熱流密度作為邊界條件直接分析傳熱方程求解熱防護系統(tǒng)的溫度場，而忽略氣動熱與結(jié)構(gòu)傳熱之間的耦合效應，這必然會影響熱防護系統(tǒng)熱控性能的分析精度。

由于空天飛行器再入過程中產(chǎn)生巨大的氣動加熱會造成結(jié)構(gòu)溫度急劇升高，而結(jié)構(gòu)溫度升高后，激波層和邊界層內(nèi)氣體與熱防護系統(tǒng)外表面之間的溫度梯度將會減小，造成壁面熱流密度的降低，即氣動加熱與結(jié)構(gòu)溫度場存在強烈的耦合效應，應給予足夠的重視。

耦合分析的關(guān)鍵在于氣動熱的準確計算、耦合迭代的策略以及耦合面之間的數(shù)據(jù)插值。以往的研究常常采用Eckert[8]提出的參考焓法等工程算法進行熱流密度的分析，其分析精度有限，現(xiàn)代氣動熱分析采用計算流體力學(Computational fluid dynamics,CFD)方法求解Navier-Stokes方程[9-10]，可借助CFD++或FLUENT等軟件或自編程序準確分析駐點及其他區(qū)域的熱流密度分布。而耦合分析常采用迭代求解法，其核心技術(shù)是要保證瞬態(tài)耦合的時間精度和迭代的協(xié)調(diào)性。此外由于結(jié)構(gòu)傳熱網(wǎng)格單元的尺度遠大于流體單元，耦合面節(jié)點非一一對應，故耦合面之間的數(shù)據(jù)交換需要采用插值法完成，其關(guān)鍵技術(shù)是要保證插值精度和耦合面熱流通量的守恒性。

本文在充分考慮氣動熱與結(jié)構(gòu)傳熱耦合的基礎(chǔ)上，提出了熱防護系統(tǒng)熱控性能分析的分區(qū)協(xié)調(diào)耦合推進方法，對比了圓管算例的計算結(jié)果與試驗結(jié)果，最后進行了防熱瓦與主動冷卻系統(tǒng)同時存在時的熱防護系統(tǒng)熱控性能及其影響因素分析。

1 流體和結(jié)構(gòu)傳熱的數(shù)值算法

1.1 流體分析數(shù)值算法

不考慮體積力和內(nèi)熱源情況下，直角坐標系下的流體動力學N-S方程的積分形式為：

(1)

式中：W為守恒向量，F(xiàn)c為對流通量，F(xiàn)v為黏性通量，dS為控制體V的邊界面，n為邊界面dS的外法線單位向量。將式(1)按有限體積法(Finite volume method,FVM)進行空間離散可得：

(2)

式中：Wi和Vi分別為控制體i的守恒向量和體積，NF為控制體邊界面的數(shù)目，ΔSN為第N個邊界面的面積。由于對流通量Fc具有高度非線性特點，并且集中體現(xiàn)了流場的對流特征，采用具有TVD(Total Variation Diminishing)性質(zhì)的NND格式[11]對其進行空間離散，其中半離散化的上風型NND格式為：

(3)

(4)

(5)

(6)

式中：sgn為符號函數(shù)，min表示取較小值。

控制方程中的黏性項采用中心差分格式離散，湍流模型采用Menter’s SST兩方程模型[12]，而針對非定常問題的時間離散，在n+1時刻采用時間二階精度的隱式三點后差離散可得到二階精度的離散方程：

(7)

(8)

式中：Δτ和Δt分別為虛擬時間步長和物理時間步長，稱為雙時間步長法[13]；p和n分別為虛擬時間迭代步和物理時間迭代步。虛擬時間步上的內(nèi)迭代可采用LU-SGS格式[14]求解，當p→∞時虛擬時間項趨近于零，式(8)的定常解即為二階精度的非定常解。

1.2 結(jié)構(gòu)傳熱數(shù)值算法

在無體積熱源的假設下結(jié)構(gòu)瞬態(tài)熱傳導的控制方程為：

(9)

式中：ρ0為結(jié)構(gòu)材料密度，c為材料比熱容，kx，ky和kz分別為材料三個方向的導熱系數(shù)。其中比熱容和導熱系數(shù)一般為溫度的函數(shù)。

針對本文熱防護系統(tǒng)的熱分析問題，其外表面邊界條件為壁面熱流密度Qaero和壁面熱輻射量Qrad，其表達式分別為：

(10)

(11)

式中：?T/?n為壁面法向溫度梯度，ε為壁面熱輻射率，玻爾茲曼常數(shù)σ=5.67×10-8W/(m2·K4)，Twall為壁面溫度，Tat為大氣環(huán)境溫度。

對式(9)進行有限元離散可得到總體合成矩陣求解方程：

(12)

(13)

求解式(13)即可得到結(jié)構(gòu)各個時刻節(jié)點的瞬時溫度。

2 耦合分析模型和策略

空天飛行器再入過程中熱防護系統(tǒng)外壁面巨大的氣動加熱效應會造成結(jié)構(gòu)溫度急劇升高，而結(jié)構(gòu)溫度升高后，激波層和邊界層內(nèi)氣體與壁面之間的溫度梯度將會減小，造成壁面熱流密度的降低，即氣動加熱與結(jié)構(gòu)溫度場存在強烈的相互耦合效應，如圖1所示。

本文采用分區(qū)協(xié)調(diào)耦合推進方法進行氣動熱與結(jié)構(gòu)傳熱的分析，其耦合推進策略如圖2所示，即將流場與結(jié)構(gòu)分為兩個不同的區(qū)域，分開建模與求解，其中流場采用FVM求解，而結(jié)構(gòu)熱傳導采用有限元法(Finite element method,FEM)求解，其基本假設和特點為：

1) 在任意n到n+1時間步內(nèi)FVM(或FEM)求解過程中壁面溫度(或壁面熱流密度)不變，即各子學科求解時邊界條件凍結(jié)。

2) 各子學科在相同的時間節(jié)點進行數(shù)據(jù)交換，以保證耦合分析的協(xié)調(diào)性以及耦合時間精度。

根據(jù)圖2所示的耦合推進策略，本文提出了適用于瞬態(tài)和穩(wěn)態(tài)的熱防護系統(tǒng)熱分析流程(見圖3)，主要步驟為：

1) 首先建立相應的FVM和FEM數(shù)值分析模型，定義來流條件和初始溫度T0，進行定常氣動熱分析，結(jié)果作為耦合分析的熱流密度初始條件。

2) 進行第n到n+1時間步的流場求解，將計算獲得的熱流密度Qn+1通過插值方法傳遞給結(jié)構(gòu)FEM模型外表面。

3) 進行第n到n+1時間步的結(jié)構(gòu)FEM傳熱分析，輸出分析獲得的結(jié)構(gòu)溫度場Tn+1。

4) 判斷時間是否達到分析的總時間ttotal(針對瞬態(tài)分析)，或壁面熱流密度和結(jié)構(gòu)溫度場是否收斂(針對穩(wěn)態(tài)分析)，若達到ttotal或收斂，結(jié)束耦合分析，否則進行熱流密度和壁面溫度的數(shù)據(jù)傳遞。

5) 返回2)，進行下一個迭代步的求解，直到分析時間達到ttotal或壁面熱流密度和結(jié)構(gòu)溫度場收斂，結(jié)束分析。

3 耦合面數(shù)據(jù)傳遞

以上氣動熱和結(jié)構(gòu)傳熱的耦合分析中涉及到熱流密度和壁面溫度的數(shù)據(jù)傳遞，若這兩個學科之間的數(shù)值模型網(wǎng)格節(jié)點一一對應，那么通過節(jié)點之間的對應關(guān)系就可以方便地進行熱流密度和壁面溫度的數(shù)據(jù)傳遞。但實際分析中結(jié)構(gòu)網(wǎng)格尺度通常遠大于流體網(wǎng)格，即FVM和FEM耦合面節(jié)點不一致，如圖4所示，故需要在耦合面進行插值以實現(xiàn)數(shù)據(jù)傳遞。

本文的分析涉及到熱流密度和結(jié)構(gòu)壁面溫度的數(shù)據(jù)插值，其中的關(guān)鍵技術(shù)是要保證耦合面上熱流密度通量的守恒性。

(14)

式中：Q為FVM網(wǎng)格面的熱流密度，q為與FVM網(wǎng)格節(jié)點對應的結(jié)構(gòu)傳熱FEM網(wǎng)格面的熱流密度。當Scoup為耦合面時，式(14)體現(xiàn)了熱流密度通量的總體守恒性。本文采用如下基于控制面的雙向映射插值方法進行熱流密度的數(shù)據(jù)傳遞，算法主要步驟為：

1) 將FVM和FEM耦合面上節(jié)點從物理空間(x,y,z)通過坐標轉(zhuǎn)換映射到參數(shù)空間(u,v)中，即u=u(x,y,z)和v=v(x,y,z)，它將三維曲面上的空間節(jié)點轉(zhuǎn)換到二維的控制面上。

2) 在物理空間中采用Kdtree算法搜索任意FEM網(wǎng)格節(jié)點ζi(x,y,z)附近的FVM網(wǎng)格節(jié)點ηi(x,y,z)，并將其轉(zhuǎn)換到二維的控制面得到ζi(u,v)和ηi(u,v)。

3) 將FVM網(wǎng)格節(jié)點坐標ηi(u,v)和相應的熱流Qi(u,v)代入到如下三次插值函數(shù)Q(u,v)，采用最小二乘法求解插值函數(shù)的系數(shù)bj(j=1,2,…,10)。

Q(u,v)=b1u3+b2v3+b3u2v+b4v2u+

b5u2+b6v2+b7uv+b8u+b9v+b10

(15)

4) 將所有FEM網(wǎng)格節(jié)點ζi(u,v)代入到已知系數(shù)bj(j=1,2,…,10)的插值函數(shù)Q(u,v)中，即可求得FEM網(wǎng)格節(jié)點插值熱流密度qi(u,v)，并通過式(14)進行熱流密度通量的守恒性檢驗。

4 算例1：耦合方法的驗證

本文選取NASA高超聲速圓管風洞試驗模型[15]進行算例分析。圓管內(nèi)徑R1=25.4 mm，外徑R2=38.1 mm，圓管材料為不銹鋼，其密度ρ=8030 Kg/m3，導熱系數(shù)k=16.72 W/(m·K)，比熱容c=502.48 J/(kg·K)，此外圓管初始壁面溫度T0=294.4 K。高超聲速來流馬赫數(shù)Ma=6.47，來流溫度T=241.5 K，來流壓強P=648.1 Pa，攻角α=0°。建立了二維分析模型，如圖5所示，其中流場網(wǎng)格量為60×100，壁面第一層網(wǎng)格高度Δh=1×10-5m，并且在激波處進行了局部加密。圓管結(jié)構(gòu)有限元網(wǎng)格量為40×25，采用四節(jié)點平面單元模擬。此外分析類型為瞬態(tài)分析，時間步長Δt=1×10-4s，總時間ttotal=2 s，并且將定常流場作為耦合分析的初始條件。

分析獲得了圓管駐點溫度隨時間的變化曲線，如圖6所示，從圖中可觀察到曲線斜率在初始階段很大，隨時間的推移逐漸減小，說明結(jié)構(gòu)開始階段升溫很快，隨后溫度上升速度逐漸減緩，最后趨于一個穩(wěn)定的溫度值(穩(wěn)態(tài)解)。

圖7為2 s時刻圓管外壁面的相對溫度分析情況,從圖中可知壁面耦合分析結(jié)果分布曲線與試驗結(jié)果分布曲線吻合良好，其中駐點處溫度耦合分析結(jié)果Tstag為442 K，而試驗結(jié)果為465 K，相對誤差為4.95%。從而驗證了分區(qū)協(xié)調(diào)耦合推進方法的正確性與分析精度，其可用于分析熱防護系統(tǒng)氣動熱與結(jié)構(gòu)傳熱的耦合問題。

5 算例2：空天飛行器頭錐熱防護系統(tǒng)耦合分析

現(xiàn)代空天飛行器同時采用防熱瓦及主動冷卻系統(tǒng)以保證蒙皮內(nèi)表面溫度穩(wěn)定在設定值，如圖8所示。本文進行了空天飛行器頭錐熱防護系統(tǒng)的穩(wěn)態(tài)耦合分析。其中防熱瓦外表面涂層熱輻射率ε=0.85，LI-900防熱瓦厚度h1=60 mm，應變隔離墊(SIP)厚度h2=4.5 mm，鋁合金蒙皮厚度h3=1.5 mm，蒙皮導熱系數(shù)k=180 W/(m·K)。防熱瓦和SIP的導熱系數(shù)均隨溫度呈非線性變化趨勢[16]。飛行馬赫數(shù)Ma=7，高度H=60 km，攻角α=25°，初始結(jié)構(gòu)溫度T0=300 K，大氣環(huán)境溫度Tat=247.1 K。建立了三維分析模型，流場和結(jié)構(gòu)網(wǎng)格如圖9所示，結(jié)構(gòu)熱分析采用八節(jié)點六面體單元模擬，蒙皮內(nèi)表面溫度約束在325 K，為此主動冷卻系統(tǒng)需要持續(xù)的提供功率以保證蒙皮內(nèi)表面維持在此溫度。

利用本文的耦合方法進行了空天飛行器頭錐熱防護系統(tǒng)的分析。圖10為非耦合和耦合方法分析獲得的對稱壁面處的熱流密度分布曲線，可知非耦合分析方法獲得的駐點熱流密度比耦合方法高了53.8%。這是由于非耦合方法未考慮壁面溫度升高造成激波層和邊界層與壁面溫度梯度減小。

圖11和和12分別為非耦合和耦合方法分析得到的防熱瓦及SIP外表面溫度云圖。其中非耦合和耦合方法獲得的防熱瓦最高溫度分別為1108.5 K和994.1 K，相差了114.4 K；而非耦合和耦合方法獲得的SIP最高溫度分別為472.9 K和440.3 K，相差了32.6 K。故非耦合方法獲得的防熱瓦和SIP溫度偏高，這是由于非耦合方法考慮壁面溫度升高對氣動熱的反饋作用造成的。

研究了飛行馬赫數(shù)Ma對防熱瓦和SIP最高溫度、主動冷卻系統(tǒng)最大功率的影響，如圖13所示，從圖中可觀察到馬赫數(shù)從5增加到9，防熱瓦最高溫度上升了69.2%，SIP最高溫度上升了36.1%，而最大冷卻功率上升了266.1%。圖14為以上響應隨熱防護系統(tǒng)外表面涂層熱輻射率ε的變化情況，熱輻射率從0.25增加到1，防熱瓦最高溫度降低了24.5%，SIP最高溫度降低了18.2%，而最大冷卻功率降低了50.1%，故應盡量采用高輻射率的涂層以降低熱防護系統(tǒng)的厚度和主動冷卻系統(tǒng)的功率，達到減重的目的。

在進行熱防護系統(tǒng)的設計時，防熱瓦的導熱系數(shù)和厚度的選取直接影響到熱控性能，是熱防護系統(tǒng)最重要的設計參數(shù)。圖15和圖16分別為防熱瓦和SIP最高溫度、主動冷卻系統(tǒng)最大功率隨防熱瓦導熱系數(shù)比k/k0以及厚度h1的變化曲線，k0代表防熱瓦原始導熱系數(shù)。從圖中可觀察到防熱瓦導熱系數(shù)比從0.128增加到8時，防熱瓦最高溫度降低了3.2%，而防熱瓦厚度從20 mm增加到80 mm，防熱瓦最高溫度上升了1.5%，即防熱瓦導熱系數(shù)和厚度對其最高溫度影響很小。隨著防熱瓦導熱系數(shù)的增加和厚度的降低，外部氣動熱能量能夠更容易地傳遞到SIP及蒙皮，造成SIP溫度升高以及主動冷卻系統(tǒng)功率增加，并且影響顯著。故在防熱瓦重量約束條件下應盡量采取較厚和較低導熱系數(shù)的防熱瓦，以提高熱防護系統(tǒng)的隔熱性能和降低主動冷卻系統(tǒng)的重量。

6 結(jié) 論

1) 本文針對空天飛行器熱防護系統(tǒng)熱控性能研究，提出了分區(qū)協(xié)調(diào)耦合推進分析方法，其中氣動熱和結(jié)構(gòu)傳熱分別采用FVM和FEM求解，且在耦合面進行數(shù)據(jù)插值以實現(xiàn)數(shù)據(jù)傳遞。

2) 進行了圓管算例分析，圓管壁面熱流密度和溫度分布與試驗結(jié)果吻合良好，從而驗證了本文的耦合推進方法的正確性與分析精度。

3) 研究了空天飛行器頭錐熱防護系統(tǒng)的熱控性能，非耦合方法獲得的防熱瓦和SIP溫度相對耦合方法偏高，這是由于非耦合方法未考慮壁面溫度升高對氣動熱的反饋作用，而耦合方法充分考慮了此影響。

[1] Olynick D. Trajectory-based thermal protection system sizing for an X-33 winged vehicle concept [J]. Journal of Spacecraft and Rockets, 1998, 35(3): 249-257.

[2] Shideler J L, Webb G L, Pittman C M. Verification tests of durable thermal protection system concepts [J]. Journal of Spacecraft and Rockets, 1985, 22(6): 598-604.

[3] Oscar A M, Anurag S, Bhavani V S, et al. Thermal force and moment determination of an integrated thermal protection system [J]. AIAA Journal, 2010, 48(1): 119-128.

[4] 劉雙, 張博明. 發(fā)汗式主動冷卻金屬熱防護系統(tǒng)主動冷卻效率研究[J]. 宇航學報, 2011, 32(2):433-438. [Liu Shuang, Zhang Bo-ming. Investigation on transpiration active cooling metallic thermal protection systems [J]. Journal of Astronautics, 2011, 32(2): 433-438.]

[5] 孟松鶴, 楊強, 霍施宇, 等. 一體化熱防護技術(shù)現(xiàn)狀和發(fā)展趨勢[J]. 宇航學報, 2013, 34(10):1295-1302. [Meng Song-he, Yang Qiang, Huo Shi-yu, et al. State-of-arts and trend of integrated thermal protection system [J]. Journal of Astronautics, 2013, 34(10): 1295-1302.]

[6] 魏東, 杜雁霞, 石友安, 等. 非均勻氣動加熱下隔熱層結(jié)構(gòu)的優(yōu)化設計[J]. 宇航學報, 2015, 36(10):1108-1113. [Wei Dong, Du Yan-xia, Shi You-an, et al. Optimization design of heat insulation layer with non-uniform heat flow loads [J]. Journal of Astronautics, 2015, 36(10): 1108-1113.]

[7] 朱言旦, 劉偉, 曾磊, 等. 飛行器結(jié)構(gòu)部件導熱/輻射耦合傳熱特性預測方法[J]. 宇航學報, 2016, 37(11):1371-1377. [Zhu Yan-dan, Liu Wei, Zeng Lei, et al. Method of predicting conduction-radiation coupled heat transfer characteri-stics for vehicle structural component [J]. Journal of Astron-autics, 2016, 37(11): 1371-1377.]

[8] Eckert E R G. Engineering relations for friction and heat transfer to surfaces in high velocity flow [J]. Journal of Aerospace Sciences, 1955, 22(8): 585-587.

[9] Ju Y. Lower-upper scheme for chemically reacting flow with finite rate chemistry [J]. AIAA Journal, 1985, 33(8): 1418-1425.

[10] 楊肖峰, 唐偉, 桂業(yè)偉, 等. 火星環(huán)境高超聲速催化加熱特性[J]. 宇航學報, 2017, 38(2):205-211. [Yang Xiao-feng, Tang Wei, Gui Ye-wei, et al. Hypersonic catalytic aeroheating characteristics for Mars entry process [J]. Journal of Astronautics, 2017, 38(2): 205-211.]

[11] 沈清, 顧鋼民, 高樹椿, 等. NND格式在航天飛機頭部段N-S方程求解中的應用[J]. 空氣動力學學報, 1989, 7(2):145-155. [Shen Qing, Gu Gang-min, Gao Shu-chun, et al. Applications of the NND-scheme to the solving of Navier-Stokes equations on the fore-head of space-shuttle-like body [J]. Atca Aerodynamica Sinica, 1989, 7(2): 145-155.]

[12] Menter F R. Two-equation eddy-viscosity turbulence models for engineering applications [J]. AIAA Journal, 1994, 32(8): 1598-1605.

[13] 王保國, 劉淑艷, 張雅, 等. 雙時間步長加權(quán)ENO-強緊致高分辨率格式及在葉輪機械非定常流動中的應用[J]. 航空動力學報, 2005, 20(4):534-539. [Wang Bao-guo, Liu Shu-yan, Zhang Ya, et al. Pseudo-time method of weighted ENO-strongly compact high resolution schemes and its application to turbomachinery unsteady flow [J]. Journal of Aerospace Power, 2005, 20(4): 534-539.]

[14] Yoon S, Jameson A. Low-upper Gauss-Sediel method for the Euler and Navier-Stokes equations [J]. AIAA Journal, 1988, 26(9): 1025-1026.

[15] Wieting A R. Experimental study of shock wave interference heating on a cylindrical leading edge [R]. NASA-TM-100484, 1987.

[16] Ng W H, Friedmann P P, Waas A. Thermomechanical analysis of a thermal protection system with defects and heat shorts [R]. AIAA-2006-2212, 2006.