王 強，徐 濤，姚永濤

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

引言

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

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

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

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

1 數(shù)學(xué)物理模型

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

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

固體導(dǎo)熱微分方程為

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

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

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

在定常計算中，當(dāng)參數(shù)的平均殘差小于10-4或殘差趨于水平后，即終止迭代；對非定常計算，當(dāng)?shù)锢頃r間達到總的物理時間后計算終止.

2 高超聲速流動與換熱數(shù)值仿真

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

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

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

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

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

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

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

氣動力與氣動熱是高超聲速飛行器設(shè)計的重要參數(shù)，圖4、圖5 分別給出了本程序預(yù)測的壁面壓力與壁面熱流參數(shù)，并與文獻中的計算結(jié)果、試驗結(jié)果進行了對比.在后臺階拐點之后，出現(xiàn)膨脹波系，加速流動導(dǎo)致局部的壓強、密度等參數(shù)驟降，壁面壓力值降低約70%，直到再附點附近，強壓縮流動導(dǎo)致局部流動參數(shù)急劇提升，與再附激波前壓強值相比，約增加2 倍.與壓強變化相比，壁面熱流密度變化具有相同的趨勢，但是變化的程度更大，從試驗結(jié)果來看，臺階拐點前后、再附激波前后等區(qū)域熱流密度值變化均在10 倍以上.從對比結(jié)果來看，在后臺階角點之前的入口段、分離流動再附之后的區(qū)域，本程序計算的壁面壓力與試驗結(jié)果、文獻計算結(jié)果吻合較好，誤差約9.6%，在10%以內(nèi)；但是臺階角點后、再附點之前的區(qū)域（50 mm ＜x＜ 80 mm）內(nèi)，文獻計算結(jié)果和本程序結(jié)果均與試驗值存在較大的誤差，該區(qū)域的預(yù)測也是目前湍流模型研究的一個難點問題，但是程序計算結(jié)果與文獻結(jié)果相比，誤差在10%以內(nèi).從熱流分布結(jié)果來看，本程序計算的熱流較文獻計算結(jié)果更接近于試驗結(jié)果，在50 mm ＜x＜ 80 mm 區(qū)域內(nèi)，預(yù)測的熱流與試驗值吻合很好，誤差在5%之內(nèi)，較大的誤差出現(xiàn)在再附激波以后的區(qū)域，該區(qū)域的文獻結(jié)果與本程序預(yù)測結(jié)果均低于試驗值，但是本程序預(yù)測結(jié)果更接近于試驗值，但最大誤差仍然超過了25%，熱流密度預(yù)測的精度對現(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]為計算方案，驗證程序?qū)Ω叱曀倏p隙流動與換熱的預(yù)測能力.縫隙幾何參數(shù)如圖6 所示，縫隙的寬深比為0.383.來流參數(shù)為：Mach 數(shù)6.94，Reynolds 數(shù)7.45 × 106，壓強3 527 Pa，溫度154.6 K，攻角0°；出口設(shè)置為超音速出口，壁溫設(shè)置為300 K.

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

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

圖7 縫隙流動網(wǎng)格拓撲結(jié)構(gòu)示意圖Fig.7 The mesh topology for the hypersonic gap flow

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

圖8 縫隙內(nèi)流線Fig.8 Streamlines in the gap

圖10 對比了計算得到的縫隙后壁面熱流分布，橫坐標(biāo)定義與圖9 相同，縱坐標(biāo)表示壁面熱流，采用該位置無縫隙時平板的熱流qref做無量綱化.從圖中可以看出，由于縫隙頂部存在強剪切流動，對流換熱效應(yīng)顯著，局部熱流密度較大，隨著縫深的增加，對流換熱效用逐漸減弱，局部熱流密度迅速降低；預(yù)測的壁面熱流與試驗測量值的最大誤差約為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]的研究中相同的計算算例，來流參數(shù)為[21]：Mach 數(shù)6.47，壓強648.1 Pa，溫度241.5 K，Reynolds 數(shù)1.22 × 105，攻角0°，圓管導(dǎo)熱系數(shù)與密度參數(shù)按照不銹鋼材料進行選取.圓管的內(nèi)外半徑分別為25.4 mm 和38.1 mm.采用結(jié)構(gòu)化網(wǎng)格離散計算域，流體區(qū)域壁面第一層網(wǎng)格離開壁面距離設(shè)置為1.0 × 10-5m，網(wǎng)格增長率設(shè)置為1.1；固體區(qū)域沿著徑向均勻分布；最終流體與固體部分網(wǎng)格單元數(shù)目分別為64 × 168，32 × 168；計算網(wǎng)格如圖11 所示.

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

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

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

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

圖13、圖14 給出了在t=0 s 時刻流體壁面熱流密度與壓力分布曲線（均采用駐點處參數(shù)做歸一化處理），可以看出所開發(fā)的計算程序得到的計算結(jié)果與試驗值吻合程度較好，所得到的流場結(jié)果與熱流結(jié)果是可信的；同時計算得到的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 結(jié)論

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

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

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

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

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

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