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

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

0 引 言

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

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

1 數(shù)值計(jì)算方法

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

(1)

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

(2)

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

(3)

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

(4)

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

結(jié)構(gòu)溫度變化會(huì)引起結(jié)構(gòu)體積發(fā)生變化，同時(shí)，內(nèi)部溫度場(chǎng)分布不均勻會(huì)導(dǎo)致結(jié)構(gòu)產(chǎn)生的熱變形不協(xié)調(diào)，因此會(huì)產(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)克符號(hào)。

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

(7)

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

2 多物理場(chǎng)耦合策略

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

Tfluid=Tsolid

(8)

(9)

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

3 多場(chǎng)耦合方法驗(yàn)證算例

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

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

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

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

Relocal=ρUΔn0/μ

(10)

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

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

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

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

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

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

圖2 計(jì)算得到的流場(chǎng)密度與試驗(yàn)結(jié)果對(duì)比Fig.2 Comparison of calculated flow field density and test result

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.2.1 流動(dòng)傳熱耦合分析

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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