吳王浩，段 旭，張 鑫，陳 丹，徐振東

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

0 引 言

隨著飛行器飛行速度的不斷提高，高馬赫數(shù)下的氣動(dòng)加熱逐漸成為一個(gè)不得不考慮的問題。在高馬赫數(shù)繞流的邊界層內(nèi)，空氣的大部分動(dòng)能都會(huì)向內(nèi)能轉(zhuǎn)化，從而產(chǎn)生氣動(dòng)加熱效應(yīng)，導(dǎo)致飛行器表面的空氣溫度非常高?？諝庠跇O端的高溫條件下會(huì)發(fā)生化學(xué)反應(yīng)，隨著溫度的進(jìn)一步升高，還會(huì)發(fā)生離解甚至電離。長時(shí)間的氣動(dòng)加熱會(huì)使結(jié)構(gòu)內(nèi)部溫度發(fā)生變化，引起溫度分布不均勻，從而會(huì)產(chǎn)生結(jié)構(gòu)的變形，甚至發(fā)生斷裂。

對(duì)于一些簡單外形，一般是在其表面鋪設(shè)燒蝕材料來進(jìn)行防熱處理，這種方法簡單而且有較高的裕度，對(duì)于氣動(dòng)熱預(yù)測的精度要求不高。隨著飛行器外形變得越來越復(fù)雜，飛行環(huán)境越來越苛刻，對(duì)飛行器熱防護(hù)提出的要求也越來越高。尤其是近幾年高馬赫數(shù)飛行器的不斷發(fā)展進(jìn)步，和以前的再入飛行器不同，現(xiàn)在的高馬赫數(shù)飛行器設(shè)計(jì)要滿足長時(shí)間的臨近空間飛行，可重復(fù)使用等要求，這就導(dǎo)致熱防護(hù)技術(shù)面臨新的問題，長時(shí)間的飛行所產(chǎn)生的熱累積效應(yīng)非常嚴(yán)重，所以對(duì)氣動(dòng)加熱的預(yù)測精度有了更高的要求。

根據(jù)經(jīng)驗(yàn)，在高馬赫數(shù)飛行狀態(tài)下，一般飛行器受氣動(dòng)加熱最嚴(yán)重的位置往往位于前緣駐點(diǎn)處，所以駐點(diǎn)處氣動(dòng)加熱效應(yīng)需要著重關(guān)注。受限于計(jì)算能力，早期對(duì)于飛行器整體熱流求解困難，一般使用經(jīng)驗(yàn)或者半經(jīng)驗(yàn)的工程估算方法計(jì)算駐點(diǎn)處加熱狀況。

傳統(tǒng)方法將氣動(dòng)加熱與結(jié)構(gòu)傳熱分兩步進(jìn)行計(jì)算，不考慮結(jié)構(gòu)溫度變化對(duì)流場的影響。使用氣動(dòng)-傳熱耦合計(jì)算考慮了結(jié)構(gòu)傳熱與氣動(dòng)加熱之間的相互作用，更接近實(shí)際情況。為了準(zhǔn)確地預(yù)測氣動(dòng)加熱對(duì)飛行器結(jié)構(gòu)所產(chǎn)生的影響，需要考慮氣動(dòng)加熱和結(jié)構(gòu)傳熱耦合分析。

對(duì)于高馬赫數(shù)鈍頭體的氣動(dòng)加熱問題，Wieting等進(jìn)行了實(shí)驗(yàn)和仿真計(jì)算；Dechaumphai等研究了氣動(dòng)、傳熱和結(jié)構(gòu)應(yīng)力的耦合計(jì)算；黃唐等完成了二維的流場、熱、結(jié)構(gòu)一體化數(shù)值模擬；張智超等對(duì)高馬赫數(shù)氣動(dòng)熱數(shù)值計(jì)算中壁面網(wǎng)格劃分尺度進(jìn)行了討論；夏剛等對(duì)迭代耦合和單向耦合求解結(jié)構(gòu)傳熱進(jìn)行了對(duì)比。

