沈恩楠, 陸志良, 郭同慶, 周 迪

(南京航空航天大學 航空學院, 江蘇 南京 210016)

0 引 言

高溫熱密封結構是各種航天運載火箭及空天飛行器的關鍵零部件之一[1]，比如在美國的Apollo 航天器[2]與大力神洲際導彈[3]就大量采用了此類結構。目前，國內(nèi)外已有較多的研究關注了飛行器外部高超聲速流場與結構溫度場的耦合問題。當前最常用的兩類方法為耦合方法[4-7]與同步方法[8-9]。耦合方法與同步計算方法的區(qū)別在于：耦合方法一般忽略特征時間較小一方隨時間變化的細節(jié)，具有較高的計算效率；而同步計算方法考慮連續(xù)變化的過程，理論上比耦合方法更符合實際物理過程，對流場的非定常效應模擬更加精準。耦合方法的每一個時間步需要在邊界處交換物理信息。因而當計算流場(結構溫度場)的時候，需要“凍結”結構溫度場(流場)。同步方法則通過構建統(tǒng)一的控制方程，避免了耦合過程中邊界條件的反復迭代，實現(xiàn)了流場和結構溫度場同步求解。以上研究只包含高超聲速流場與結構溫度場，未涉及空腔結構。實際上，高溫熱密封結構在外部高溫氣流的作用下會受熱導致溫度不均，進而引起密封結構內(nèi)氣體的流動。

目前數(shù)值計算空腔內(nèi)氣體熱對流的方法主要分為兩種：一種是求解基于Bossinesq假設[10-11]的不可壓Navier-Stokes (N-S)方程，另一種是通過引入預處理矩陣[12-16]求解可壓縮N-S方程。當冷熱源溫差低于30 K的時候，采用Bossinesq假設計算得到的結果與真實情況相符[11]。然而，在高溫差情況下，采用求解不可壓方程的方法將帶來數(shù)值計算上的困難[12]。Paillere等人[13]采用預處理方法對一個溫差比較小的二維封閉空間熱對流問題進行模擬，結果表明，采用該方法得到的結果與求解不可壓方程得到的結果基本一致，但是對于溫差較大的情形，結果差異較為明顯。Weiss和Simth[14]采用預處理方法對一個二維同心圓環(huán)熱對流問題進行了模擬，結果表明采用該方法能夠精確高效地得到流場特性。

綜合已有文獻來看，目前國內(nèi)外對同時考慮高超聲速氣動環(huán)境、結構傳熱以及空腔內(nèi)低速流動的瞬態(tài)同步數(shù)值模擬研究仍然較少。鑒于此，本文基于已發(fā)展的高超聲速外流場/結構溫度場同步數(shù)值模擬方法，進一步將其拓展至能夠模擬包含空腔熱對流的多區(qū)域耦合問題。通過在相鄰場交界面引入虛擬單元實現(xiàn)不同計算區(qū)域之間的物理量交換。通過與文獻對比，驗證了方法在求解單獨氣動熱/結構熱傳導問題中的有效性，同時也驗證了引入預處理矩陣方法在空腔自然對流問題中的精度。分別對封閉和帶開孔的兩種高超聲速運動圓環(huán)進行多區(qū)域同步數(shù)值模擬，最終對計算得到的空腔內(nèi)流動、氣體輸運特性以及結構溫度分布隨時間變化歷程進行了分析研究。

1 控制方程與離散方法

1.1 控制方程

內(nèi)部空腔流動從無到有的產(chǎn)生過程是內(nèi)壁面溫度連續(xù)變化導致的。隨著內(nèi)壁面溫度升高，空腔內(nèi)氣體流動速度加快，輸運特性發(fā)生變化，特征時間尺度也發(fā)生變化。同時為了能夠精確模擬空腔內(nèi)部流動導致的結構溫度場分布變化，本文采用同步計算方法。流場(包含外流場與空腔內(nèi)流場)與結構溫度場控制方程可以寫成直角坐標系下的統(tǒng)一形式：

(1)

