李宗陽，李 煜，程廣益，竇怡彬，張曉宏

(上海機電工程研究所，上海 201109)

0 引 言

臨近空間飛行器關(guān)鍵技術(shù)包括飛行器總體設(shè)計技術(shù)、多學(xué)科優(yōu)化設(shè)計技術(shù)、超燃沖壓發(fā)動機技術(shù)以及長時間高溫?zé)岱雷o技術(shù)[1-2]。其中，長時間高溫?zé)岱雷o技術(shù)涉及到熱防護材料的研發(fā)、熱防護結(jié)構(gòu)的設(shè)計以及熱防護結(jié)構(gòu)的試驗驗證等方面。臨近空間飛行器熱環(huán)境的確定是進行熱防護設(shè)計的前提和依據(jù)，數(shù)值仿真和地面試驗是獲取臨近空間飛行器熱環(huán)境以及進行熱防護設(shè)計和驗證的主要手段。由于多場耦合效應(yīng)的多場耦合分析技術(shù)可以更為準(zhǔn)確地分析飛行器的熱環(huán)境，因此，國內(nèi)外學(xué)者在多場耦合領(lǐng)域做了大量的工作。1987年，美國學(xué)者ALLAN等針對飛行器前緣氣動熱環(huán)境以及激波干擾分布情況等問題，開展了圓柱氣動加熱試驗，詳細(xì)研究激波相互干擾的流場分布以及結(jié)構(gòu)傳熱情況[3-4]。2010年，趙曉利等考慮流動傳熱和結(jié)構(gòu)傳熱之間的相互影響，對圓管氣動加熱試驗進行了三維數(shù)值模擬，得到的結(jié)果與試驗值符合，實現(xiàn)了流動傳熱耦合的準(zhǔn)確計算[5]。2011年，CULLER等根據(jù)多場耦合分析中各物理場的特征時間，采用準(zhǔn)靜態(tài)多場耦合分析方法對C/C復(fù)合材料進行多場耦合分析，并比較了單/雙向耦合分析、時間步長以及迭代步數(shù)對結(jié)果的影響[6]。2014年，德國學(xué)者FRAUHOLZ等采用基于有限元求解器的多場耦合分析方法，模擬了二維進氣道流動傳熱耦合問題，并針對進氣道壁面為不同材料以及不同內(nèi)壁溫度的情況，展開對進氣道性能影響的分析[7]。

飛行器頭罩是整個飛行器氣動熱最為嚴(yán)酷的部位之一，因此頭罩的熱防護結(jié)構(gòu)設(shè)計也成為飛行器熱防護工作研究的重點[8]。本文在充分考慮各物理場之間的強弱耦合關(guān)系的基礎(chǔ)上，應(yīng)用將考慮流動傳熱全耦合的有限體積法和熱結(jié)構(gòu)弱耦合的有限元法相結(jié)合的多物理場耦合分析方法，對不同形式的飛行器頭罩熱防護結(jié)構(gòu)進行數(shù)值分析，得到了飛行器頭罩的熱環(huán)境以及頭罩熱防護結(jié)構(gòu)內(nèi)部的溫度和熱應(yīng)力分布，可為飛行器頭罩的熱防護設(shè)計提供參考。

1 數(shù)值計算方法

基于有限體積法的流動傳熱全耦合計算方法是將流體區(qū)域流動傳熱控制方程和固體區(qū)域傳熱控制方程寫成統(tǒng)一的格式，采用統(tǒng)一的離散格式進行離散，并在同一求解器中進行求解，從而實現(xiàn)流體區(qū)域流動傳熱和固體區(qū)域的傳熱全耦合計算。在進行流體區(qū)域計算時，連續(xù)性方程采用可壓縮流動方程，動量方程采用雷諾平均的Navier-Stokes方程，不考慮內(nèi)部熱源、質(zhì)量力以及熱輻射，流體區(qū)域流動傳熱控制方程的微分形式可以統(tǒng)一表述為

(1)

式中：W是守恒變量；Fi是無粘通量；Gi是粘性通量；Xi是坐標(biāo)分量，表達(dá)式分別為