本文求解高馬赫數(shù)鈍頭體流場，驗(yàn)證了流場求解和高馬赫數(shù)下熱邊界確定的準(zhǔn)確性；使用迭代耦合實(shí)現(xiàn)了氣動(dòng)-傳熱的一體化計(jì)算；同時(shí)對(duì)比了迭代耦合計(jì)算和單向求解計(jì)算得到的結(jié)果之間的差異。

1 計(jì)算方法

1.1 流動(dòng)控制方程

使用CFD 數(shù)值求解方法進(jìn)行高馬赫數(shù)飛行狀態(tài)下的氣動(dòng)加熱模擬，使用的基本方程為N-S方程，積分形式的N-S方程如下：

式中：=(，，，，)為守恒向量；、(，，)和分別是密度、直角坐標(biāo)系下的速度分量和單位質(zhì)量氣體的總能量；?是某一固定區(qū)域的邊界；是邊界的外法向量；矢通量可分解成對(duì)流矢通量和粘性矢通量兩部分。

對(duì)于完全氣體，其狀態(tài)方程為

式中：為空氣壓強(qiáng)；為氣體密度；為氣體常數(shù)；為氣體溫度；為氣體焓值；為氣體定壓比熱。

1.2 熱傳導(dǎo)方程

熱傳導(dǎo)是由于溫度梯度導(dǎo)致的能量轉(zhuǎn)移而產(chǎn)生的。直角坐標(biāo)系下，基于各向同性假設(shè)的二維瞬態(tài)熱傳導(dǎo)方程可以寫成以下形式：

式中：為固體密度；為固體的定壓比熱；為固體溫度；為固體的熱傳導(dǎo)系數(shù)。

本文由于結(jié)構(gòu)溫度變化不大，溫度變化對(duì)熱傳導(dǎo)系數(shù)的影響很小，所以不考慮溫度變化對(duì)材料特性所產(chǎn)生的影響。

1.3 耦合計(jì)算方法

通過迭代耦合的方法實(shí)現(xiàn)了氣動(dòng)-傳熱的耦合計(jì)算。圖1為迭代耦合數(shù)據(jù)傳遞示意圖，其中：為瞬時(shí)表面熱流分布；為瞬時(shí)表面溫度分布。在流場和結(jié)構(gòu)傳熱的瞬態(tài)計(jì)算的每個(gè)時(shí)間步上，由于時(shí)間步長很小，所以對(duì)于流場來說可以認(rèn)為每個(gè)時(shí)間步上結(jié)構(gòu)表面的溫度是不變的，對(duì)于結(jié)構(gòu)傳熱來說可以認(rèn)為每個(gè)時(shí)間步上的熱流邊界也是固定不變的。

圖1 迭代耦合數(shù)據(jù)傳遞示意圖Fig.1 Data transfer of iterative coupling method

迭代耦合的步驟如下：

1）假定物體初始時(shí)刻結(jié)構(gòu)溫度是一個(gè)固定值，在計(jì)算流場時(shí)將物體表面看作等溫壁，求解穩(wěn)態(tài)流場得到一個(gè)結(jié)構(gòu)表面的初始熱流邊界。

2）將流場與結(jié)構(gòu)的交界面作為耦合邊界，在邊界上進(jìn)行流場與結(jié)構(gòu)之間數(shù)據(jù)的傳遞，進(jìn)行非定常的模擬。

3）在時(shí)間步Δt中，假定結(jié)構(gòu)溫度不變，可以求得該時(shí)間步上流場的結(jié)構(gòu)表面熱流。

4）將流場求得的熱流傳遞到結(jié)構(gòu)傳熱瞬態(tài)計(jì)算中作為該時(shí)間步的熱流邊界，假定熱流邊界在該時(shí)間步內(nèi)固定不變，進(jìn)行傳熱計(jì)算，得到下一個(gè)時(shí)間步的結(jié)構(gòu)初始溫度。

5）將結(jié)構(gòu)溫度傳遞給流場計(jì)算中，重復(fù)步驟3～4，如此迭代。

2 計(jì)算結(jié)果及分析

2.1 計(jì)算模型