式中，W為守恒變量。f、q為對流通量項。R、T為黏性通量項，Q為源項，這里特指空腔內(nèi)熱對流計算時出現(xiàn)的重力項。各量定義如下：

(2)

其中，p為氣體壓強，ρ為氣體密度，ρ0為參考氣體密度，E為單位質(zhì)量混合氣體的總能，k為氣體導熱系數(shù)，T為溫度，u、v、w為氣體速度在直角坐標系下的速度分量，τ為應力張量，g為重力加速度。對于理想氣體介質(zhì)有：

p=ρRT

(3)

(4)

氣體介質(zhì)的黏性系數(shù)與熱傳導系數(shù)由溫度決定：

(5)

其中μ0為海平面空氣黏度值，Pr為普朗特數(shù)。

對于固體介質(zhì)有：

ρsEs=ρsCsTs

(6)

其中，ρs為固體密度，Cs為固體比熱容，ks為固體導熱系數(shù)。在固體區(qū)域中，顯然有u=w=0，因而方程組(2)中的前三個方程自然滿足。

1.2 數(shù)值方法

整個計算空間包含三個區(qū)域，分別為外部高超聲速流動區(qū)域、結構熱傳導區(qū)域與空腔內(nèi)氣體熱對流區(qū)域。

對于外部高超聲速流動區(qū)域，對流通量采用AUSM+格式進行離散，可以保證計算的穩(wěn)定性與精度。黏性通量采用中心格式離散。結構熱傳導方程本質(zhì)與黏性擴散方程相同，因此也采用中心格式離散。對于空腔內(nèi)部熱對流區(qū)域，通過求解帶預處理矩陣的可壓縮N-S方程來模擬。時間離散采用雙時間步長推進的方法。氣體介質(zhì)內(nèi)迭代時間步長的確定方法為：

(7)

(8)

(9)

其中，μL和μT分別表示層流和湍流黏性系數(shù)，PrL和PrT分別表示層流和湍流普朗特數(shù)。對于氣體的熱傳導系數(shù)有：

(10)

由式(10)可以得到：

(11)

(12)

根據(jù)黏性譜半徑的定義與式(11)，確定了固體介質(zhì)熱傳導譜半徑：

(13)

采用熱傳導譜半徑的方法確定當?shù)貢r間步長，相當于在固體介質(zhì)中的定義了熱傳導的“黏性”，能夠保證氣體介質(zhì)和固體介質(zhì)在計算熱傳導的方法上保持一致。氣體和固體介質(zhì)邊界信息在每個虛擬迭代時間步都保持同步更新，不存在耦合時間，無需凍結流場和溫度場。

在進行同步計算之前，需要規(guī)定全部區(qū)域物理場的初始條件與邊界條件。對于外部高超聲速流動區(qū)域，采用定常方法計算得到的結果作為初始流場，將等溫壁作為交界面的邊界條件；對于結構溫度場，規(guī)定一個均勻的溫度場作為初始條件，同時規(guī)定交界面處溫度相等、熱流密度連續(xù)；空腔內(nèi)區(qū)域初場速度為0，初始溫度場與結構溫度場相同，將等溫壁作為交界面的邊界條件。

由于氣體介質(zhì)和固體介質(zhì)在邊界處的物理性質(zhì)相差較大，本文采用引入虛擬單元的方法來高效處理固體域和流體域之間的交界面。虛擬單元物理屬性與對應的真實單元物理屬性相同，單元內(nèi)的溫度采用插值方法求得：

T0=2Tw-T1

(14)

式中T0、T1分別為虛擬單元與第一層單元的溫度，Tw為交界面溫度。交界面的溫度Tw根據(jù)熱流密度連續(xù)公式確定：

(15)

式中df和ds分別為氣體和固體介質(zhì)第一層單元中心到交界面之間的距離。Tf為氣體介質(zhì)第一層單元溫度，Ts為固體介質(zhì)第一層單元溫度。ks和kf分別為固體和氣體介質(zhì)的熱傳導系數(shù)。在交界面處，氣體介質(zhì)還需滿足物面法向壓力梯度為0以及無滑移邊界條件。