(2)

式中：ρ表示流體密度；u、v、w分別表示空間坐標(biāo)3個方向上的流動速度；i、j、k分別表示空間坐標(biāo)3個方向的單位矢量；E表示流體控制體的總能量；p表示流體壓力；τ表示流體粘性應(yīng)力張量；q表示熱流密度。同樣，在不考慮內(nèi)部熱源以及熱輻射的前提下，固體區(qū)域傳熱控制方程的微分形式可以描述為

(3)

式中：ρs表示固體材料的密度；k是材料的熱傳導(dǎo)系數(shù)；表示拉普拉斯算子；T代表溫度；h代表焓，本文計算所采用的氣體模型是熱完全氣體模型，其焓值是隨溫度變化的函數(shù)，可表示為

(4)

式中：Cp為定壓比熱容；Tref表示參考溫度。

結(jié)構(gòu)溫度變化會引起結(jié)構(gòu)體積發(fā)生變化，同時，內(nèi)部溫度場分布不均勻會導(dǎo)致結(jié)構(gòu)產(chǎn)生的熱變形不協(xié)調(diào)，因此會產(chǎn)生熱應(yīng)力。結(jié)構(gòu)的總應(yīng)變由彈性應(yīng)變和熱應(yīng)變組成，總應(yīng)變張量可表示為

(5)

胡克定律闡述了物體應(yīng)力和應(yīng)變之間的關(guān)系，即彈性體的本構(gòu)關(guān)系。假設(shè)固體材料是各向同性材料，彈性應(yīng)變張量可以表示為

(6)

式中：ET為材料的彈性模量，其值是隨溫度變化的；σij、σkk為應(yīng)力張量；υ為材料的泊松比；δij表示克羅內(nèi)克符號。

假設(shè)材料熱膨脹是線性的，其熱應(yīng)變張量為

(7)

式中，αT為材料的熱膨脹系數(shù)，其值是隨溫度變化的。

2 多物理場耦合策略

本文所采用的多物理場耦合計算方法是涉及流體和結(jié)構(gòu)的熱流固耦合計算，其中流體與固體區(qū)域的流動傳熱耦合計算是基于有限體積法并采用統(tǒng)一形式的方程進行求解的，同時需要滿足流體和固體交界面上的熱通量和溫度的守恒，可以表示為

Tfluid=Tsolid

(8)

(9)

結(jié)構(gòu)熱變形的計算采用單向耦合的計算方法，即將流動傳熱耦合計算出的結(jié)構(gòu)表面壓力以及結(jié)構(gòu)內(nèi)部溫度作為初始變量，展開結(jié)構(gòu)熱力耦合分析。本文假設(shè)結(jié)構(gòu)發(fā)生的變形是小量，忽略結(jié)構(gòu)變形對流場的影響。

3 多場耦合方法驗證算例

3.1 風(fēng)洞試驗概述

圓管前緣氣動加熱試驗是多場耦合領(lǐng)域的一個非常經(jīng)典的試驗。該試驗是1987年ALLAN等在NASA蘭利研究中心的高焓風(fēng)洞中完成的，目的是研究二元進氣道前緣的結(jié)構(gòu)傳熱特性。試驗給出了圓管前緣壁面的壓力和熱通量分布數(shù)據(jù)，常作為驗證多場耦合計算方法準(zhǔn)確性的算例。試驗中不銹鋼圓管的外徑是38.1 mm，內(nèi)徑是25.4 mm，文獻[3]中詳細(xì)列出了不銹鋼圓管壁面的壓力分布和冷壁熱流分布。

3.2 計算模型及網(wǎng)格

數(shù)值模擬的來流條件在表1中給出，其中，p∞、T∞分別指自由來流的壓力和溫度。不銹鋼的材料參數(shù)在表2中給出，其中，λ是材料的熱導(dǎo)率。湍流模型對壁面冷壁熱流的計算精確度有很大的影響。本文采用MENTER提出并改進過的剪切應(yīng)力傳輸(shear stress transfer, SST)k-ω湍流模型[9]。該湍流模型綜合了k-ε湍流模型和k-ω湍流模型的優(yōu)點，可以使SSTk-ω湍流模型同時適用于分離流和附著流，并且具有較高的精確度。近壁面網(wǎng)格的分布對近壁面溫度以及壁面冷壁熱流的影響比較大，文獻[10]指出可以采用當(dāng)?shù)鼐W(wǎng)格雷諾數(shù)來確定近壁面第一層網(wǎng)格高度，其計算公式為