1987年，Wieting 等為了研究氣動(dòng)-熱-結(jié)構(gòu)耦合效應(yīng)，對(duì)圓柱模型分別進(jìn)行了實(shí)驗(yàn)和仿真數(shù)值計(jì)算。該實(shí)驗(yàn)被后來多次用于驗(yàn)證耦合傳熱仿真的可靠性。

實(shí)驗(yàn)對(duì)象為一長直圓管，對(duì)該實(shí)驗(yàn)進(jìn)行簡化后，使用二維圓管進(jìn)行仿真數(shù)值計(jì)算，圓管的外形及網(wǎng)格劃分如圖2所示。由于流場是上下對(duì)稱的，所以只取上半部分流場進(jìn)行計(jì)算。本文只關(guān)心圓管迎風(fēng)面處的氣動(dòng)加熱情況，因此對(duì)背風(fēng)面處的流動(dòng)狀況不進(jìn)行討論，且圓管背風(fēng)面對(duì)計(jì)算結(jié)果沒有影響，所以只取前1/4 圓管進(jìn)行計(jì)算。圓管的內(nèi)徑為24.5 mm，外徑為38.1 mm。網(wǎng)格中，為了準(zhǔn)確計(jì)算壁面處的熱流，對(duì)靠近壁面處的網(wǎng)格進(jìn)行了加密處理，使壁面處的網(wǎng)格雷諾數(shù)接近于1，流場中靠近壁面的網(wǎng)格第一層高度為3.86×10m。

圖2 流場和結(jié)構(gòu)網(wǎng)格Fig.2 The grid of flow field and structure

2.2 計(jì)算條件

二維圓柱繞流計(jì)算中，來流溫度=241.5 K，馬赫數(shù)=6.47，雷諾數(shù)=1.31×10m；壓力=648.1 Pa，使用層流進(jìn)行計(jì)算。

對(duì)于結(jié)構(gòu)部分，實(shí)驗(yàn)中材料為321 不銹鋼，其熱力學(xué)參數(shù)為：密度=7 920 kg/m，比熱=500 J/(kg·K)，熱傳導(dǎo)系數(shù)=16.3 W/(m·K)，固體結(jié)構(gòu)的初始溫度為294.5 K。

2.3 穩(wěn)態(tài)計(jì)算

圖3為穩(wěn)態(tài)流場中計(jì)算得到的熱流以駐點(diǎn)熱流值為基準(zhǔn)進(jìn)行歸一化后的曲線與實(shí)驗(yàn)數(shù)據(jù)之間的對(duì)比。從圖中可以看出，熱流繞圓柱的分布趨勢(shì)與實(shí)驗(yàn)數(shù)據(jù)吻合很好，但是駐點(diǎn)熱流與實(shí)驗(yàn)值有一定偏差。本文得到的仿真駐點(diǎn)熱流為5.057×10W/m，文獻(xiàn)［7］中仿真計(jì)算得到的駐點(diǎn)熱流為4.82×10W/m，文獻(xiàn)［11］計(jì)算值為5.41×10W/m。從與其他文獻(xiàn)的結(jié)果對(duì)比中可以看出，文獻(xiàn)中的數(shù)值模擬和本章的模擬均使用了層流進(jìn)行計(jì)算，本章和其他文獻(xiàn)計(jì)算結(jié)果接近，相差6.5%。

圖3 鈍頭體表面穩(wěn)態(tài)熱流比較Fig.3 The steady state heat flux comparison of the blunt body

2.4 迭代耦合

耦合傳熱得到的結(jié)果如圖4～圖8所示。由于高馬赫數(shù)氣動(dòng)加熱這類問題中，流場的特征時(shí)間遠(yuǎn)遠(yuǎn)小于結(jié)構(gòu)傳熱的特征時(shí)間，所以先將固體壁面設(shè)置為等溫壁，計(jì)算得到穩(wěn)態(tài)流場和壁面上的熱流分布，再以此為初始條件進(jìn)行流場和固體結(jié)構(gòu)之間的傳熱耦合計(jì)算。

