高超聲速?gòu)?fù)雜氣動(dòng)問題數(shù)值方法研究進(jìn)展

王江峰， 伍貽兆， 季衛(wèi)棟， 樊孝峰， 趙法明， 呂偵軍

南京航空航天大學(xué) 航空宇航學(xué)院， 南京 210016

高超聲速?gòu)?fù)雜氣動(dòng)問題數(shù)值方法研究進(jìn)展

王江峰*， 伍貽兆， 季衛(wèi)棟， 樊孝峰， 趙法明， 呂偵軍

南京航空航天大學(xué) 航空宇航學(xué)院， 南京 210016

高超聲速流場(chǎng)具有復(fù)雜流動(dòng)特征，其中真實(shí)氣體效應(yīng)、磁流體干擾效應(yīng)和力熱結(jié)構(gòu)耦合效應(yīng)等對(duì)氣動(dòng)力分析產(chǎn)生了重要影響。將流體力學(xué)研究擴(kuò)展到分子動(dòng)力學(xué)、電磁流體力學(xué)以及流固耦合等交叉學(xué)科領(lǐng)域，這給數(shù)值模擬方法帶來了巨大挑戰(zhàn)。針對(duì)高超聲速氣動(dòng)力/熱分析的熱點(diǎn)問題，重點(diǎn)關(guān)注高溫效應(yīng)與低密度流動(dòng)效應(yīng)、磁流體干擾效應(yīng)和力熱結(jié)構(gòu)耦合效應(yīng)等，結(jié)合算例分析了相應(yīng)的數(shù)值求解技術(shù)；在氣動(dòng)熱方面主要比較了3類求解方法(純工程方法、純數(shù)值方法和基于Prandtl邊界層理論的方法)，并給出了相應(yīng)算例；對(duì)于氣動(dòng)力/熱/結(jié)構(gòu)耦合問題，從耦合模型及耦合計(jì)算方法兩方面開展了分析。最后指出了高超聲速?gòu)?fù)雜氣動(dòng)問題數(shù)值求解技術(shù)未來需重點(diǎn)關(guān)注的幾個(gè)方面。

高超聲速； 數(shù)值模擬； 氣動(dòng)加熱； 磁流體力學(xué)； 氣動(dòng)熱彈性力學(xué)； 氣動(dòng)力/熱/結(jié)構(gòu)耦合

近年來，高超聲速飛行器的戰(zhàn)略價(jià)值受到了各科技大國(guó)的高度重視，成為各國(guó)近期的研究熱點(diǎn)。相比于超聲速流動(dòng)，高超聲速流動(dòng)具有以下幾個(gè)顯著特征：①薄激波層；②熵層；③黏性干擾效應(yīng)；④高溫效應(yīng)；⑤低密度流動(dòng)效應(yīng)[1-2]。這些特點(diǎn)使得應(yīng)用于一般超聲速流動(dòng)的數(shù)值模擬技術(shù)已不再適用。高超聲速飛行中出現(xiàn)的各種特殊流動(dòng)現(xiàn)象，如激波與邊界層的相互干擾、邊界層傳熱傳質(zhì)、化學(xué)反應(yīng)以及燒蝕等，受到了國(guó)內(nèi)外學(xué)者的高度重視，藉此，高超聲速氣動(dòng)問題的數(shù)值求解技術(shù)得到了快速的發(fā)展。求解的控制方程已經(jīng)從最初的Euler方程、各種簡(jiǎn)化Navier-Stokes方程到全Navier-Stokes方程；湍流尺度從雷諾平均方法、大渦模擬到直接數(shù)值模擬；氣體熱力學(xué)狀態(tài)模型也已經(jīng)從早期的量熱完全氣體模型發(fā)展到熱力學(xué)化學(xué)非平衡模型；化學(xué)反應(yīng)模型從最基本的單步基元反應(yīng)到復(fù)雜的多步復(fù)合反應(yīng)，包含數(shù)百個(gè)組元以及數(shù)千個(gè)反應(yīng)式；計(jì)算網(wǎng)格也從單塊結(jié)構(gòu)網(wǎng)格發(fā)展到多塊、非結(jié)構(gòu)網(wǎng)格和重疊網(wǎng)格等。數(shù)值格式的研究和應(yīng)用更是蓬勃發(fā)展，從傳統(tǒng)的中心格式發(fā)展到現(xiàn)在的多種高效迎風(fēng)格式。

高超聲速的復(fù)雜氣動(dòng)問題涉及眾多學(xué)科專業(yè)，如高超聲速空氣動(dòng)力學(xué)、計(jì)算流體力學(xué)(CFD)、化學(xué)反應(yīng)流體動(dòng)力學(xué)、電磁流體力學(xué)以及氣動(dòng)光學(xué)等，其數(shù)值求解方法也隨著所關(guān)注問題的重點(diǎn)而異。本文就高超聲速流動(dòng)中所涉及的氣動(dòng)力模擬、化學(xué)反應(yīng)效應(yīng)、磁流體力學(xué)、氣動(dòng)加熱以及流固耦合問題開展分析，旨在與高超聲速?gòu)?fù)雜流動(dòng)數(shù)值模擬技術(shù)領(lǐng)域的研究人員共同探討。

1 高超聲速?gòu)?fù)雜流動(dòng)數(shù)值模擬方法

1.1 控制方程

高超聲速連續(xù)流的主要流動(dòng)控制方程是從質(zhì)量守恒、動(dòng)量守恒和能量守恒三大定律導(dǎo)出的。

下面給出流體力學(xué)中最常用的從空間位置固定的有限控制體模型導(dǎo)出的流體控制方程——積分形式的Navier-Stokes方程：

(1)

式中：Ω為控制體； ?Ω為控制體單元的邊界； dS為控制體單元面；W為守恒量；Fc為對(duì)流通量;Fv為黏性通量;Q為源項(xiàng)。三維守恒量、對(duì)流通量、黏性通量及源項(xiàng)的表達(dá)式為

(2)

(3)

對(duì)于真實(shí)氣體，還需要考慮空氣中不同組分之間的化學(xué)反應(yīng)效應(yīng)。考慮N種不同組分的混合氣體，則需要在Navier-Stokes方程中額外添加N-1個(gè)輸運(yùn)方程。向量W、Fc、Fv和Q的新表達(dá)式分別為

(4)

其中：

(5)

當(dāng)飛行器飛行高度增高，空氣密度下降到一定程度后，連續(xù)流的假設(shè)不再滿足真實(shí)的流動(dòng)環(huán)境，這時(shí)候需要考慮稀薄氣體的影響。高超聲速稀薄流的流動(dòng)控制方程為Boltzmann方程，由于Boltzmann方程本身的高度非線性，傳統(tǒng)數(shù)值求解方法在實(shí)際運(yùn)用時(shí)遇到了極大的困難。而基于統(tǒng)計(jì)思想的直接模擬蒙特卡羅(DirectSimulationMonteCarlo，DSMC)方法，從微觀角度建立模擬分子的碰撞、遷移和能量交換等過程，是目前最有能力模擬現(xiàn)實(shí)情況下三維高超聲速稀薄氣體流的方法。下面給出Boltzmann方程：

(6)

磁流體力學(xué)(Magnetohydrodynamics，MHD)作為電磁學(xué)與流體力學(xué)的交叉學(xué)科，其基本控制方程同樣由三大守恒定律給出。MHD的控制方程在電磁介質(zhì)各向同性、電中性非相對(duì)流動(dòng)等假設(shè)下可推導(dǎo)出全MHD方程組，形式如下：

(7)式中：B為磁感應(yīng)強(qiáng)度;I為單位矩陣；τ為黏性應(yīng)力張量；μe為磁導(dǎo)率；ve為磁擴(kuò)散率。

氣動(dòng)熱彈性力學(xué)作為流體力學(xué)、傳熱學(xué)和結(jié)構(gòu)(動(dòng))力學(xué)的耦合交叉學(xué)科，除流體力學(xué)的控制方程組外，還引入了傳熱學(xué)和結(jié)構(gòu)(動(dòng))力學(xué)方程：

(8)

(9)

1.2 離散方法

對(duì)于Navier-Stokes方程的時(shí)間導(dǎo)數(shù)項(xiàng)，有限體積解法中常采用的離散格式分為顯式格式和隱式格式。顯式格式主要有多步龍格-庫(kù)塔(Runge-Kutta)迭代法[3]和混合多步迭代法[4]。隱式格式主要有ADI(Alternating Direction Implicit)方法[5]、LU-SGS(Lower-Upper Symmetric Gauss-Seidel)方法[6-7]和克雷洛夫子空間迭代法(Krylov-subspace Method)[8]。

空間離散分為中心格式(Central Scheme)[3]、迎風(fēng)格式(Upwind Scheme)[9-13]、TVD(Total Variation Diminishing)格式[14]。迎風(fēng)格式又可分為矢通量分裂格式(Flux-vector Splitting Scheme)、通量差分分裂格式(Flux-difference Splitting Scheme)。

