高超聲速流動(dòng)與換熱數(shù)值仿真研究*

王 強(qiáng)，徐 濤，姚永濤

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

引言

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

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

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

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

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

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

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

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

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

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

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

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

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

2.1 高超聲速后臺(tái)階流動(dòng)與換熱

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

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

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

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

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

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

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

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

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

2.2 高超聲速縫隙流動(dòng)與換熱

以Allan 的縫隙流動(dòng)與換熱試驗(yàn)[17]為計(jì)算方案，驗(yàn)證程序?qū)Ω叱曀倏p隙流動(dòng)與換熱的預(yù)測(cè)能力.縫隙幾何參數(shù)如圖6 所示，縫隙的寬深比為0.383.來(lái)流參數(shù)為：Mach 數(shù)6.94，Reynolds 數(shù)7.45 × 106，壓強(qiáng)3 527 Pa，溫度154.6 K，攻角0°；出口設(shè)置為超音速出口，壁溫設(shè)置為300 K.

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

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

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

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

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

圖10 對(duì)比了計(jì)算得到的縫隙后壁面熱流分布，橫坐標(biāo)定義與圖9 相同，縱坐標(biāo)表示壁面熱流，采用該位置無(wú)縫隙時(shí)平板的熱流qref做無(wú)量綱化.從圖中可以看出，由于縫隙頂部存在強(qiáng)剪切流動(dòng)，對(duì)流換熱效應(yīng)顯著，局部熱流密度較大，隨著縫深的增加，對(duì)流換熱效用逐漸減弱，局部熱流密度迅速降低；預(yù)測(cè)的壁面熱流與試驗(yàn)測(cè)量值的最大誤差約為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 無(wú)限長(zhǎng)圓管高超聲速繞流非定常流動(dòng)與傳熱

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

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

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

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