圖4為不同時(shí)刻下鈍頭體表面熱流分布曲線。隨著時(shí)間的增加，靠近駐點(diǎn)處溫度隨著時(shí)間的推進(jìn)上升較快，導(dǎo)致駐點(diǎn)附近溫度梯度下降，從而導(dǎo)致歸一化熱流有下降的趨勢(shì)。由于圓柱后半部分溫度上升不明顯，所以50°之后的熱流曲線基本沒有產(chǎn)生變化。

圖4 不同時(shí)刻表面熱流曲線Fig.4 Surface heat flow at different times

圖5為熱流在進(jìn)行耦合傳熱后駐點(diǎn)熱流的下降趨勢(shì)，從圖中看出，在剛進(jìn)行耦合傳熱的極短時(shí)間內(nèi)，熱流有較明顯的下降趨勢(shì)。之后隨著時(shí)間的增加，駐點(diǎn)溫度不斷升高，駐點(diǎn)熱流下降趨于平緩?？梢灶A(yù)見，經(jīng)過足夠長的時(shí)間后，熱流會(huì)達(dá)到一個(gè)相對(duì)穩(wěn)定的值。

圖5 駐點(diǎn)熱流隨時(shí)間變化曲線Fig.5 Stagnation heat flow varying with time

圖6為駐點(diǎn)溫度隨時(shí)間變化的曲線，剛開始由于駐點(diǎn)熱流較大，溫度上升比較快，隨著駐點(diǎn)溫度上升，熱流下降，駐點(diǎn)溫度上升速度逐漸降低，同時(shí)熱流也逐漸降低，最后熱流和圓柱表面溫度分布都會(huì)達(dá)到一個(gè)平衡狀態(tài)。本文計(jì)算得到的2 s 時(shí)刻下駐點(diǎn)溫度為395 K，文獻(xiàn)［7］計(jì)算得到的2 s 時(shí)刻下駐點(diǎn)溫度為388.9 K，誤差僅為1.6%，說明本文耦合傳熱計(jì)算結(jié)果是可靠的。

圖6 駐點(diǎn)溫度變化曲線Fig.6 Stagnation point temperature change

圖7為圓柱駐點(diǎn)處沿半徑溫度分布隨時(shí)間的變化，在0 時(shí)刻即初始條件下，物面溫度為294.5 K。隨著傳熱的進(jìn)行，駐點(diǎn)溫度不斷升高，且從圖中可以看出溫度沿半徑逐漸向內(nèi)傳播的過程。

圖7 圓柱駐點(diǎn)處沿半徑溫度變化Fig.7 Temperature variation along the radius at the stagnation point

圖8為不同時(shí)刻下圓柱表面溫度分布，表面初始溫度為等溫壁294.5 K，在進(jìn)行耦合傳熱計(jì)算后，物體不斷受到流場的氣動(dòng)加熱效應(yīng)，表面溫度不斷上升。由于駐點(diǎn)處受到的氣動(dòng)加熱最嚴(yán)重，所以駐點(diǎn)處溫度上升最快，但50°之后的溫度曲線變化比較緩慢。這是由于表面角度的增加使得物面受到的氣流加熱效應(yīng)大幅下降，所以溫度上升緩慢。

圖8 不同時(shí)刻表面溫度分布Fig.8 Surface temperature distribution at different times

2.5 單向求解

通過分別求解N-S 方程和結(jié)構(gòu)熱傳導(dǎo)方程，實(shí)現(xiàn)了氣動(dòng)加熱-傳熱的單向耦合計(jì)算。由于高馬赫數(shù)飛行下壁面附近流場溫度很高，短時(shí)間內(nèi)結(jié)構(gòu)表面溫度的上升對(duì)于壁面處溫度梯度影響不大，所以可以假定在短時(shí)間內(nèi)結(jié)構(gòu)表面熱流即為以等溫壁為條件的穩(wěn)態(tài)流場下求得的表面熱流。然后再將該熱流作為結(jié)構(gòu)傳熱計(jì)算的熱流邊界條件進(jìn)行傳熱分析，不考慮結(jié)構(gòu)溫度上升對(duì)流場的反作用。