Relocal=ρUΔn0/μ

(10)

式中：ρ為流體的密度；U為流體相對于壁面的切向速度；μ為動力粘性系數(shù)；Δn0為第一層網(wǎng)格高度。為了保證計算結(jié)果的精確度，當(dāng)?shù)鼐W(wǎng)格雷諾數(shù)應(yīng)不大于3，因此近壁面第一層網(wǎng)格的高度為1×10-6m。采用ICEM軟件劃分網(wǎng)格，并用自適應(yīng)網(wǎng)格技術(shù)對激波發(fā)生區(qū)域的網(wǎng)格做加密處理，網(wǎng)格模型如圖1所示。

表1 自由來流條件Tab.1 Free stream condition

表2 材料物理屬性Tab.2 Physical properties of material

圖1 流場及結(jié)構(gòu)域網(wǎng)格Fig.1 Grids of flow field and structural domain

3.3 數(shù)值模擬結(jié)果分析

在流動傳熱耦合計算時，設(shè)定固體區(qū)域初始溫度為294.4 K，時間步長為1×10-4s。圖2是采用本文的計算方法得到的流場密度云圖與試驗結(jié)果紋影圖的對比。圖2中，流體密度突變的位置就是圓管前緣脫體激波發(fā)生的位置，試驗紋影中圓管前端黑色線條即脫體激波。通過對比可以發(fā)現(xiàn)，數(shù)值模擬的結(jié)果中脫體激波的位置以及形狀都與試驗結(jié)果吻合。

圖2 計算得到的流場密度與試驗結(jié)果對比Fig.2 Comparison of calculated flow field density and test result

圓管繞流試驗的數(shù)據(jù)處理方法是采用駐點的壓力p0和熱流密度值q0對圓管壁面的數(shù)據(jù)進行歸一化處理。駐點為坐標(biāo)原點，橫坐標(biāo)為圓管壁面各點和圓心連線與駐點和圓心連線之間的夾角θ。圖3是圓管壁面壓力p分布曲線，可以看出本算例數(shù)值計算的結(jié)果與試驗值是吻合的，圓管前緣駐點處壓力值是35 636.6 Pa，與試驗結(jié)果的誤差在6%以內(nèi)，這說明采用本文多場耦合計算方法計算出的氣動壓力數(shù)據(jù)是可信的。

圖4(a)是圓管壁面冷壁熱流計算結(jié)果與試驗結(jié)果的對比曲線圖，圖4(b)是圓管壁面熱壁熱流計算結(jié)果與試驗結(jié)果的對比曲線圖,q為壁面熱流值。從圓管壁面冷壁熱流與試驗結(jié)果的對比曲線中可以看出，二者吻合得很好，其中數(shù)值計算的駐點熱流值是653 320.9 W/m2，與試驗結(jié)果誤差在2.55%以內(nèi)。從不同時刻圓管壁面熱壁熱流與試驗結(jié)果的對比曲線中可以發(fā)現(xiàn)，熱壁熱流變化明顯的位置在駐點附近，這是因為駐點附近結(jié)構(gòu)表面的溫度變化較大。

圖3 圓管壁面壓力分布曲線Fig.3 Curve of pressure distribution on the cylinder wall

圖5表示的是不同時刻壁面溫度分布曲線，其中實點代表的是5.0 s時刻試驗得到的壁面溫度，這與數(shù)值計算出的壁面溫度非常符合，整體誤差在2.35%以內(nèi)。

(a)冷壁熱流分布 (b)熱壁熱流分布 圖4 壁面熱流分布曲線Fig.4 Distribution curve of the wall heat flux