1.3 湍流模型

根據(jù)模擬的湍流尺度的不同，Navier-Stokes方程的求解可以分為以下3種：雷諾平均Navier-Stokes(Reynolds-averaged Navier-Stokes, RANS)方程、大渦模擬(Large Eddy Simulation, LES)[15]及直接數(shù)值模擬(Direct Numerical Simulation, DNS)[16]。RANS方法是目前應(yīng)用最廣的方法，通過合理選用湍流模型，能夠很好地模擬大多數(shù)流動(dòng)現(xiàn)象，且計(jì)算效率高。由于RANS方程的不封閉性，人們引入了湍流模型來封閉方程組。湍流模型一般分為以下幾種：零方程模型[17-18]、一方程模型[19-20]、二方程模型[21-24]及雷諾應(yīng)力模型。

2 氣動(dòng)力計(jì)算

氣動(dòng)力計(jì)算是高超聲速氣動(dòng)問題數(shù)值求解的一個(gè)核心問題，本節(jié)結(jié)合算例，給出高超聲速飛行器氣動(dòng)力分析及考慮高溫效應(yīng)(化學(xué)非平衡流、化學(xué)動(dòng)力學(xué)模型)與低密度流動(dòng)效應(yīng)，磁流體干擾效應(yīng)等影響的求解技術(shù)。

2.1 高超聲速飛行器氣動(dòng)力流場(chǎng)分析

以吸氣式?jīng)_壓發(fā)動(dòng)機(jī)為動(dòng)力的高超聲速飛行器流場(chǎng)特性計(jì)算為例，其流場(chǎng)包含前體壓縮激波、激波-附面層干擾、隔離段激波串以及后體-尾噴羽流等復(fù)雜的流動(dòng)現(xiàn)象，這對(duì)數(shù)值計(jì)算方法的精度與效率提出了很高的要求。

X-43A飛行成功后，國(guó)內(nèi)針對(duì)其氣動(dòng)力特性開展了許多研究。圖1給出了類X-43A通氣外形在馬赫數(shù)Ma=6.976 5、迎角α=1.5°情況下的前體壓縮和后體尾噴流場(chǎng)結(jié)構(gòu)。采用基于混合網(wǎng)格技術(shù)的三維Navier-Stokes方程數(shù)值求解方法，空間離散采用VanLeer迎風(fēng)格式，時(shí)間推進(jìn)為多步顯示格式，Spalart-Allmaras湍流模型，采用了分布式并行計(jì)算技術(shù)用以提高計(jì)算效率。圖中展示了前體壓縮流場(chǎng)和后體尾噴流流場(chǎng)的幾個(gè)典型截面的等馬赫數(shù)線分布。

圖1 類X-43A飛行器流場(chǎng)結(jié)構(gòu)(馬赫數(shù)Ma=6.976 5、迎角α=1.5°)Fig.1 Flow field structure around X-43A type aircraft (Mach mumber Ma=6.976 5, angle of attack α=1.5°)

2.2 化學(xué)非平衡流

在真實(shí)氣體效應(yīng)的CFD研究中，空氣化學(xué)非平衡流動(dòng)不僅因其機(jī)理與現(xiàn)象的復(fù)雜性[1-2]，更因?yàn)樗菬g、壁面催化和推進(jìn)劑燃料燃燒等現(xiàn)象的基礎(chǔ)而受到高度重視，國(guó)內(nèi)外廣泛開展了這一方面的研究。Blottner[25]針對(duì)尖錐外形的再入飛行器作了關(guān)于電離空氣非平衡層流邊界層流動(dòng)的研究，指出在靠近尖錐區(qū)域的流場(chǎng)有強(qiáng)烈的化學(xué)非平衡效應(yīng)。國(guó)內(nèi)外多位學(xué)者將成熟的計(jì)算格式擴(kuò)展到高超聲速化學(xué)非平衡流中，取得了良好的效果。張向洪等[26]將HLLE+格式(Harten-Lax-van Leer-Einfeldt plus scheme)擴(kuò)展應(yīng)用在高超聲速多組分化學(xué)非平衡流動(dòng)的數(shù)值模擬，保留了HLLE+格式數(shù)值精度的同時(shí)，并且具有較高的穩(wěn)定性，計(jì)算結(jié)果與實(shí)驗(yàn)值符合良好。劉晨等[27]將AUFS格式(Artificially Upstream Fluy vector Splitting scheme)應(yīng)用于鈍頭體激波誘導(dǎo)超聲速爆轟燃燒流場(chǎng)進(jìn)行了數(shù)值模擬，反應(yīng)模型采用經(jīng)典的氫氧7組元8反應(yīng)模型，得到了與實(shí)驗(yàn)相符合的計(jì)算結(jié)果(如圖2所示，球頭圓柱體直徑為15 mm，馬赫數(shù)Ma=6.976 5，溫度T=250 K，壓力p=426 62 Pa，來流為氫氣與空氣的混合氣體：H2∶O2∶N2=2∶1∶3.76，圖中x為沿駐點(diǎn)線坐標(biāo)，R為球頭圓柱體半徑)，說明AUFS格式能較好捕捉超聲速爆轟流場(chǎng)中的激波和燃燒波等間斷。

圖2 球頭圓柱體算例[27]Fig.2 Case of cylinder model[27]

柳軍等[28]將迎風(fēng)型NND格式(Non-oscillatory and Non-free-parameter Dissipation difference scheme)應(yīng)用于三維高超聲速熱化學(xué)非平衡數(shù)值模擬中，采用合適的熵修正方法消除了高超聲速流數(shù)值模擬中的“Carbuncle”現(xiàn)象，對(duì)高超聲速圓球飛行流場(chǎng)進(jìn)行的數(shù)值模擬結(jié)果表明了計(jì)算方法的有效性。

本文采用基于混合網(wǎng)格的三維多組Navier-Stotes方程求解方法，對(duì)簡(jiǎn)化航天飛機(jī)頭部的雙橢球外形進(jìn)行了數(shù)值計(jì)算，如圖3所示。計(jì)算條件為馬赫數(shù)Ma=20、飛行高度h=60 km，考慮了空氣11組元化學(xué)反應(yīng)流效應(yīng)。計(jì)算網(wǎng)格為混合網(wǎng)格(內(nèi)層為六面體結(jié)構(gòu)網(wǎng)格、外層為四面體非結(jié)構(gòu)網(wǎng)格)?？臻g離散采用VanLeer迎風(fēng)格式，時(shí)間推進(jìn)為考慮了流動(dòng)時(shí)間與化學(xué)反應(yīng)時(shí)間的多步顯示格式，采用分布式并行計(jì)算技術(shù)用以提高計(jì)算效率。圖3中給出了典型組元(O2,N2,NO,NO+,e-)的密度分布。

圖3 雙橢球化學(xué)非平衡繞流(Ma=20,α=1.5°,飛行高度h=60 km)Fig.3 Chemical non-equilibrium flow around double ellipsoid (Ma=20,α=1.5°,flight altitude h=60 km)

2.3 化學(xué)動(dòng)力學(xué)模型

在經(jīng)典熱力學(xué)和可壓縮流動(dòng)研究中，通常假定氣體的比熱為常數(shù)，由此，比熱比也是常數(shù)，在這些假定下的氣體為量熱完全氣體。然而，在高超聲速飛行條件下，飛行器的動(dòng)能使得其周圍的空氣達(dá)到很高的溫度，高溫使得氣體分子的振動(dòng)能被激發(fā)，甚至引起氣體分子的離解、電離等復(fù)雜的化學(xué)反應(yīng)，這些反應(yīng)會(huì)影響空氣的特性并使之偏離量熱完全氣體或者熱完全氣體。此時(shí)的比熱不再是常數(shù)，它不僅取決于氣體的溫度和壓強(qiáng)，在非平衡情況下還依賴于組元濃度，這種情況下的氣體為化學(xué)反應(yīng)完全氣體混合物。為了準(zhǔn)確描述高超聲速流動(dòng)下的流場(chǎng)結(jié)構(gòu)，觀察流動(dòng)現(xiàn)象，化學(xué)動(dòng)力學(xué)模型的選擇就顯得尤為重要。Herzberg[29]和Moore[30]從光譜數(shù)據(jù)分析，確定了空氣分子以及原子的能級(jí)，是今后研究平衡熱力學(xué)特性參數(shù)和輸運(yùn)參數(shù)的基礎(chǔ)。Hansen和Heims[31]總結(jié)前人的研究成果，針對(duì)高溫空氣化學(xué)動(dòng)力學(xué)模型開展了一系列研究并與飛行試驗(yàn)結(jié)果作了對(duì)比。結(jié)果表明采用迭代算法時(shí)，平衡熱力學(xué)特性參數(shù)的誤差在0.5%之內(nèi)，采用半經(jīng)驗(yàn)公式計(jì)算時(shí)，其誤差在10%～20%之間；而針對(duì)輸運(yùn)特性以及化學(xué)反應(yīng)速率的研究并沒有得到令人滿意的結(jié)果。Gupta等[32]考慮11組元空氣熱化學(xué)非平衡模型，總結(jié)了空氣部分電離時(shí)，熱力學(xué)特性和輸運(yùn)特性的近似和精確計(jì)算公式，同時(shí)指出這些公式使用時(shí)的限制，還指出在化學(xué)反應(yīng)速率的計(jì)算中，仍然存在較大的誤差。不同的化學(xué)動(dòng)力學(xué)模型對(duì)化學(xué)反應(yīng)流場(chǎng)的數(shù)值模擬起最直接的影響。