實(shí)際情況下，隨著氣動(dòng)加熱和結(jié)構(gòu)傳熱的進(jìn)行，結(jié)構(gòu)表面溫度會(huì)上升，從而導(dǎo)致結(jié)構(gòu)表面熱流降低。從圖4可以看出，隨著時(shí)間的增加，表面溫度上升還是會(huì)導(dǎo)致表面熱流下降，從而會(huì)影響結(jié)構(gòu)傳熱。由于使用單向耦合不考慮表面熱流隨時(shí)間的下降，所以計(jì)算得到的結(jié)構(gòu)溫度相較于實(shí)際情況和耦合傳熱會(huì)偏高。

圖9和圖10分別為2 s 和10 s 時(shí)單向求解結(jié)構(gòu)瞬態(tài)傳熱得到的溫度場，從圖中可以看出，隨著傳熱的進(jìn)行，溫度不斷向結(jié)構(gòu)內(nèi)部進(jìn)行傳導(dǎo)，同時(shí)結(jié)構(gòu)表面溫度也隨時(shí)間不斷升高。

圖9 2 s時(shí)結(jié)構(gòu)溫度分布Fig.9 Structure temperature distribution at 2 s

圖10 10 s時(shí)結(jié)構(gòu)溫度分布Fig.10 Structure temperature distribution at 10 s

2.6 計(jì)算結(jié)果對(duì)比

圖11為使用單向求解和迭代耦合分別計(jì)算得到的駐點(diǎn)熱流變化曲線。單向求解進(jìn)行到2 s時(shí)，駐點(diǎn)溫度為396.7 K，使用迭代耦合計(jì)算得到2 s 時(shí)的駐點(diǎn)溫度為395 K，兩者差距僅為0.43%；在10 s時(shí)，迭代耦合和單向求解兩種方法計(jì)算得到的駐點(diǎn)溫度分別為536.12 K 和517.2 K，差距為3.66%。這說明兩種計(jì)算方法在傳熱進(jìn)行后的較短時(shí)間內(nèi)，計(jì)算得到的結(jié)構(gòu)溫度之間的差距非常小，隨著加熱時(shí)間的增長，兩者計(jì)算溫度差距逐漸增大，迭代耦合計(jì)算得到的溫度低于單向求解的溫度。相比于迭代耦合傳熱計(jì)算，使用單向求解的方法能提高計(jì)算的效率，并且在短時(shí)間內(nèi)能保持較高的精度，但是隨著傳熱時(shí)間的增加，單向求解的誤差會(huì)逐漸增加，而迭代耦合計(jì)算的結(jié)構(gòu)更加接近實(shí)際值。

圖11 迭代耦合和單向求解不同時(shí)刻駐點(diǎn)溫度比較Fig.11 The Stagnation point temperatures of iterative coupling and one-way solving methods at different times

3 結(jié) 論

本文通過CFD 數(shù)值模擬方法求解N-S方程，得到了鈍頭體表面穩(wěn)態(tài)熱流分布，分別用迭代耦合和單向求解的方法模擬了結(jié)構(gòu)內(nèi)部傳熱過程，并將兩者的結(jié)果進(jìn)行了對(duì)比，結(jié)論如下：

1）單向求解得到的結(jié)構(gòu)溫度要高于迭代耦合計(jì)算結(jié)果。使用迭代耦合計(jì)算方法相比于單向耦合可以更加準(zhǔn)確地預(yù)測氣動(dòng)加熱和結(jié)構(gòu)傳熱的過程，但迭代耦合需要在每一個(gè)時(shí)間步內(nèi)都對(duì)流場進(jìn)行求解，所以迭代耦合計(jì)算量會(huì)遠(yuǎn)大于單向耦合。

2）在相同的來流條件下，傳熱開始后的短時(shí)間內(nèi)，迭代耦合和單向耦合計(jì)算的結(jié)果相差很小，隨著時(shí)間的增加，誤差逐漸增大。所以在一定的誤差許可下，使用單向耦合能極大地提高計(jì)算效率，同時(shí)得到相對(duì)比較準(zhǔn)確的結(jié)果，但隨著時(shí)間加長，單向求解誤差將逐漸增大。