2 方法驗證

2.1 高超聲速流場/結構溫度場同步計算方法

模型為一厚度為12.7 mm、外徑為76.2 mm的不銹鋼半圓環(huán)。流場參數(shù)為：Ma∞=6.47，α=0°，T∞=241.5 K，Tw=294.4 K，p∞=701.8 Pa。流場網(wǎng)格包含12 000個單元，半圓環(huán)網(wǎng)格包含3000個單元，在交界面沿周向有著相同的分布規(guī)律。流體域與固體域內(nèi)法向第一層網(wǎng)格高度均取為1.016×10-5m。圓環(huán)的內(nèi)壁采用絕熱壁邊界條件。

在進行數(shù)值計算之前，需要對初始時刻的流場進行計算。本文計算得到的駐點熱流密度為494.030 kW/m2，與文獻[5]符合較好。

表面壓強分布與文獻[17]對比如圖1所示，結果符合良好。最終同步計算方法計算得到的駐點溫度為390.227 K (702.408 R),文獻[5]計算得到的駐點溫度為388.889 K (約700 R)，計算結果與文獻符合較好。

圖1 同步計算方法計算得到的表面壓強分布與試驗值對比

2.2 熱對流預處理方法

通過與文獻對比溫度在半高度位置沿水平方向分布(圖2)可知，引入預處理矩陣的可壓縮N-S方程得到的數(shù)值模擬結果與文獻符合很好。圖3、圖4給出了不同瑞利數(shù)溫度云圖與文獻[10]對比。

(a)Ra=1×103

(b)Ra=1×105

圖3 Ra=1×103溫度云圖對比

圖4 Ra=1×105溫度云圖對比

3 多區(qū)域同步數(shù)值模擬方法應用

3.1 計算模型

模型為一厚度為3.81 mm，外徑為76.2 mm的不銹鋼圓環(huán)。流場參數(shù)為Ma∞=7，α=0°，T∞=216.65 K，Tw=294.4 K，p∞=5.529×103Pa，Tin=294.4 K，pin=1.01325×105Pa。

Tin為空腔內(nèi)氣體初始溫度，pin為空腔內(nèi)氣體初始壓強。圖5為計算網(wǎng)格，其中外部高超聲速流動區(qū)域包含48 000個單元，結構熱傳導區(qū)域包含12 000個單元，空腔內(nèi)部熱對流區(qū)域包含15 116個單元，邊界處網(wǎng)格沿周向分布一致。

(a)整體 (b)局部

3.2 氣動熱/結構熱傳導引起封閉空腔內(nèi)氣體流動數(shù)值模擬

在高超聲速氣動加熱的作用下，外部流動前緣駐點處的溫度會急劇上升，此處的熱流密度達到最大，結構傳熱效果最明顯。隨著結構傳熱過程的進行，結構內(nèi)壁前緣最先升溫，造成空腔內(nèi)流場溫度分布不均勻，從而產(chǎn)生熱對流現(xiàn)象。

圖6、圖7顯示，當t=5 s時，由于空腔前后緣內(nèi)壁溫度均有增加，高溫氣流上升，低溫氣流下降，氣流分別在前后緣處產(chǎn)生了旋渦；當t=15 s時，腔內(nèi)空氣的流動逐漸顯著，流速增加，氣體的輸運特性也越發(fā)明顯，前緣高溫氣流不斷上升，后緣處旋渦消失，旋渦中心距離前緣較近；在t=25 s后的流動中，旋渦中心向空腔中心移動，旋渦形狀與渦核位置逐漸穩(wěn)定。

