王 強，徐 濤，姚永濤

（1.中北大學 能源動力工程學院，太原 030051；2.特種環(huán)境復合材料技術國家級重點實驗室，哈爾濱 150001）

引言

高超聲速飛行器是21 世紀航空航天的一個主要發(fā)展方向，其總體技術涵蓋氣動力與氣動熱數值模擬、一體化設計、多學科優(yōu)化設計、熱防護與熱管理技術等方面[1].氣動熱效應顯著是高超聲速飛行器的重要特征，準確預測飛行器的氣動力與熱載荷是進行氣動設計、熱防護系統(tǒng)設計的前提.高超聲速飛行器表面通常存在縫隙，如防熱瓦之間、控制舵與機體之間；由于安裝、局部受熱不均等影響，防熱瓦之間還存在一定的階差.在飛行過程中，高速氣流會流入防熱瓦縫隙中，導致防熱瓦表面的流動特性與傳熱方式發(fā)生變化[2]，階差使防熱瓦表面的流動與換熱情況進一步復雜化.此外飛行器頭部一般為鈍體結構，在高超聲速來流作用下，局部氣動熱現(xiàn)象嚴重，準確預測該區(qū)域溫度分布對飛行器的設計至關重要.

Hinderks 等[3]采用弱耦合方法實現(xiàn)了自編程序TAU 以及商業(yè)求解器ANSYS 的耦合，進行了縫隙流動、防熱瓦傳熱以及熱變形的氣-熱-彈耦合仿真分析，其研究揭示了縫隙內的旋渦結構，同時由于防熱瓦熱變形的緣故，縫隙幾何形狀會發(fā)生變化，進而影響縫隙內壁面的熱流密度分布.應用多場耦合技術，沈淳等[4]研究了縫隙-腔體密封結構在高速氣流沖擊下的整體流動、傳熱特征；殷超等[5]采用FLUENT 商業(yè)求解器對高超聲速飛行器縫隙流動傳熱問題開展了數值仿真，研究了不同狀態(tài)參數對縫隙流動與傳熱的影響.邱波等[2]進行了高超聲速飛行器表面橫縫旋渦結構及氣動熱環(huán)境數值模擬，對縫隙表面熱流分布、縫隙內流動狀況等進行了較系統(tǒng)的研究.

聶濤等[6]采用有限體積法與有限元法分別求解了非定常Navier-Stokes 方程與非穩(wěn)態(tài)導熱方程，基于準穩(wěn)態(tài)假設對固體結構問題進行了分析，實現(xiàn)了高超聲速飛行器前緣流固耦合的計算.基于熱流計算的可靠性問題，李邦明等[7]研究了高超聲速飛行器前駐點熱流數值模擬的物理準則.張昊元等[8]對一種開縫前緣的簡化模型進行了數值仿真，研究了縫隙誘導形成的三維旋渦的空間分布特征和旋渦運動對物面氣動加熱的影響規(guī)律，發(fā)現(xiàn)縫隙內主旋渦的再附導致了側壁存在“非常規(guī)”的高熱流區(qū).

針對高超聲速流動與換熱問題，本課題組基于有限差分法開發(fā)了仿真求解器，并將該求解器分別用于高超聲速臺階、縫隙流動以及非定常鈍體繞流與換熱問題的仿真求解，分析其流動與傳熱特點，并對求解器的計算能力進行了驗證.

1 數學物理模型

任意曲線坐標系下，引入預處理技術的守恒變量形式的無量綱Reynolds 平均Navier-Stokes（Reynolds averaged Navier-Stokes，RANS）方程為

其中Q=J·(ρ ρuρvρw e)，ρ為密度，u，v，w為速度矢量在直角坐標系下的三個分量；E，F(xiàn)，G和Ev，F(xiàn)v，Gv分別為曲線坐標系的無黏通量與黏性通量；t為時間；ξ，η，ζ為自然曲線坐標系的坐標分量.方程中的無黏通量采用AUSM + - up 格式[9]離散，該差分格式為AUSM 類差分格式最新發(fā)展的類型，該格式在構造過程中，考慮了低Mach 數效應的影響，可以應用于全速域流動的仿真計算，具有擴展性好的優(yōu)點；黏性通量采用中心差分格式離散.采用q-ω低Reynolds 數二方程模型[10]對RANS 方程進行封閉，該模型本質上是一種低Reynolds 數湍流模型，在壁面網格滿足要求的條件下，不需要壁面函數，對邊界層流動的預測精度較高，同時該模型對高超聲速流動預測也較令人滿意.采用預處理技術提高求解器對不同流域流動問題的求解能力[11]，離散后的代數方程組采用目前在計算流體力學領域廣泛應用的LU-SGS 隱式方法[12]求解.