近二十年來，隨著實(shí)驗(yàn)測(cè)量技術(shù)以及計(jì)算機(jī)技術(shù)的發(fā)展，國(guó)內(nèi)外在燃燒化學(xué)動(dòng)力學(xué)模型的研究中也有不小的進(jìn)展。Eklund等[33]針對(duì)膨脹型噴射斜坡，采用Jachimowski提出的7組元7反應(yīng)有限速率氫氧化學(xué)反應(yīng)動(dòng)力學(xué)模型，數(shù)值求解了化學(xué)反應(yīng)流場(chǎng)結(jié)構(gòu)并與實(shí)驗(yàn)結(jié)果相互印證，得到了相符合的結(jié)果。Clutter等[34]針對(duì)飛行馬赫數(shù)Ma=6.46的鈍頭彈體，考慮氫氣空氣混合，研究了不同化學(xué)動(dòng)力學(xué)模型對(duì)激波誘導(dǎo)燃燒流場(chǎng)激波結(jié)構(gòu)的影響，研究指出壓強(qiáng)是影響數(shù)值模擬結(jié)果的關(guān)鍵因素。Ingram等[35]采用簡(jiǎn)化過的燃燒反應(yīng)模型模擬了甲烷在喉道的燃燒，研究指出合理簡(jiǎn)化的燃燒反應(yīng)模型可以反映包含重要中間反應(yīng)在內(nèi)的整個(gè)燃燒過程。劉晨等[36]在對(duì)鈍體激波誘導(dǎo)燃燒模擬中發(fā)現(xiàn)不同反應(yīng)模型在靠近對(duì)稱軸位置會(huì)產(chǎn)生類似紅玉現(xiàn)象的數(shù)值異常，指出采用經(jīng)驗(yàn)證較為可靠的模型可有效改善這類數(shù)值異常。

2.4 稀薄流效應(yīng)

DSMC方法由澳大利亞學(xué)者Bird于1963年在分子動(dòng)力學(xué)方法基礎(chǔ)上提出[37]，該方法直接以分子為研究對(duì)象，最早用來模擬均勻氣體中的松弛問題，后來發(fā)展到模擬二維、三維復(fù)雜的真實(shí)氣體流動(dòng)，并且成功應(yīng)用于美國(guó)航天飛機(jī)升阻比計(jì)算，所得結(jié)果與飛行測(cè)量數(shù)據(jù)相符[38]。

國(guó)內(nèi)最早由中國(guó)科學(xué)院力學(xué)研究所的沈青[39]、樊菁等[40-41]針對(duì)DSMC方法開展相關(guān)研究，在基礎(chǔ)理論和工程實(shí)踐應(yīng)用方面均取得較為顯著的成果。此后，伴隨我國(guó)航天事業(yè)的飛速發(fā)展，尤其是近來臨近空間飛行器概念的提出，過渡流域飛行器氣動(dòng)力、氣動(dòng)熱預(yù)估等相關(guān)問題的解決迫在眉睫，DSMC方法愈受國(guó)內(nèi)學(xué)者的重視，相關(guān)研究迅速開展，研究重點(diǎn)涉及新的分子碰撞模型的建立[42]、不同網(wǎng)格技術(shù)的發(fā)展、DSMC與連續(xù)流流場(chǎng)求解技術(shù)的耦合[43]等方面。王學(xué)德等[44-45]基于非結(jié)構(gòu)網(wǎng)格技術(shù)，提出了面積元/體積元坐標(biāo)搜索算法與交替二叉樹搜索算法結(jié)合的分子搜索算法，設(shè)計(jì)GSS-3(Generalized Soft Sphere)分子碰撞模型，提出動(dòng)態(tài)局部時(shí)間步長(zhǎng)技術(shù)提高程序效率，加入化學(xué)反應(yīng)模型，并引入并行計(jì)算技術(shù)，最終成功模擬了多個(gè)類型高超聲速飛行器算例。圖4給出了X-37B外形在稀薄流區(qū)多組元?dú)怏w模型下的氣動(dòng)特性。計(jì)算網(wǎng)格采用非結(jié)構(gòu)網(wǎng)格，采用基于PCs-Cluster的分布式并行計(jì)算技術(shù)用以提高計(jì)算效率。圖中給出了飛行器表面壓力分布特性。

圖4 X-37B表面壓力分布(α=25°， h=90 km，來流速度V∞=7 km/s， N2∶O2∶O=0.788∶0.209∶0.003)Fig.4 Pressure distribution on surface of X-37B(α=25°， h=90 km， free stream velocity V∞=7 km/s, N2∶O2∶O=0.788∶0.209∶0.003)

2.5 MHD數(shù)值模擬方法

MHD主要研究導(dǎo)電流體與磁場(chǎng)的相互作用，是流體力學(xué)與電磁學(xué)的一門交叉學(xué)科，其在高超聲速流動(dòng)控制領(lǐng)域的應(yīng)用是當(dāng)前國(guó)內(nèi)外研究的熱點(diǎn)，眾多國(guó)內(nèi)外學(xué)者對(duì)MHD控制系統(tǒng)在高超聲速飛行器減阻與熱防護(hù)、超燃沖壓發(fā)動(dòng)機(jī)進(jìn)氣道流場(chǎng)控制、磁流體能量分流以及發(fā)動(dòng)機(jī)燃燒室燃料促混等方面的應(yīng)用作了深入研究，取得了大量成果。

MHD數(shù)值模擬包含了低磁雷諾數(shù)MHD方程數(shù)值模擬與全MHD方程的數(shù)值模擬，前者的計(jì)算方法與Euler方程類似，但是后者卻要復(fù)雜的多。這主要是由于[46-47]：①多維全MHD方程是奇異的，并且其特征雅可比矩陣具有非齊次性，因此無法應(yīng)用傳統(tǒng)的迎風(fēng)格式進(jìn)行通量分裂；②全MHD方程在求解過程中會(huì)由于數(shù)值誤差產(chǎn)生額外的磁場(chǎng)散度，最終可能導(dǎo)致計(jì)算發(fā)散。針對(duì)上述數(shù)值困難，國(guó)內(nèi)外的研究者們作了很多嘗試，主要可以歸結(jié)為如下的幾類處理方法：

1)Brackbill和Barnes[48]提出了投影方法，該方法通過求解磁場(chǎng)勢(shì)函數(shù)的泊松方程，達(dá)到控制磁場(chǎng)散度的目的。

2)Evans和Hawley[49]以交錯(cuò)網(wǎng)格方法為基礎(chǔ)，提出了約束傳遞(Constrained Transport)磁場(chǎng)清理技術(shù)。

3)Powell等[50]修正了MHD方程雅可比矩陣奇異性問題，提出了八波形式的MHD方程組模型，計(jì)算結(jié)果表明該方法能提高計(jì)算的穩(wěn)定性，但是對(duì)磁場(chǎng)散度的抑制效果并不太好，計(jì)算域中的磁場(chǎng)散度的峰值仍然較大。

4)Dedner等[51]通過引入拉格朗日算子與約束變量，構(gòu)造了雙曲型EGLM(Extended Generalized Lagrange Multiplier)-MHD方程，在EGLM-MHD方程求解過程中，磁場(chǎng)散度可以得到有效的控制。

近幾年國(guó)內(nèi)也有部分研究人員做了相關(guān)的工作，潘勇[52]以非結(jié)構(gòu)網(wǎng)格有限體積算法為基礎(chǔ)，應(yīng)用雙曲型散度清除技術(shù)計(jì)算了高超聲速鈍頭體MHD繞流流場(chǎng)，得到了非常理想的結(jié)果，如圖5所示；田正雨等[53]對(duì)投影法作了深入研究，發(fā)現(xiàn)投影法會(huì)對(duì)平滑流場(chǎng)區(qū)域帶來相對(duì)較大的誤差，并以此為據(jù)提出了局部投影法。張向洪[47]針對(duì)三維MHD流動(dòng)模擬技術(shù)開展了研究，分析了Rotor問題及磁場(chǎng)對(duì)壓縮管道的影響特性，如圖6和圖7所示。圖6中左圖為密度云圖，右圖為馬赫數(shù)云圖；圖7中的參考數(shù)據(jù)參見文獻(xiàn)[50]。

2.6 小 結(jié)