圖5 不同時刻壁面溫度分布曲線Fig.5 Temperature distribution curve at different time

基于von Mises應(yīng)力準(zhǔn)則，采用有限元法進行結(jié)構(gòu)熱力響應(yīng)分析，分別給出2 s時刻和5 s時刻的圓管等效應(yīng)力云圖，如圖6所示。圖6(a)為2 s時刻的等效應(yīng)力云圖，圖6(b)為5 s時刻的等效應(yīng)力云圖?？梢钥闯觯畲蟮刃?yīng)力出現(xiàn)在圓管的前緣部分，并且最大等效應(yīng)力隨著時間推移逐漸增大，在2 s時刻最大等效應(yīng)力是275.4 MPa，在5 s時刻最大等效應(yīng)力是330.78 MPa，這是因為圓管結(jié)構(gòu)的溫度逐漸升高。

(a) 2 s時刻 (b) 5 s時刻圖6 結(jié)構(gòu)等效應(yīng)力分布云圖Fig. 6 Equivalent stress distribution of structure

4 飛行器頭罩?jǐn)?shù)值分析

針對飛行器頭罩展開多場耦合計算分析，按金屬結(jié)構(gòu)和多層熱防護結(jié)構(gòu)兩種結(jié)構(gòu)類型進行分析，得到頭罩結(jié)構(gòu)溫度場以及應(yīng)力應(yīng)變場，為飛行器頭罩結(jié)構(gòu)熱防護設(shè)計提供參考。

4.1 幾何模型及網(wǎng)格

飛行器頭罩的幾何外形和網(wǎng)格模型如圖7所示，其中多層熱防護結(jié)構(gòu)頭罩外形尺寸與金屬結(jié)構(gòu)一致，由外向內(nèi)第一層是耐燒蝕層，第二層是隔熱層，第三層是金屬層。采用ICEM軟件劃分流體域和結(jié)構(gòu)域一體化的結(jié)構(gòu)網(wǎng)格，并在激波發(fā)生區(qū)域進行網(wǎng)格加密處理，第一層網(wǎng)格高度是2×10-6m，網(wǎng)格增長率是1.15，網(wǎng)格總數(shù)是19萬，其中流體區(qū)域15萬，固體區(qū)域4萬。

4.2 數(shù)值計算與分析

本文對飛行器頭罩進行多場耦合分析，分別對金屬頭罩和防熱頭罩進行數(shù)值分析，得到飛行器頭罩的溫度場、應(yīng)力場以及熱變形情況，并對兩種結(jié)構(gòu)形式頭罩的數(shù)據(jù)進行對比分析。

圖7 飛行器頭罩網(wǎng)格模型Fig.7 Grid model of the aircraft hood

4.2.1 流動傳熱耦合分析

圖8是流場特性云圖，分別是流場的速度、壓力和溫度分布云圖。在飛行器頭罩前部形成脫體激波，受到脫體激波的作用，駐點前部區(qū)域氣流的速度降低，壓力和溫度升高，駐點處的壓力為4.557 MPa，駐點處的溫度是1 610 K，并在此處形成一處低速高溫高壓區(qū)域，會產(chǎn)生嚴(yán)酷的氣動加熱效應(yīng)。

(a) 速度 (b) 壓力 (c) 溫度圖8 流場特性云圖Fig.8 Contours of flow field characteristics

圖9是300 s時刻兩種飛行器頭罩結(jié)構(gòu)溫度分布云圖，從圖中可以看出金屬結(jié)構(gòu)的飛行器頭罩在300 s時刻內(nèi)壁駐點溫度達(dá)到1 505 K，內(nèi)壁面平均溫度達(dá)到1 390 K，結(jié)構(gòu)整體內(nèi)外溫度梯度較小，整體溫度較高。防熱結(jié)構(gòu)的飛行器頭罩在300 s時刻外壁面駐點溫度是1 602 K,內(nèi)壁面駐點溫度是576 K，內(nèi)壁面平均溫度達(dá)到516 K，結(jié)構(gòu)整體內(nèi)外溫度梯度較大，整體溫度較低，熱量主要集中在耐燒蝕層和隔熱層。