固體導熱微分方程為

其中C為固體比熱容，T為溫度，λ為固體導熱系數.由于導熱的擴散性質，在計算中采用中心差分離散，并采用Douglas-Rachford（D-R）交替隱式迭代法[13]求解離散后的固體導熱代數方程組.

對氣熱耦合仿真需要保證交接面處流固兩側壁溫Tw相同，熱流密度qw連續(xù)，即涉及到熱流密度計算，本文采用直接耦合方法實現(xiàn)流體與固體交界面處的數據傳遞，該方法避免了熱流密度qw的計算，通過流固兩側單元的溫度、兩側網格單元距離交接面的距離Δn計算壁面溫度，易于數據傳遞的實施，流固耦合交接面上溫度計算如下：

非定常氣熱耦合計算包括流動控制方程以及固體導熱方程的非定常計算，必須要考慮到流場以及固體溫度場的時間同步推進.在計算中，對流動控制方程，采用雙時間步法[14]實現(xiàn)非定常隱式求解；對非定常固體導熱方程，將D-R 隱式迭代法中的時間替換成物理時間，進行顯式推進；考慮到氣動加熱問題的強耦合屬性，為了保證計算的精度，在計算中對兩個物理場采用統(tǒng)一的物理時間步長.

在定常計算中，當參數的平均殘差小于10-4或殘差趨于水平后，即終止迭代；對非定常計算，當迭代物理時間達到總的物理時間后計算終止.

2 高超聲速流動與換熱數值仿真

2.1 高超聲速后臺階流動與換熱

采用與Grotowsky 與Ballmann 相同的后臺階算例[15]進行后臺階流動與傳熱的驗證，該算例試驗數據來自于Jesson 等[16]的研究，該算例文獻中的模型如圖1 所示.

圖1 后臺階算例模型尺寸Fig.1 Geometry of the back-step model

來流條件為：Mach 數7.898，溫度122.0 K，壓強613 Pa，Reynolds 數3.716 × 106，攻角-15°.臺階高度H為6 mm、入口段長度L為50 mm、出口位于臺階下游120 mm 處，計算中設置為超音速出口.計算中壁面溫度設置為300 K.出于對比目的，壁面第一層網格單元距離為1.0 × 10-5m，與文獻[15]保持一致.網格單元數為26 496 個，與文獻接近（25 836 個），計算網格如圖2 所示.

圖2 高超聲速后臺階流動計算網格Fig.2 The computational mesh for hypersonic flow over backward facing step

圖3 給出了預測的密度等值線云圖，圖中參數用參考值0.012 51 kg/m3進行無量綱化處理.以-15°攻角流過入口段，由于流動面積突變，在入口段最前方產生一道斜激波，氣體流過斜激波，由于流動突擴，后臺階位置拐點處形成膨脹波系，同時流動分流在拐角形成漩渦流動，氣流在臺階下游，x≈ 82 mm 處（文獻給出的位置為81.87 mm）再附，并在該位置處產生一道再附激波.

圖3 計算得到的密度等值線圖Fig.3 Predicted density contours