近幾十年來，高溫化學(xué)反應(yīng)流場(chǎng)的數(shù)值模擬技術(shù)已經(jīng)取得了可觀的進(jìn)展，并已經(jīng)在工程實(shí)際中得到了應(yīng)用，取得了令人滿意的成果，但仍有許多問題亟待解決。

1)對(duì)于化學(xué)非平衡流動(dòng)，除了傳統(tǒng)的求解大型非線性方程組所遇到的困難外，組元連續(xù)方程中的化學(xué)源項(xiàng)進(jìn)一步增加了計(jì)算的難度。當(dāng)組元的化學(xué)反應(yīng)特征時(shí)間和流場(chǎng)特征時(shí)間相差幾個(gè)量級(jí)時(shí)，不可避免的會(huì)出現(xiàn)剛性問題。另外，數(shù)值格式的穩(wěn)定性很大程度上依賴于格式的耗散性，然而數(shù)值耗散是非物理的，過多的數(shù)值耗散很可能會(huì)掩蓋真實(shí)的物理現(xiàn)象，而當(dāng)數(shù)值精度提高時(shí)，會(huì)造成穩(wěn)定性降低，這二者之間往往是相互矛盾的。如何協(xié)調(diào)數(shù)值精度與數(shù)值穩(wěn)定性是另一研究重點(diǎn)。

圖5 磁流體力學(xué)(MHD)流場(chǎng)散度分布[52]Fig.5 Distribution of divergence in magnetohydrodynamics (MHD) flowfield[52]

圖7 三維高超聲速M(fèi)HD壓縮管道[47]

2)化學(xué)動(dòng)力學(xué)模型對(duì)化學(xué)反應(yīng)流場(chǎng)的數(shù)值模擬結(jié)果往往起到?jīng)Q定性的作用，而模型中各類參數(shù)的確定往往是通過統(tǒng)計(jì)學(xué)理論或者半經(jīng)驗(yàn)公式，缺乏統(tǒng)一的標(biāo)準(zhǔn)，因此統(tǒng)一建模是目前的一大難題。

3)高超聲速M(fèi)HD流場(chǎng)屬于復(fù)雜的多物理耦合場(chǎng)，數(shù)值計(jì)算的復(fù)雜性很大，還有很多工作需要深入開展，如強(qiáng)電場(chǎng)對(duì)高超聲速流場(chǎng)的影響，全MHD方程求解的效率、穩(wěn)定性以及MHD流場(chǎng)控制的應(yīng)用等。

3 氣動(dòng)熱方面

飛行馬赫數(shù)越大，高溫氣流對(duì)飛行器表面的加熱程度就越嚴(yán)重，改變飛行器結(jié)構(gòu)的強(qiáng)度和剛度，可能導(dǎo)致飛行器氣動(dòng)外形發(fā)生變化，對(duì)飛行器的正常飛行產(chǎn)生嚴(yán)重影響。氣動(dòng)加熱現(xiàn)象極為復(fù)雜，既與飛行器的飛行速度、高度和附面層狀態(tài)有關(guān)，又涉及到激波和附面層干擾、滑流、駐點(diǎn)熱交換等問題。目前，氣動(dòng)熱的計(jì)算方法主要有3類：①純工程方法[54-57]；②純數(shù)值方法[58-62]，直接求解Navier-Stokes方程及其近似形式；③基于Prandtl邊界層理論的邊界層外無黏外流解與邊界層內(nèi)工程方法相結(jié)合的方法[63]。這3種方法在效率和精度方面各有優(yōu)缺點(diǎn)。

3.1 純工程方法

純工程方法對(duì)簡(jiǎn)單外形的氣動(dòng)熱求解具有很高的效率及可靠的精度，但對(duì)復(fù)雜外形的氣動(dòng)熱問題適應(yīng)性較差。國(guó)內(nèi)外學(xué)者在拓展純工程方法的使用范圍上付出了極大努力，現(xiàn)在純工程方法已基本能實(shí)現(xiàn)帶迎角的全機(jī)外形的氣動(dòng)熱預(yù)測(cè)。

de Jarnette 和Hamilton[54]發(fā)展了一套適用于理想氣體和空氣平衡氣體的計(jì)算高超聲速流中任意三維鈍頭體層流、轉(zhuǎn)捩和湍流加熱率的方法。該方法采用軸對(duì)稱比擬法，在已知表面流線的情況下，通過求解軸對(duì)稱邊界方程的方法得到層流和湍流加熱率。美國(guó)國(guó)家航空航天局(NASA)Langley 研究中心開發(fā)的平衡多功能性、計(jì)算效率以及精度的一款氣動(dòng)熱預(yù)測(cè)程序——LMINIVER，后來為了將LMINIVER整合到AVID系統(tǒng)中，將LMINIVER升級(jí)為ANMIN程序[55]。LMINIVER程序在駐點(diǎn)區(qū)域采用了經(jīng)典的Fay-Riddle公式，層流采用Blasius表面摩阻公式，湍流用Schultz-Grunow 表面摩阻公式，并用參考焓方法來考慮高速層流壓縮性的影響。Zoby和Simmonds[56]發(fā)展了一種用于傳熱預(yù)測(cè)和衡量熱傳遞中的可變熵邊界條件影響的近似無黏流方法——LATCH。該程序可以計(jì)算零迎角的雙曲面、橢球、拋物面和球錐的無黏流場(chǎng)與對(duì)流加熱。有迎角情況下，通過考慮球錐迎風(fēng)面和背風(fēng)面的流線分布影響及圓周熱流率來近似。李建林等[57]采用工程算法對(duì)升力體和乘波體構(gòu)型飛行器開展了氣動(dòng)熱計(jì)算，得到的結(jié)果與數(shù)值計(jì)算的結(jié)果相近，表明工程算法能滿足工程估算的精度需求。

3.2 純數(shù)值方法

高超聲速氣動(dòng)熱的數(shù)值模擬向來是CFD研究中的難點(diǎn)，首先離散數(shù)值方法在熱-化學(xué)模型及流動(dòng)物理模型(層流-湍流轉(zhuǎn)捩、湍流和湍流分離)方面存在不足，其次氣動(dòng)熱的計(jì)算受數(shù)值格式、網(wǎng)格分布、收斂過程甚至是熱流的后處理計(jì)算等各方面因素的影響。這些因素導(dǎo)致氣動(dòng)熱的數(shù)值模擬十分困難且對(duì)計(jì)算資源需求巨大。

連續(xù)流區(qū)數(shù)值求解Navier-Stokes方程的方法通常包括直接求解黏性激波層(Viscous Shock Layer，VSL)方程[58]、拋物線化的Navier-Stokes(Parabolized Navier-Stokes，PNS)方程[59-60]和Navier-Stokes 方程[61]幾種。稀薄流區(qū)一般采用DSMC方法[62]。

對(duì)于連續(xù)流區(qū)氣動(dòng)熱計(jì)算的數(shù)值格式、網(wǎng)格分布、收斂過程方面的問題，國(guó)內(nèi)一些學(xué)者也開展了相關(guān)研究。潘沙等[64]開展了氣動(dòng)熱數(shù)值模擬中的網(wǎng)格相關(guān)性及收斂性研究，結(jié)果表明網(wǎng)格是氣動(dòng)熱數(shù)值模擬中的關(guān)鍵因素，壁面附近法向網(wǎng)格最為敏感，熱流結(jié)果隨網(wǎng)格間距不同會(huì)出現(xiàn)數(shù)倍乃至數(shù)量級(jí)的差異；以方程殘值的下降量級(jí)和壓力收斂作為熱流收斂的判別標(biāo)準(zhǔn)是不合適的，應(yīng)直接觀察熱流數(shù)據(jù)的收斂，確保得到真正的收斂解。閻超等[65-66]開展了熱流CFD計(jì)算中格式和網(wǎng)格效應(yīng)的研究，結(jié)果表明不同格式對(duì)熱流計(jì)算精度有影響，AUSM+格式在熱流計(jì)算精確性方面占優(yōu)(見圖8，圖中Q/Qsph為熱流與駐點(diǎn)熱流的比值，L為鈍雙錐的長(zhǎng)度)；與格式效應(yīng)相比，網(wǎng)格影響更顯著，法向網(wǎng)格方面只有第一層網(wǎng)格高度對(duì)熱流影響巨大，之后的法向網(wǎng)格如何分布對(duì)熱流并無顯著影響。圖9為本文針對(duì)稀薄流區(qū)的氣動(dòng)熱問題開展的DSMC方法研究。基于非結(jié)構(gòu)網(wǎng)格及分布式并行計(jì)算技術(shù)，對(duì)X-37B外形在稀薄流區(qū)的氣動(dòng)熱問題進(jìn)行了求解，得到了飛行器表面熱流密度分布特性。

圖8 鈍雙錐表面流向熱流分布[66](Ma=9.86，α=10°，雷諾數(shù)Re=2.115×105)Fig.8 Heat flux distribution on surface of blunt biconical[66](Ma=9.86， α=10°， Reynolds number Re=2.115×105)