圖8為不同時刻是否考慮空腔內(nèi)氣體流動的結構溫度差，即將不考慮空腔氣體流動情形下的結構單元溫度減去考慮空腔氣體流動情形下相應單元的溫度。當t=5 s時，最大差別點出現(xiàn)在近空腔前緣附近，原因是對于后者，該處結構與空腔內(nèi)氣體溫差最大，結構向空腔氣體傳熱最明顯。隨著時間推進，由重力驅(qū)動的高溫氣體將能量從駐點向其他地方輸送，高溫氣體上升，對上半圓環(huán)進行加熱，上半圓環(huán)溫度升高，結構溫差減小；而低溫氣體下降，在渦流作用下經(jīng)過后半圓環(huán)向下運動，吸收下半圓環(huán)的熱量，下半圓環(huán)溫度降低，導致下半圓環(huán)結構溫差增加。

(a)t=5 s

(b)t=15 s

(c)t=25 s

(d)t=35 s

(a)t=5 s

(b)t=15 s

(c)t=25 s

(d)t=35 s

(a)t=5 s

(b)t=15 s

(c)t=25 s

(d)t=35 s

3.3 氣動熱/結構熱傳導引起開孔空腔內(nèi)氣體流動數(shù)值模擬

在上文模型的基礎上對開孔空腔模型進行同步數(shù)值模擬進行研究。如圖9所示，開孔位置位于來流方向逆時針旋轉90°方向，開孔內(nèi)壁兩端點與圓心連線夾角為3.31°。初始外部流場和初始空腔內(nèi)部流場與上一節(jié)相同。

(a)整體

(b)開孔處

空腔內(nèi)流場不同時刻流線如圖10所示。通過不同時刻流線圖可以看出，在流動初期，由于流速較快，在開孔處形成旋渦，堵塞流動，外部氣體沿著右側孔壁流入腔內(nèi)。受到開孔壁面影響，氣體成一定角度流入空腔，空腔內(nèi)氣流被流入氣流分成兩個部分，形成兩個旋渦。經(jīng)過一段時間后，空腔內(nèi)氣流流速增加，與流入的氣流形成干擾?？鬃髠鹊男郎u逐漸減小，在孔右側形成了一個新的旋渦。此外，如圖11所示，開孔的存在對于結構溫度場影響較大，主要體現(xiàn)在孔附近的壁面。t=0.1 s時，外部流場駐點溫度為361.29 K，開孔邊緣溫度為508.11 K；t=0.3 s時，駐點溫度為413.01 K，開孔邊緣溫度為468.53 K；t=0.5 s時，駐點溫度為446.99 K，開孔邊緣溫度為451.86 K?？梢钥闯觯虝r間內(nèi)開孔邊緣溫度可以超過外部流場駐點溫度。由此可見，當密封結構存在開孔時會出現(xiàn)局部過熱的現(xiàn)象，其附近結構瞬時溫度甚至會超過駐點溫度，可能造成結構破壞。

(a)t=0.1 s

(b)t=0.3 s

(c)t=0.5 s

(a)t=0.1 s

(b)t=0.3 s

(c)t=0.5 s

4 結 論

本文在已發(fā)展的高超聲速外流場/結構溫度場同步計算方法基礎上，進一步發(fā)展了考慮了空腔內(nèi)低速流動的同步數(shù)值模擬方法。通過對不同時刻空腔內(nèi)流動特征和結構溫度場分布進行對比分析，研究了結構溫度場不均勻分布引起的空腔內(nèi)對流流動反作用于結構溫度場的影響，以及開孔對結構溫度分布的影響?？偨Y得出以下結論：

1)封閉空腔內(nèi)氣體由于結構溫度分布不均勻產(chǎn)生了流動。在空腔內(nèi)氣體流動的影響下，封閉圓環(huán)前緣溫度在35 s內(nèi)最多下降0.8%左右。

2)開孔空腔的開孔處邊緣溫度在0.5 s內(nèi)能夠超過外流駐點溫度。

由此可以看出，高超聲速飛行器封閉結構內(nèi)空腔氣體流動對結構溫度場影響較小，而結構間隙對結構溫度場的影響不容忽視。未來可以進一步深入研究縫隙大小對結構溫度場影響。所發(fā)展的方法為未來薄壁、帶空腔的高超聲速飛行器空腔內(nèi)部流動數(shù)值模擬，以及高溫熱密封結構與熱防護設計提供一種有效手段。