氣動力與氣動熱是高超聲速飛行器設計的重要參數，圖4、圖5 分別給出了本程序預測的壁面壓力與壁面熱流參數，并與文獻中的計算結果、試驗結果進行了對比.在后臺階拐點之后，出現(xiàn)膨脹波系，加速流動導致局部的壓強、密度等參數驟降，壁面壓力值降低約70%，直到再附點附近，強壓縮流動導致局部流動參數急劇提升，與再附激波前壓強值相比，約增加2 倍.與壓強變化相比，壁面熱流密度變化具有相同的趨勢，但是變化的程度更大，從試驗結果來看，臺階拐點前后、再附激波前后等區(qū)域熱流密度值變化均在10 倍以上.從對比結果來看，在后臺階角點之前的入口段、分離流動再附之后的區(qū)域，本程序計算的壁面壓力與試驗結果、文獻計算結果吻合較好，誤差約9.6%，在10%以內；但是臺階角點后、再附點之前的區(qū)域（50 mm ＜x＜ 80 mm）內，文獻計算結果和本程序結果均與試驗值存在較大的誤差，該區(qū)域的預測也是目前湍流模型研究的一個難點問題，但是程序計算結果與文獻結果相比，誤差在10%以內.從熱流分布結果來看，本程序計算的熱流較文獻計算結果更接近于試驗結果，在50 mm ＜x＜ 80 mm 區(qū)域內，預測的熱流與試驗值吻合很好，誤差在5%之內，較大的誤差出現(xiàn)在再附激波以后的區(qū)域，該區(qū)域的文獻結果與本程序預測結果均低于試驗值，但是本程序預測結果更接近于試驗值，但最大誤差仍然超過了25%，熱流密度預測的精度對現(xiàn)有程序來說仍然是一個嚴峻的挑戰(zhàn).

圖4 壁面壓力分布Fig.4 Distribution of the wall pressure

圖5 壁面熱流分布Fig.5 Distribution of the wall heat flux

2.2 高超聲速縫隙流動與換熱

以Allan 的縫隙流動與換熱試驗[17]為計算方案，驗證程序對高超聲速縫隙流動與換熱的預測能力.縫隙幾何參數如圖6 所示，縫隙的寬深比為0.383.來流參數為：Mach 數6.94，Reynolds 數7.45 × 106，壓強3 527 Pa，溫度154.6 K，攻角0°；出口設置為超音速出口，壁溫設置為300 K.

圖6 縫隙幾何參數（單位：mm）Fig.6 The gap geometry （unit:mm）

研究表明[18]，計算網格對熱流密度預測精度影響極大，第一層網格單元的網格Reynolds 數Rec保持在8 左右時可以保證計算的收斂性以及熱流的準確性，Rec=Re∞× Δn，其中Re∞與Δn分別表示來流Mach 數與壁面第一層網格至壁面距離；在計算中，采用如圖7 所示的拓撲結構進行網格劃分，壁面網格第一層單元距離壁面距離保持1.0 × 10-6m，網格增長率設置為1.1，最終網格單元數目為142 256 個.

圖7 縫隙流動網格拓撲結構示意圖Fig.7 The mesh topology for the hypersonic gap flow

如圖8 所示，縫隙入口段上方區(qū)域存在一個漩渦，在該漩渦的誘導作用下，縫隙下方區(qū)域又形成兩個強度較弱的漩渦，該分布符合試驗測量的結果.圖9 為縫隙中心線上流向速度沿縫深變化的曲線，橫坐標x表示距離縫隙入口的距離，即縫隙局部深度，該距離采用縫深參數d做無量綱化，縱坐標表示流向速度.從圖中可以看出：與縫隙內漩渦分布相對應，流動速度存在波動，在頂部區(qū)域速度波動較大，從170 m/s 迅速降至-60 m/s，但是在x/d約0.5 左右，流動速度的波動已經很小，不超過10 m/s 量級，隨著深度的增加，速度波動趨向于0，對流換熱對該區(qū)域的影響已經很小，傳熱以導熱方式為主.

圖8 縫隙內流線Fig.8 Streamlines in the gap

圖10 對比了計算得到的縫隙后壁面熱流分布，橫坐標定義與圖9 相同，縱坐標表示壁面熱流，采用該位置無縫隙時平板的熱流qref做無量綱化.從圖中可以看出，由于縫隙頂部存在強剪切流動，對流換熱效應顯著，局部熱流密度較大，隨著縫深的增加，對流換熱效用逐漸減弱，局部熱流密度迅速降低；預測的壁面熱流與試驗測量值的最大誤差約為14%.

圖9 縫隙中線流向速度沿縫深變化Fig.9 The velocity distribution along the gap central line

圖10 縫隙后壁面熱流分布Fig.10 The heat flux on the rear surface of the gap

2.3 無限長圓管高超聲速繞流非定常流動與傳熱