圖9 X-37B表面熱流密度分布(α=25°，h=90 km, N2∶O2∶O=0.788∶0.209∶0.003，V∞=7 km/s)Fig.9 Heat flux density distribution on surface of X-37B(α=25°， h=90 km，N2∶O2∶O=0.788∶0.209∶0.003, V∞=7 km/s)

3.3 基于Prandtl邊界層理論的數(shù)值與工程相結(jié)合的方法

與純數(shù)值解法相比，該方法對(duì)計(jì)算機(jī)要求并不是很苛刻，可節(jié)省大量計(jì)算時(shí)間和計(jì)算資源，同時(shí)又能提供比較準(zhǔn)確的氣動(dòng)熱解。而比起純工程算法，其適用范圍更廣，能應(yīng)用于復(fù)雜外形。

此方法中通常采用修正牛頓理論、活塞理論和Euler方程等獲得邊界層外緣參數(shù)分布。邊界層內(nèi)求解熱流密度主要基于局部相似法和參考焓法。由于邊界層方程是拋物型方程，沿著物面法向的信息傳播速度要快得多，因此上游信息的影響很快就衰減了。這樣，決定某個(gè)位置的邊界層特性主要是該處的壁面條件和邊界層外緣條件。對(duì)于可壓縮流，工程實(shí)踐證明，對(duì)不可壓流的摩阻公式進(jìn)行修正適用于絕大多數(shù)可壓縮流，其中應(yīng)用最廣的是?？颂?Eckert)參考溫度(參考焓)法。

呂麗麗等[67]通過求解三維Euler方程獲得邊界層外緣參數(shù)，利用局部相似性解的方法計(jì)算了鈍錐和純雙錐有迎角的再入的表面熱流，并與國(guó)內(nèi)外文獻(xiàn)的Navier-Stokes方程數(shù)值計(jì)算結(jié)果和風(fēng)洞試驗(yàn)結(jié)果進(jìn)行了比較，3種方法的結(jié)果吻合很好。Ji和Wang[68]通過求解Euler方程獲得邊界層外緣參數(shù)，邊界層內(nèi)采用平板模型，同時(shí)耦合了結(jié)構(gòu)傳熱，考慮了化學(xué)反應(yīng)效應(yīng)及轉(zhuǎn)捩模型，對(duì)復(fù)雜帶翼飛行器進(jìn)行了計(jì)算，結(jié)果表明，考慮化學(xué)非平衡后得到的熱流密度值減小(見圖10，圖中Tsurf為表面溫度)。

圖10 帶翼飛行器外表面溫度分布(Ma=5,α=6°, t=0～1 000 s)Fig.10 Temperature distribution on outer surface of winged aircraft(Ma=5,α=6°, t=0-1 000 s)

圖11 多體氣動(dòng)加熱(Ma=5,α=3°, h=40 km, time=150 s)Fig.11 Multi-body aerodynamics heating(Ma=5, α=3°, h=40 km, time=150 s)

圖11和圖12為本文對(duì)高超聲速組合體及組合體分離過程的氣動(dòng)加熱數(shù)值模擬方法進(jìn)行的研究(圖中Qtotal為總熱流)。圖11給出了2個(gè)飛行器的組合體在Ma=5、α=3°、h=40 km高空飛行150 s時(shí)的表面溫度與熱流分布(圖11(a)為表面溫度分布，圖11(b)為熱流分布)，其防熱結(jié)構(gòu)為2 mm的鈦及20 mm的二氧化硅。圖12給出了組合體在第150 s開始分離、釋放掛載3 s內(nèi)的載機(jī)及掛載的溫度與熱流分布，研究了高超聲速多體分離過程中的氣動(dòng)加熱干擾特性。

圖12 分離過程中的氣動(dòng)加熱(Ma=5, α=3°, h=40 km)Fig.12 Aerodynamic heating during separation(Ma=5,α=3°, h=40 km)

3.4 小 結(jié)

高超聲速飛行器的迅速發(fā)展對(duì)氣動(dòng)熱預(yù)測(cè)精度的要求越來越高，雖然經(jīng)過多年的努力，在氣動(dòng)熱預(yù)測(cè)方面有了一定的進(jìn)步，但還存在著大量的問題亟待解決，主要有：

1)激波交匯區(qū)、縫隙和凸起等局部氣動(dòng)加熱及分離區(qū)氣動(dòng)加熱的計(jì)算方法。

2)純數(shù)值氣動(dòng)熱計(jì)算的效率及精度方面，新的加速收斂措施及收斂判據(jù)。

4 氣動(dòng)力/熱/結(jié)構(gòu)耦合問題數(shù)值求解方法

氣動(dòng)熱彈性力學(xué)主要研究氣動(dòng)力、氣動(dòng)熱與飛行器的力學(xué)行為及其耦合效應(yīng)，是流體力學(xué)、傳熱學(xué)和結(jié)構(gòu)(動(dòng))力學(xué)的一門交叉學(xué)科，目前是國(guó)內(nèi)外在高超聲速領(lǐng)域的另一個(gè)研究熱點(diǎn)。

圖13 氣動(dòng)熱彈性問題耦合關(guān)系Fig.13 Coulping relationship of aerothermoelasticity problem

流體與結(jié)構(gòu)的耦合求解模型是氣動(dòng)力/熱/結(jié)構(gòu)耦合數(shù)值求解方法研究的重點(diǎn)之一。圖13[69]表示的是氣動(dòng)加熱形成的熱環(huán)境與構(gòu)成氣動(dòng)彈性問題的氣動(dòng)力、慣性力和彈性力形成的強(qiáng)弱不同的耦合關(guān)系。

目前采用的耦合模型[70]一般認(rèn)為熱傳導(dǎo)的特征時(shí)間遠(yuǎn)大于氣動(dòng)彈性系統(tǒng)的特征時(shí)間，因此重點(diǎn)考慮氣動(dòng)熱對(duì)結(jié)構(gòu)彈性的影響。該影響主要表現(xiàn)為兩方面：一是材料力學(xué)性能隨溫度的升高而變化；二是氣動(dòng)加熱是結(jié)構(gòu)產(chǎn)生額外的熱應(yīng)力，從而改變結(jié)構(gòu)靜、動(dòng)力學(xué)行為。

流體與結(jié)構(gòu)的耦合求解方法是氣動(dòng)力/熱/結(jié)構(gòu)耦合數(shù)值求解方法研究的另一個(gè)重點(diǎn)。氣動(dòng)熱彈性問題是一個(gè)多學(xué)科問題，一次性完全求解這樣的多學(xué)科耦合問題在相當(dāng)長(zhǎng)的時(shí)間內(nèi)難以在工程中實(shí)現(xiàn)。目前一般采用松耦合方法求解或者分層求解。美國(guó)20世紀(jì)90年代初的NASP計(jì)劃，90年代末期的Hyper-X等項(xiàng)目均開展了高超聲速飛行器氣動(dòng)熱彈性的研究，建立了高超聲速飛行器氣動(dòng)彈性問題較為成熟的工程分析框架，取得了大量的研究成果。Loehner等[71]提出了以CFD、計(jì)算結(jié)構(gòu)動(dòng)力學(xué)(CSD)和計(jì)算熱動(dòng)力學(xué)(CTD)的軟件解決氣動(dòng)、傳熱和結(jié)構(gòu)之間的耦合問題的分層求解思路。Thornton和Dechaumphai[72]用有限元法把氣動(dòng)、傳熱和結(jié)構(gòu)分析集成綜合代碼。對(duì)不銹鋼翼段的計(jì)算結(jié)果表明，在高超聲速下，由于氣動(dòng)熱和動(dòng)壓導(dǎo)致的翼段變形使氣流產(chǎn)生了激波、膨脹波以及回流區(qū)，加熱速率分布被嚴(yán)重改變了。國(guó)內(nèi)在此方面的研究尚處于起步階段，對(duì)于僅考慮氣動(dòng)加熱對(duì)于結(jié)構(gòu)彈性力部分的影響，吳志剛等[73]采用了分層求解的分析方法，即先計(jì)算高超聲速定常氣動(dòng)力與氣動(dòng)熱，然后進(jìn)行結(jié)構(gòu)傳熱分析，得到溫度分布，最好進(jìn)行顫振分析。張偉偉等[74]采用松耦合方法建立了氣動(dòng)彈性仿真模型，將問題分成了氣動(dòng)熱計(jì)算、結(jié)構(gòu)的溫度場(chǎng)分布計(jì)算、熱結(jié)構(gòu)計(jì)算和氣動(dòng)彈性計(jì)算，在時(shí)間域內(nèi)實(shí)現(xiàn)了高超聲速熱氣動(dòng)彈性的仿真。圖14為本文針對(duì)機(jī)翼氣動(dòng)力、氣動(dòng)熱和結(jié)構(gòu)耦合效應(yīng)開展的研究，計(jì)算給出了溫度及結(jié)構(gòu)熱應(yīng)力分布特性，分析了力熱結(jié)構(gòu)耦合特性在機(jī)翼流動(dòng)中的影響特性。圖中：F1、F2、F3和F4分別為機(jī)翼前四階模態(tài)的頻率。