(a) 金屬頭罩 (b) 防熱頭罩圖9 300 s時刻飛行器頭罩溫度分布Fig.9 Temperature distribution of the aircraft hood at 300 s

為了研究不同結(jié)構(gòu)形式的飛行器頭罩結(jié)構(gòu)溫度隨時間的變化情況，在圖10中給出了兩種結(jié)構(gòu)形式的頭罩外壁面駐點溫度隨時間變化曲線，可以看出防熱頭罩的外壁面溫升比金屬頭罩外壁溫升更快，且在300 s時刻溫度基本一致。圖11中給出了300 s時刻不同位置處的內(nèi)壁面溫度變化曲線，駐點處的X坐標(biāo)是0 m，飛行器頭罩根部X坐標(biāo)是0.20 m，整體上看金屬頭罩內(nèi)壁面溫度比防熱頭罩內(nèi)壁面溫度高約900 K，可見防熱頭罩的熱防護效果明顯。

4.2.2 結(jié)構(gòu)熱流耦合分析

將流動傳熱耦合計算得到的氣動壓力載荷和溫度載荷作為初始條件，采用有限元方法研究飛行器頭罩結(jié)構(gòu)熱力響應(yīng)情況。圖12是兩種結(jié)構(gòu)形式的飛行器頭罩在300 s時刻的等效應(yīng)力分布云圖。由圖12可見，兩種結(jié)構(gòu)形式的應(yīng)力場分布規(guī)律基本一致，金屬頭罩前緣內(nèi)部位置最大應(yīng)力為870 MPa，防熱頭罩前緣內(nèi)部位最大應(yīng)力為342 MPa，可見防熱結(jié)構(gòu)可以降低結(jié)構(gòu)的應(yīng)力。

圖10 飛行器頭罩外壁面駐點溫度曲線Fig.10 Temperature distribution curve of the stagnation point at the outside wall

圖11 飛行器頭罩內(nèi)壁面不同位置處溫度曲線Fig.11 Temperature distribution curve of the stagnation point at the inside wall

(a) 金屬頭罩 (b) 防熱頭罩圖12 飛行器頭罩等效應(yīng)力分布Fig.12 Equivalent stress distribution of the aircraft hood

圖13是300 s時刻兩種結(jié)構(gòu)形式的飛行器頭罩熱變形分布云圖，可見金屬頭罩發(fā)生最大變形的位置在駐點處，最大變形為4.2 mm。這是由于飛行器頭罩根部約束導(dǎo)致結(jié)構(gòu)受熱后發(fā)生膨脹。防熱頭罩前緣區(qū)域的熱變形也是指向結(jié)構(gòu)外部的，外部駐點處的最大變形為0.52 mm，內(nèi)部駐點處的變形為0.49 mm，遠(yuǎn)小于金屬頭罩的熱變形。

5 結(jié) 論

本文對流動傳熱耦合有限體積法和熱力耦合有限元法相結(jié)合的多物理場耦合計算方法進行了探索，并應(yīng)用該多場耦合計算方法對圓管繞流試驗進行仿真分析。結(jié)論如下：

1) 本文所提出的多場耦合計算方法可用于對圓管繞流試驗進行仿真分析，得到的圓管壁面壓力和熱流值與試驗結(jié)果相符，驗證了本文方法的準(zhǔn)確性和可靠性。

2) 經(jīng)過對金屬頭罩和防熱頭罩進行多場耦合仿真分析，得到了飛行器頭罩的溫度場、應(yīng)力場以及熱變形情況。由仿真數(shù)據(jù)可知，金屬頭罩內(nèi)壁面溫度比防熱頭罩內(nèi)壁面溫度高出約900 K，金屬頭罩內(nèi)壁駐點處的應(yīng)力比防熱頭罩內(nèi)壁駐點處的應(yīng)力高出529 MPa，金屬頭罩最大熱變形為4.2 mm，防熱頭罩最大熱變形為0.52 mm，該結(jié)果可以為飛行器頭罩的防熱及結(jié)構(gòu)設(shè)計提供數(shù)據(jù)支持。