采用與Dechaumphai 等[19]、李佳偉等[20]的研究中相同的計算算例，來流參數為[21]：Mach 數6.47，壓強648.1 Pa，溫度241.5 K，Reynolds 數1.22 × 105，攻角0°，圓管導熱系數與密度參數按照不銹鋼材料進行選取.圓管的內外半徑分別為25.4 mm 和38.1 mm.采用結構化網格離散計算域，流體區(qū)域壁面第一層網格離開壁面距離設置為1.0 × 10-5m，網格增長率設置為1.1；固體區(qū)域沿著徑向均勻分布；最終流體與固體部分網格單元數目分別為64 × 168，32 × 168；計算網格如圖11 所示.

圖11 無限長圓管高超聲速繞流計算網格Fig.11 Computational grids for the unsteady coupled heat transfer simulation of hypersonic flow over infinite-length pipe

首先對流體域進行定常計算，計算中設置固體溫度恒定為294.4 K，收斂后流場參數以及恒定固體溫度294.4 K 作為非定常氣熱耦合計算的初始條件.模擬飛行時間為2 s，由于時間較短，圓管內壁面與氣體對流引起的熱量交換還比較弱，故在計算中將內壁面設置為絕熱壁面.文獻[20]的非定常耦合計算對本算例統(tǒng)一采用1.0 ×10-3s 的物理時間步長隱式推進，由于本文在固體域采用顯式推進，考慮到推進的穩(wěn)定性，流場與固體溫度場統(tǒng)一采用1.0 × 10-5s 的物理時間步長推進求解.

圖12 為非定常氣熱耦合計算得到的對稱線上不同時刻的溫度分布曲線.程序預測的激波厚度約為2 mm，在激波內，流體溫度從來流值241.5 K 迅速增加到2 100 K 量級，并對固體進行加熱.從曲線圖來看，溫度邊界層厚度約為1 mm；在邊界層以外區(qū)域，不同時刻預測的溫度分布差異很小，在邊界層內，溫度梯度很大，固體表面（對應駐點區(qū)域）溫度從初始溫度294.4 K 增加至390.8 K，而試驗測得的2 s 時刻駐點處溫度為388.72 K，誤差不超過0.53%，說明得到的固體溫度結果是可信的.

圖12 對稱線上溫度分布Fig.12 The temperature distribution along the centerline of the cylinder

圖13、圖14 給出了在t=0 s 時刻流體壁面熱流密度與壓力分布曲線（均采用駐點處參數做歸一化處理），可以看出所開發(fā)的計算程序得到的計算結果與試驗值吻合程度較好，所得到的流場結果與熱流結果是可信的；同時計算得到的t=0 s 時刻的駐點熱流密度為68.4 W/cm2，試驗值為69 W/cm2，誤差不超過1%.

圖13 t=0 s 時刻，流固交界面壓力分布Fig.13 The pressure on the fluid-solid interface at t=0 s

圖14 t=0 s 時刻，流固交界面熱流密度分布Fig.14 The heat flux density on the fluid-solid interface at t=0 s

3 結論

筆者基于有限差分法開發(fā)了高超聲速流動與換熱氣熱耦合求解器，運用該求解器分別對三個典型的高超聲速流動與換熱算例進行數值仿真，并將計算結果與試驗值進行了對比，得到如下結論：

1）以負攻角流過后臺階，在臺階拐點處形成膨脹波，臺階下游附近存在漩渦流動，氣體壓強、溫度等參數降低，熱流密度強度減弱，流動再附后形成再附激波，局部壓強、溫度提高，熱流密度驟升.

2）高超聲速氣流流過縫隙后，在縫隙入口處存在強烈的剪切運動，在頂部漩渦誘導作用下，縫隙底部也存在漩渦流動，但隨著縫深的增加，局部流動速度減小，對流換熱效應減弱.

3）高超聲速氣流在鈍體前方形成脫體激波，激波后氣體溫度迅速增加；流場內邊界層外氣動參數隨時間變化差異很小，溫度邊界層內存在較大的溫度梯度，壁面溫度隨時間持續(xù)升高.

4）三個算例仿真得到的氣動參數、壁面熱流密度的分布與試驗測量結果吻合得較好，本文開發(fā)的仿真求解器計算能力得到一定的驗證；由于激波與附面層的相互干擾，后臺階流動再附激波下游區(qū)域得到的熱流密度偏離實驗較大，熱流密度的預測精度有待提高.

致謝本文作者衷心感謝“特種環(huán)境復合材料技術國家級重點實驗室”基金（JCKYS2019603C003）對本文的資助.