目前的耦合模型忽略了結(jié)構(gòu)變形和振動(dòng)對(duì)熱環(huán)境的影響及結(jié)構(gòu)熱生成與彈性變形之間的熱力學(xué)耦合，因此精確的氣動(dòng)力/熱/結(jié)構(gòu)耦合分析模型還有待建立?；诰€性有限元、全位勢(shì)流和線性熱傳導(dǎo)的耦合方法較為成熟，在工程分析中得到了廣泛應(yīng)用，而基于非線性有限元，Euler或Navier-Stokes方程和非線性熱傳導(dǎo)方程的耦合分析方程則是當(dāng)前熱氣動(dòng)彈性力學(xué)的重要研究方向之一。

圖14 典型機(jī)翼氣動(dòng)力/熱/結(jié)構(gòu)耦合計(jì)算Fig.14 Coupled flow/thermal/structural analysis of typical wing

5 總 結(jié)

在過去的幾十年間，高超聲速?gòu)?fù)雜氣動(dòng)問題的數(shù)值求解技術(shù)有了長(zhǎng)足的進(jìn)展，取得了顯著的研究成果，并且部分成果已成功的應(yīng)用于工程設(shè)計(jì)。近些年來，在高超聲速飛行器設(shè)計(jì)方面，各科技大國(guó)越來越關(guān)注臨近空間(地球大氣層20～100 km范圍內(nèi))的開發(fā)和利用，高超聲速飛行器技術(shù)需求不斷提高，越來越多的涉及交叉學(xué)科與多場(chǎng)耦合的氣動(dòng)問題被廣泛關(guān)注，如氣動(dòng)/熱/結(jié)構(gòu)耦合的高超聲速熱氣動(dòng)彈性力學(xué)、氣動(dòng)/化學(xué)反應(yīng)/電磁場(chǎng)耦合的高超聲速電磁流體力學(xué)與氣動(dòng)光學(xué)、氣/液/化學(xué)反應(yīng)相耦合的超聲速燃燒流動(dòng)等。這些問題的有效解決，將在高超聲速飛行器先進(jìn)布局設(shè)計(jì)、高效動(dòng)力系統(tǒng)等技術(shù)方面獲得重大突破。

結(jié)合當(dāng)前技術(shù)需求與發(fā)展，在以下幾個(gè)方面值得關(guān)注：

1)跨流域、全速域以及多場(chǎng)耦合氣動(dòng)物理問題的認(rèn)知。

2)復(fù)雜氣動(dòng)物理問題的數(shù)學(xué)建模，包括控制方程、壁面條件、湍流模型、化學(xué)反應(yīng)模型以及氣/液/固/電磁多場(chǎng)耦合模型等。

3)復(fù)雜流動(dòng)高效準(zhǔn)確的求解技術(shù)與結(jié)果后處理，包括先進(jìn)網(wǎng)格技術(shù)、時(shí)/空離散格式、多場(chǎng)耦合算法以及復(fù)雜問題的結(jié)果分析等。

4)數(shù)值求解技術(shù)與風(fēng)洞試驗(yàn)技術(shù)的相互校驗(yàn)，及其在工程設(shè)計(jì)中的應(yīng)用。

相信隨著人們對(duì)物理問題的認(rèn)知越來越深、求解技術(shù)的越來越精細(xì)，以及高性能計(jì)算設(shè)備的迅速發(fā)展，高超聲速?gòu)?fù)雜流動(dòng)問題的數(shù)值模擬技術(shù)會(huì)取得更多的具有自主知識(shí)產(chǎn)權(quán)的研究成果。

致 謝

感謝國(guó)家自然科學(xué)基金、航空科學(xué)基金、南京航空航天大學(xué)創(chuàng)新基金及國(guó)家“863”計(jì)劃等對(duì)本文研究工作的大力支持與資助。

[1] Anderson J D. Hypersonic and high temperature gas dynamics[M]. Reston: AIAA, 2000: 13-23.

[2] Bian Y G, Xu L G. Aerothermodynamics[M]. Hefei: Press of University of Science and Technology of China, 1997: 2-8 (in Chinese). 卞蔭貴, 徐立功. 氣動(dòng)熱力學(xué)[M]. 合肥: 中國(guó)科學(xué)技術(shù)大學(xué)出版社, 1997: 2-8.

[3] Jameson A, Schmidt W, Turkel E. Numerical solutions of the Euler equations by finite volume methods using Runge-Kutta time-stepping schemes, AIAA-1981-1259[R]. Reston: AIAA, 1981.

[4] Mavriplis D J, Jameson A. Multigrid solution of the Navier-Stokes equations on triangular meshes[J]. AIAA Journal, 1990, 28(8): 1415-1425.

[5] Briley W R, McDonald H. Solution of the multidimensional compressible Navier-Stokes equations by a generalized implicit method[J]. Journal of Computational Physics, 1977, 24(4): 372-397.

[6] Jameson A, Turkel E. Implicit schemes and LU-decompositions[J]. Mathematics of Computation, 1981, 37(156): 385-397.

[7] Yoon S, Jameson A. Lower-upper symmetric-Gauss-Seidel method for the Euler and Navier-Stokes equations[J]. AIAA Journal, 1988, 26(9): 1025-1026.

[8] Saad Y, Schultz M H. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[J]. SIAM Journal on Scientific and Statistical Computing, 1986, 7(3): 856-869.

[9] van Leer B. Flux-vector splitting for the Euler equations[J]. Lecture Notes on Physics, 1982, 170: 507-512.

[10] Roe P L. Approximate Riemann solvers, parameter vectors, and difference schemes[J]. Journal of Computational Physics, 1981, 43(2): 357-372.

[11] Liou M S, Steffen C J. A new flux splitting scheme[J]. Journal of Computational Physics, 1993, 107(1): 23-39.

[12] Liou M S. Ten years in the making: AUSM-family, AIAA-2001-2521[R]. Reston: AIAA, 2001.

[13] Jameson A. Analysis and design of numerical schemes for gas dynamics, 1: artificial diffusion, upwind biasing, limiters and their effect on accuracy and multigrid convergence[J]. International Journal of Computational Fluid Dynamics, 1995, 4(3): 171-218.

