高馬赫數(shù)鈍頭體氣動/傳熱一體化計算方法研究

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

（上海機電工程研究所，上海 201109）

0 引 言

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

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

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

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

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

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

1 計算方法

1.1 流動控制方程

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

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

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

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

1.2 熱傳導(dǎo)方程

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

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

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

1.3 耦合計算方法

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

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

迭代耦合的步驟如下：

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

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

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

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

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

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

2.1 計算模型

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

實驗對象為一長直圓管，對該實驗進行簡化后，使用二維圓管進行仿真數(shù)值計算，圓管的外形及網(wǎng)格劃分如圖2所示。由于流場是上下對稱的，所以只取上半部分流場進行計算。本文只關(guān)心圓管迎風(fēng)面處的氣動加熱情況，因此對背風(fēng)面處的流動狀況不進行討論，且圓管背風(fēng)面對計算結(jié)果沒有影響，所以只取前1/4 圓管進行計算。圓管的內(nèi)徑為24.5 mm，外徑為38.1 mm。網(wǎng)格中，為了準確計算壁面處的熱流，對靠近壁面處的網(wǎng)格進行了加密處理，使壁面處的網(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 計算條件

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

對于結(jié)構(gòu)部分，實驗中材料為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)計算

圖3為穩(wěn)態(tài)流場中計算得到的熱流以駐點熱流值為基準進行歸一化后的曲線與實驗數(shù)據(jù)之間的對比。從圖中可以看出，熱流繞圓柱的分布趨勢與實驗數(shù)據(jù)吻合很好，但是駐點熱流與實驗值有一定偏差。本文得到的仿真駐點熱流為5.057×10W/m，文獻［7］中仿真計算得到的駐點熱流為4.82×10W/m，文獻［11］計算值為5.41×10W/m。從與其他文獻的結(jié)果對比中可以看出，文獻中的數(shù)值模擬和本章的模擬均使用了層流進行計算，本章和其他文獻計算結(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ù)氣動加熱這類問題中，流場的特征時間遠遠小于結(jié)構(gòu)傳熱的特征時間，所以先將固體壁面設(shè)置為等溫壁，計算得到穩(wěn)態(tài)流場和壁面上的熱流分布，再以此為初始條件進行流場和固體結(jié)構(gòu)之間的傳熱耦合計算。

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

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

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

圖5 駐點熱流隨時間變化曲線Fig.5 Stagnation heat flow varying with time

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

圖6 駐點溫度變化曲線Fig.6 Stagnation point temperature change

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

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

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

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

2.5 單向求解

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

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

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

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

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

2.6 計算結(jié)果對比

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

圖11 迭代耦合和單向求解不同時刻駐點溫度比較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é)果進行了對比，結(jié)論如下：

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

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