[14] Zhang H X. Non-oscillatory and non-free-parameter dissipation difference scheme[J]. Acta Aerodynamica Sinica, 1988, 6(2): 143-165 (in Chinese). 張涵信. 無波動(dòng), 無自由參數(shù)的耗散差分格式[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 1988, 6(2): 143-165.

[15] Deardorff J W. A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers[J]. Journal of Fluid Mechanics, 1970, 41(2): 453-480.

[16] Moin P, Mahesh K. Direct numerical simulation: a tool in turbulence research[J]. Annual Review of Fluid Mechanics, 1998, 30(1): 539-578.

[17] Cebeci T, Smith A M O. Analysis of turbulent boundary layers[M].New York: Academic Press, Inc., 1974: 164-187.

[18] Baldwin B S, Lomax H. Thin layer approximation and algebraic model for separated turbulent flows, AIAA-1978-0257[R]. Reston: AIAA, 1978.

[19] Spalart P R, Allmaras S R. A one equation turbulence model for aerodinamic flows, AIAA-1992-0439[R]. Reston: AIAA, 1992.

[20] Baldwin B S, Barth T J. A one-equation turbulence transport model for high Reynolds number wall-bounded flows, NASA TM 102847[R]. Washington, D.C.: NASA, 1990.

[21] Jones W P, Launder B E. The prediction of laminarization with a two-equation model of turbulence[J]. International Journal of Heat and Mass Transfer, 1972, 15(2): 301-314.

[22] Liu J Y. An improverd SST turbulence model for hypersonic flows[J]. Acta Aeronautica et Astronautica Sinica,2012, 33(12): 2192-2201 (in Chinese). 劉景源. SST 湍流模型在高超聲速繞流中的改進(jìn)[J]. 航空學(xué)報(bào), 2012, 33(12): 2192-2201.

[23] Wilcox D C. Reassessment of the scale-determining equation for advanced turbulence models[J]. AIAA Journal, 1988, 26(11): 1299-1310.

[24] Menter F R. Two-equation eddy-viscosity turbulence models for engineering applications[J]. AIAA Journal, 1994, 32(8): 1598-1605.

[25] Blottner F G. Nonequilibrium laminar boundary-layer flow of ionized air[J]. AIAA Journal, 1964, 2(11): 1921-1927.

[26] Zhang X H, Wu Y Z, Wang J F. Numerical simulation of thermo-chemical non-equilibrium hypersonic flows using HLLE+ scheme[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2011, 43(2): 154-158 (in Chinese). 張向洪, 伍貽兆, 王江峰. HLLE+ 格式在高超聲速熱化學(xué)非平衡流場(chǎng)中的應(yīng)用[J]. 南京航空航天大學(xué)學(xué)報(bào), 2011, 43(2): 154-158.

[27] Liu C, Wang J F, Wu Y Z. Numerical simulation of multi-component reacting flows using AUFS scheme[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2006, 40(2): 191-194 (in Chinese). 劉晨, 王江峰， 伍貽兆. AUFS格式在多組分反應(yīng)流場(chǎng)模擬中的應(yīng)用[J]. 南京航空航天大學(xué)學(xué)報(bào), 2006, 40(2): 191-194.

[28] Liu J, Liu W, Zeng M, et al. Numerical simulation of 3D hypersonic thermochemical non-equilibrium flow[J]. Chinese Journal of Theoretical and Applied Mechanics， 2003， 35(6)： 730-734 (in Chinese). 柳軍, 劉偉, 曾明, 等. 高超聲速三維熱化學(xué)非平衡流場(chǎng)的數(shù)值模擬[J]. 力學(xué)學(xué)報(bào)， 2003， 35(6)： 730-734.

[29] Herzberg G. Molecular spectra and molecular structure, Vol.1[M]. New York: Reitell Press, 1957: 13-27.

[30] Moore C E. Atomic energy levels as derived from the analyses of optical spectra[M]. Washington, D. C.: US Government Printing Office, 1971: 65-98.

[31] Hansen C F, Heims S P. A review of the thermodynamic, transport, and chemical reaction rate properties of high-temperature air[M]. [s.l.]: National Advisory Committee for Aeronautics, 1958： 105-112.

[32] Gupta R N, Yos J M, Thompson R A. A review of reaction rates and thermodynamic and transport properties for the 11-species air model for chemical and thermal nonequilibrium calculations to 30 000 K, NASA-TM-101528[R]. Hampton, VA： Langley Research Center, 1989.

[33] Eklund D R, Stouffer S D, Northam G B. Study of a supersonic combustor employing swept ramp fuel injectors[J]. Journal of Propulsion and Power, 1997, 13(6): 697-704.

[34] Clutter J K, Mikolaitis D W, Shyy W. Effect of reaction mechanism in shock-induced combustion simulations, AIAA-1998-0274[R]. Reston: AIAA, 1998.

[35] Ingram D M, Jiang B, Causon D M. On methane combustion in a nozzle geometry using a reduced reaction kinetics model, AIAA-1996-0820[R]. Reston: AIAA, 1996.

[36] Liu C, Wang J F, Wu Y Z. Study on the abnormity of the axial-symmetric boundary in simulation of shock-induced combustion[J]. Journal of Aerospace Power, 2009, 24(6): 1219-1228 (in Chinese). 劉晨, 王江峰, 伍貽兆. 激波誘導(dǎo)燃燒模擬中軸對(duì)稱邊界條件數(shù)值異常研究[J]. 航空動(dòng)力學(xué)報(bào), 2009, 24(6): 1219-1228.

[37] Bird G A. Approach to translational equilibrium in a rigid sphere gas[J]. Physics of Fluids, 1963(6): 1518-1519.

[38] Bird G A. Application of the DSMC method to the full shuttle geometry, AIAA-1990-1692[R]. Reston: AIAA, 1990.

[39] Shen Q. Aerodynamics of rarefied gases[M]. Beijing: National Defense Industry Press, 2003: 17-23(in Chinese). 沈青. 稀薄氣體動(dòng)力學(xué)[M]. 北京:國(guó)防工業(yè)出版社，2003: 17-23.

[40] Fan J. Rarefied gas dynamics: advances and applications[J]. Advance in Mechanics, 2013, 43(2): 185-201 (in Chinese). 樊菁. 稀薄氣體動(dòng)力學(xué):進(jìn)展與應(yīng)用[J]. 力學(xué)進(jìn)展, 2013, 43(2): 185-201.

[41] Fan J, Liu H L, Jiang J Z, et al. Analysis and simulation of discharging residual rocket propellants in orbit[J]. Chinese Journal of Theoretical and Applied Mechanics, 2004, 36(2): 129-139 (in Chinese). 樊菁, 劉宏立, 蔣建政, 等. 火箭剩余推進(jìn)劑排放過程的分析與模擬[J]. 力學(xué)學(xué)報(bào), 2004, 36(2): 129-139.

[42] Fan J. A generalized soft sphere model for Monte Carlo simulations[J]. Physics of Fluids, 2002, 14: 4399-4405.

[43] Li Z H, Li Z H, Li H Y, et al. Research on CFD/DSMC hybrid numerical method in rarefied flows[J]. Acta Aerodynamica Sinica, 2013, 31(3): 282-287 (in Chinese). 李中華, 李志輝, 李海燕, 等. 過渡流區(qū)Navier-Stokes/DSMC耦合計(jì)算研究[J].空氣動(dòng)力學(xué)學(xué)報(bào), 2013, 31(3):282-287.

[44] Wang X D. DSMC method on unstructured grids for hypersonic rarefied gas flow and its parallelization[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2006 (in Chinese). 王學(xué)德. 高超聲速稀薄氣流非結(jié)構(gòu)網(wǎng)格 DSMC 及并行算法研究[D]. 南京： 南京航空航天大學(xué), 2006.

[45] Wang X D, Wu Y Z, Xiao J. Unstructured DSMC method and application in three dimensional thermochemical nonequilibrium flow[J]. Journal of Astronautics, 2006， 27(12): 126-131 (in Chinese). 王學(xué)德, 伍貽兆, 夏健. 三維熱化學(xué)非平衡流動(dòng)非結(jié)構(gòu)網(wǎng)格 DSMC 方法及其應(yīng)用[J]. 宇航學(xué)報(bào), 2006，27(12): 126-131.

[46] Pan Y, Wang J F, Wu Y Z. Numerical simulation of hypersonic viscous MHD flows based on unstructured meshes[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2007, 39(5): 555-559 (in Chinese). 潘勇, 王江峰, 伍貽兆. 非結(jié)構(gòu)網(wǎng)格的高超聲速粘性MHD流場(chǎng)數(shù)值模擬[J]. 南京航空航天大學(xué)學(xué)報(bào), 2007, 39(5): 555-559.

[47] Zhang X H. Numerical simulation for hypersonic flowfield with electromagnetic interference[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2013 (in Chinese). 張向洪. 高超聲速流場(chǎng)電磁干擾數(shù)值模擬研究[D]. 南京: 南京航空航天大學(xué), 2013.

[48] Brackbill J U, Barnes D C. The effect of nonzero on the numerical solution of the magnetohydrodynamic equations[J]. Journal of Computational Physics, 1980, 35: 426-430.

[49] Evans C R, Hawley J F. Simulation of magnetohydrodynamic flows: a constrained transport method[J]. The Astrophysical Journal, 1988, 332: 659-677.

[50] Powell K G, Roe P L, Myong R S. An upwind scheme for magnetohydrodynamics, AIAA-1995-1704[R]. Reston: AIAA, 1995.

[51] Dedner A, Kemm F, Kroner D, et al. Hyperbolic divergence cleaning for the MHD equations[J]. Journal of Computational Physics, 2002, 175: 645-673.

[52] Pan Y. Numerical methods for hypersonic flowfield with magnetic interference[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2007 (in Chinese). 潘勇. 高超聲速流場(chǎng)磁場(chǎng)干擾效應(yīng)數(shù)值模擬方法研究[D]. 南京: 南京航空航天大學(xué), 2007.

[53] Tian Z Y, Zhang K P, Ding G H, et al. Spurious magnetic field divergence cleaning in magnetohydrodynamic simulation[J]. Chinese Journal of Computational Physics, 2009, 26(1): 78-86 (in Chinese). 田正雨， 張康平， 丁國(guó)昊，等. MHD數(shù)值模擬中清除偽磁場(chǎng)散度方法[J]. 計(jì)算物理, 2009, 26(1): 78-86.

[54] de Jarnette F R, Hamilton H H. Inviscid surface streamlines and heat transfer on shuttle-type configurations[J]. Journal of Spacecraft and Rockets, 1973, 10(5): 314-321.

[55] Engel C D, Praharaj S C. MINIVER upgrade for the AVID system. Volume 1: LANMIN user’s manual, NASA CR-172212[R]. Washington, D.C.: NASA, 1983.

[56] Zoby E V, Simmonds A L. Engineering flowfield method with angle-of-attack applications[J]. Journal of Spacecraft and Rockets, 1985, 22(4): 398-404.

[57] Li J L, Tang Q G, Huo L, et al. The rapid engineering aero-heating calculation method for complex shaped hypersonic vehicles[J]. Journal of National University of Defense Technology, 2012, 34(6): 89-93 (in Chinese). 李建林, 唐乾剛, 霍霖, 等. 復(fù)雜外形高超聲速飛行器氣動(dòng)熱快速工程估算[J]. 國(guó)防科學(xué)技術(shù)大學(xué)學(xué)報(bào), 2012, 34(6): 89-93.

[58] Murray A L, Lewis C H. Hypersonic three-dimensional viscous shock-layer flows over blunt bodies[J]. AIAA Journal, 1978, 16(12): 1279-1286.

[59] Helliwell W S, Dickinson R P, Lubard S C. Viscous flow over arbitrary geometries at high angle of attack[J]. AIAA Journal, 1981, 19(2): 191-197.

[60] Tannehill J C, Buelow P E, Ievalts J O, et al. Three-dimensional upwind parabolized Navier-Stokes code for real gasflows[J]. Journal of Spacecraft and Rockets, 1990, 27(2): 150-159.

[61] Gnoffo P. Upwind-biased, point-implicit relaxation strategies for viscous, hypersonic flows, AIAA-1989-1972[R]. Reston: AIAA, 1989.

[62] Huang F, Zhang L, Cheng X L, et al. Effects of continuum breakdown on aerothermodynamics[J]. Journal of Astronautics, 2012, 33(2): 153-159 (in Chinese). 黃飛, 張亮, 程曉麗, 等. 稀薄氣體效應(yīng)對(duì)尖前緣氣動(dòng)熱特性的影響研究[J]. 宇航學(xué)報(bào), 2012, 33(2): 153-159.

[63] Hamilton H H, Greene F A, DeJarnette F R. Approximate method for calculating heating rates on three-dimensional vehicles[J]. Journal of Apacecraft and Rockets, 1994, 31(3): 345-354.

[64] Pan S, Feng D H, Ding G H, et al. Grid dependency and convergence of hypersonic aerothermal simulation[J]. Acta Aeronautica et Astronautica Sinica, 2010, 31(3): 493-499 (in Chinese). 潘沙, 馮定華, 丁國(guó)昊, 等. 氣動(dòng)熱數(shù)值模擬中的網(wǎng)格相關(guān)性及收斂[J]. 航空學(xué)報(bào), 2010, 31(3): 493-499.

[65] Yan C, Yu J J, Li J Z, Scheme effect and grid dependency in CFD compulations of heating transfer[J]. Acta Aerodynamica Sinica, 2006, 24(1): 125-130 (in Chinese). 閻超, 禹建軍, 李君哲. 熱流 CFD 計(jì)算中格式和網(wǎng)格效應(yīng)若干問題研究[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2006, 24(1): 125-130.

[66] Li J Z, Yan C, Ke L, et al. Research on scheme effect of computational fluid dynamics in aerothermal[J]. Journal of Beijing University of Aeronautics and Astronautics, 2003, 29(11): 1022-1025 (in Chinese). 李君哲, 閻超, 柯倫, 等. 氣動(dòng)熱CFD計(jì)算的格式效應(yīng)研究[J]. 北京航空航天大學(xué)學(xué)報(bào), 2003, 29(11): 1022-1025.

[67] Lv L L, Zhang W W, Ye Z Y. Predicting heating distributions for hypersonic reentry bodies[J]. Chinese Journal of Applied Mechanics, 2006, 23(2): 259-262(in Chinese). 呂麗麗, 張偉偉, 葉正寅. 高超聲速再入體表面熱流計(jì)算[J]. 應(yīng)用力學(xué)學(xué)報(bào), 2006, 23(2): 259-262.

[68] Ji W D, Wang J F. Calculating method of aerodynamic heating for hypersonic aircrafts[J]. Transactions of Nanjing University of Aeronautics & Astronautics, 2013, 30(3): 327-342.

[69] Cai G B, Xu D J. Hypersonic vehicle technology[M]. Beijing: Science Press, 2012: 34-45 (in Chinese). 蔡國(guó)飆, 徐大軍. 高超聲速飛行器技術(shù)[M]. 北京: 科學(xué)出版社, 2012: 34-45.

[70] Yang C, Xu Y, Xie C C. Review of studies on aeroelasticity of hypersonic vehicles[J]. Acta Aeronautica et Astronautica Sinica, 2010, 31(1): 1-11 (in Chinese). 楊超, 許赟, 謝長(zhǎng)川. 高超聲速飛行器氣動(dòng)彈性力學(xué)研究綜述[J]. 航空學(xué)報(bào), 2010, 31(1): 1-11.

[71] Loehner R, Yang C, Cerbal J, er al. Fluid-structure-thermal interaction using a loose coupling algorithm and adaptive unstructured grids, AIAA-1998-2419[R]. Reston: AIAA, 1998.

[72] Thornton E, Dechaumphai P. Coupled flow, thermal and structural analysis of aerodynamically heated panels[J]. Journal of Aircraft, 1988, 25(11): 1052-1058.

[73] Wu Z G, Hui J P, Yang C. Hypersonic aerothermoelastic analysis of wings[J]. Journal of Beijing University of Aeronautics and Astronautics, 2005, 31(3): 270-273 (in Chinese). 吳志剛, 惠俊鵬, 楊超. 高超聲速下翼面的熱顫振工程分析[J]. 北京航空航天大學(xué)學(xué)報(bào), 2005, 31(3): 270-273.

[74] Zhang W W, Xia W, Ye Z Y. A numerical method for hypersonic aerothermoelasticity[J]. Engineering Mechanics, 2006, 23(2): 41-46 (in Chinese). 張偉偉, 夏巍, 葉正寅. 一種高超音速熱氣動(dòng)彈性數(shù)值研究方法[J]. 工程力學(xué), 2006, 23(2): 41-46.

Tel: 025-84891231

E-mail: wangjf@nuaa.edu.cn

伍貽兆 男， 博士， 教授， 博士生導(dǎo)師。主要研究方向：計(jì)算流體力學(xué)，飛行器氣動(dòng)布局設(shè)計(jì)。

Tel: 025-84891231

E-mail: wyzao@nuaa.edu.cn

*Corresponding author. Tel.： 025-84891231 E-mail: wangjf@nuaa.edu.cn

Progress in numerical simulation techniques of hypersonic aerodynamic problems

WANG Jiangfeng*, WU Yizhao, JI Weidong, FAN Xiaofeng, ZHAO Faming, LYU Zhenjun

CollegeofAerospaceEngineering,NanjingUniversityofAeronauticsandAstronautics,Nanjing210016,China

Hypersonic flow field has complex flow characteristics in which real gas effects, magnetic fluid interference effects and fluid/thermal/structural coupling effects have an important impact on the aerodynamic force. They extend fluid dynamics to molecular dynamics, electromagnetic fluid dynamics, fluid/structure interaction and other interdisciplinary fields, which have brought great challenges to the numerical simulation methods. Aimed at hot issues of hypersonic aerodynamic force and aerodynamic heat, high-temperature effects, low-density flow effect, magnetic fluid interference effect and fluid/thermal/structural coupling effect have been significantly emphasized. Several examples and the corresponding numerical solution techniques are given in this paper. Three methods of aerodynamic heating are compared, i.e., pure engineering method, pure numerical method and Prandtl boundary layer theory-based method. For fluid/thermal/structural coupling problem, analyses are carried out in two aspects, i.e., coupling model and coupling calculation method. Finally, several problems of numerical simulation technologies which need to be emphasized in the future are figured out.

hypersonic; numerical simulation; aerodynamic heating; magnetohydrodynamics; aerothermoelasticity; fluid/thermal/structural coupling

2014-08-15； Revised： 2014-08-29； Accepted： 2014-09-02； Published online： 2014-09-24 14:55

National High-tech Research and Development Program of China

2014-08-15; 退修日期： 2014-08-29; 錄用日期： 2014-09-02; 網(wǎng)絡(luò)出版時(shí)間： 2014-09-24 14：55

www.cnki.net/kcms/detail/10.7527/S1000-6893.2014.0214.html

國(guó)家“863”計(jì)劃

Wang J F, Wu Y Z, Ji W D, et al. Progress in numerical simulation techniques of hypersonic aerodynamic problems[J].Acta Aeronautica et Astronautica Sinica, 2015, 36(1): 159-175. 王江峰， 伍貽兆， 季衛(wèi)棟， 等． 高超聲速?gòu)?fù)雜氣動(dòng)問題數(shù)值方法研究進(jìn)展[J]．航空學(xué)報(bào), 2015, 36(1): 159-175.

http://hkxb.buaa.edu.cn hkxb@buaa.edu.cn

10.7527/S1000-6893.2014.0214

V211.3

A

1000-6893(2015)01-0159-17

王江峰 男， 博士， 教授， 博士生導(dǎo)師。主要研究方向：高超聲速飛行器布局設(shè)計(jì)、燃燒流場(chǎng)數(shù)值模擬技術(shù)、氣動(dòng)加熱計(jì)算技術(shù)等。

*通訊作者.Tel.： 025-84891231 E-mail: wangjf@nuaa.edu.cn

URL： www.cnki.net/kcms/detail/10.7527/S1000-6893.2014.0214.html