唐志共， 張益榮， 陳堅強(qiáng)， 毛枚良， 張毅鋒， 劉化勇

中國空氣動力研究與發(fā)展中心， 綿陽 621000

更準(zhǔn)確、更精確、更高效
——高超聲速流動數(shù)值模擬研究進(jìn)展

唐志共， 張益榮， 陳堅強(qiáng)*， 毛枚良， 張毅鋒， 劉化勇

中國空氣動力研究與發(fā)展中心， 綿陽 621000

從準(zhǔn)度、精度和效率3方面回顧了近幾十年來高超聲速流動數(shù)值模擬研究的進(jìn)展。在物理模型方面，介紹了高超聲速數(shù)值模擬中高溫氣體效應(yīng)、稀薄氣體效應(yīng)以及湍流效應(yīng)的建模與模擬，基于雷諾平均Navier-Stokes(RANS)方程重點(diǎn)對現(xiàn)階段較為關(guān)注的高超聲速邊界層轉(zhuǎn)捩的模式理論研究進(jìn)行了介紹。在空間離散算法方面，主要介紹了高超聲速數(shù)值模擬中常用的二階精度迎風(fēng)格式以及高階精度格式的發(fā)展及其應(yīng)用。在時間推進(jìn)方面，主要回顧了隱式時間推進(jìn)方法的發(fā)展及其應(yīng)用。在誤差和不確定度估計方面，主要介紹了其概念、來源以及常用的分析方法，同時給出了迭代誤差估計、Richardson外插法以及敏感性導(dǎo)數(shù)方法等初步研究結(jié)果。最后，討論了高超聲速流動數(shù)值模擬中下一步需關(guān)注的問題。

高超聲速流動； 數(shù)值模擬； 物理建模； 數(shù)值方法； 誤差； 不確定度

計算流體力學(xué)(CFD)在自誕生至今的近百年時間里，隨著物理模型、數(shù)值算法和計算機(jī)軟硬件水平的不斷發(fā)展而迅速崛起，成為飛行器研制中氣動布局設(shè)計及數(shù)據(jù)庫建設(shè)的重要手段。與此同時，由于高超聲速飛行器具有的特殊軍事意義及經(jīng)濟(jì)、科技價值，世界各國紛紛將高超聲速飛行器的研制列入重點(diǎn)發(fā)展計劃[1]。以美國為例，其在技術(shù)發(fā)展水平和人員、資金投入力度等方面都走在世界最前列，先后提出了一系列高超聲速發(fā)展計劃，包括20世紀(jì)60年代的X-15驗證機(jī)，1986年的國家空天飛機(jī)計劃(NASP)，2001年提出的航空航天倡議(NAI)，2002年形成的高超聲速飛行計劃(HyFly)，2003年提出的FALCON(Force Application and Launch from CONtinental United States)計劃(HTV和HCV)，以及最近幾年發(fā)展迅速的X-37B和X-51A等，如圖1所示，在高超聲速吸氣式發(fā)動機(jī)、數(shù)值模擬物理化學(xué)建模、熱防護(hù)以及數(shù)據(jù)再利用等方面積累了豐富的經(jīng)驗，在諸多領(lǐng)域取得了舉世矚目的成就?？v觀這一系列高超聲速計劃的發(fā)展歷程，正是其中所暴露出來的問題，才不斷推動著高超聲速技術(shù)，特別是數(shù)值模擬技術(shù)的發(fā)展。例如，X-15在試驗飛行過程中由于嚴(yán)重的激波/激波和激波/邊界層干擾導(dǎo)致沖壓發(fā)動機(jī)整流罩的損毀，如圖2所示，使人們認(rèn)識到激波干擾問題的復(fù)雜性和重要性，而航天飛機(jī)機(jī)體襟翼俯仰力矩配平異常，更是直接推動了真實氣體效應(yīng)的研究，這些都促進(jìn)了數(shù)值模擬手段的發(fā)展。

隨著高超聲速飛行器的發(fā)展，其任務(wù)越來越復(fù)雜，飛行包線越來越寬，對飛行器機(jī)動性和防熱等都提出了更高要求。不管是利用CFD完成低馬赫數(shù)風(fēng)洞數(shù)據(jù)到高馬赫飛行數(shù)據(jù)的外推，還是直接利用CFD進(jìn)行高空高馬赫數(shù)飛行狀態(tài)的模擬，都對高超聲速數(shù)值模擬技術(shù)的可信度提出了更高要求。從實驗觀察到的物理現(xiàn)象到數(shù)學(xué)形式的物理模型，再到數(shù)值模擬的離散方程、邊界條件以及輔助模型等，CFD與真實情況之間存在諸多假設(shè)、簡化和近似。由此，準(zhǔn)度和精度的評估就成為可信度研究中最重要的兩個方面，對于高超聲速數(shù)值模擬，分別涉及到物理建模、邊界條件和輔助模型準(zhǔn)確性以及時間/空間離散精度等。國外對此開展了大量研究工作，如北約主持的一系列高超聲速流動數(shù)值模擬確認(rèn)研究[2-3]，美國國家航空航天局(NASA)把不確定度量化作為2014—2030年CFD研究路線圖中的重要部分[4]，提出到2030年將具備CFD誤差估計和不確定度傳播分析的能力。此外，隨著高超聲速飛行器執(zhí)行任務(wù)和幾何構(gòu)型越來越復(fù)雜，研究周期越來越短，對高超聲速流動數(shù)值模擬的求解效率提出了更高要求。對于實際應(yīng)用問題，加速收斂技術(shù)和并行計算效率是主要的影響因素，前者如多重網(wǎng)格、網(wǎng)格自適應(yīng)等常用技術(shù)并不適用于存在強(qiáng)間斷的高超聲速流動數(shù)值模擬中，后者依賴于計算機(jī)硬件水平和數(shù)值算法等，其中，隱式時間推進(jìn)方法由于優(yōu)越的魯棒性和高效性，得到了迅速發(fā)展及應(yīng)用。

圖1 高超聲速飛行器Fig.1 Hypersonic vehicle

圖2 X-15飛機(jī)Fig.2 X-15 aircraft

本文針對高超聲速流動模擬中的準(zhǔn)度、精度和效率問題，詳細(xì)介紹了近年來在物理模型、空間離散算法、時間推進(jìn)方法以及誤差和不確定度估計等方面的研究進(jìn)展，尤其對目前的熱點(diǎn)問題，如高超聲速邊界層轉(zhuǎn)捩預(yù)測、高階格式在復(fù)雜流動中的應(yīng)用以及CFD結(jié)果可信度評估等方面進(jìn)行了詳細(xì)介紹，并對下一步高超聲速流動數(shù)值模擬中擬關(guān)注的問題進(jìn)行了討論。

1 更準(zhǔn)確——物理模型研究進(jìn)展

1.1 主要的物理模型

高超聲速飛行器在飛行過程中將可能經(jīng)歷高溫真實氣體效應(yīng)、稀薄氣體效應(yīng)、黏性干擾、邊界層轉(zhuǎn)捩和分離以及熱輻射等復(fù)雜物理化學(xué)現(xiàn)象，如圖 3所示。

圖3 高超聲速飛行器典型流動特征Fig.3 Typical flow feature of a hypersonic vehicle

CFD如何真實復(fù)現(xiàn)高超聲速飛行，很大程度上取決于物理模型反映真實情況的程度。目前，工程應(yīng)用領(lǐng)域一般按高度劃分流域(通常對應(yīng)克努森數(shù))，如60 km以下連續(xù)流假設(shè)成立，該區(qū)域的數(shù)值模擬可基于連續(xù)流假設(shè)的Navier-Stokes方程組進(jìn)行，而60～90 km為滑移流區(qū)，可利用添加滑移邊界條件的Navier-Stokes方程組進(jìn)行計算，90 km以上連續(xù)假設(shè)不再成立，需要采用直接模擬蒙特卡羅(DSMC)等方法求解[5]。此外，近些年稀薄流到連續(xù)流統(tǒng)一算法的發(fā)展也較為迅速[6]，是富有希望的一種全流域求解策略。對于基于Navier-Stokes方程組求解的高超聲速流動數(shù)值模擬，高溫真實氣體效應(yīng)建模涉及到多組分化學(xué)反應(yīng)混合氣體控制方程(包括總質(zhì)量守恒、總動量守恒、總能量守恒、各組分守恒、不同內(nèi)部分子模式(轉(zhuǎn)動、振動以及電子激發(fā)，包括與平動之間的模式間相互轉(zhuǎn)換)之間的能量守恒以及熱狀態(tài)方程)以及與這些物理現(xiàn)象相對應(yīng)的邊界條件(包括允許速度和溫度滑移的氣/固界面動量輸運(yùn)，允許組分生成和消耗的氣/固界面質(zhì)量輸運(yùn)等)等。目前，國內(nèi)外主流高超聲速CFD軟件，如GASP、LAURA和VULCAN等，在模擬高溫真實氣體效應(yīng)時，普遍根據(jù)求解問題的特點(diǎn)，不同層次地采用帶5組分、7組分或11組分空氣化學(xué)反應(yīng)源項Navier-Stokes控制方程組，考慮平動、轉(zhuǎn)動、振動和電子能量激發(fā)模型，以及可能的電離效應(yīng)等影響。實驗和數(shù)值模擬的研究表明，高溫真實氣體效應(yīng)對激波脫體距離、駐點(diǎn)熱流峰值、表面摩擦阻力分布以及飛行器氣動力和力矩產(chǎn)生較大影響，但從文獻(xiàn)調(diào)研情況可知，目前普遍采用的高溫真實氣體效應(yīng)物理模型框架形成于20世紀(jì)90年代[7]，近些年卻鮮有突破性進(jìn)展。圖 4給出了真實氣體和完全氣體計算得到的頭部弓形激波脫體距離的差異。圖 5是美國航天飛機(jī)利用不同物理模型(包括GASPex軟件求解的完全氣體(Ideal Gas)、五組分空氣化學(xué)反應(yīng)模型(React. 5-sp Air)、FELISA軟件求解的平衡氣體模型(Eq. Air))的計算結(jié)果和OADB(Operational Aerodynamic Data Book)數(shù)據(jù)的比較，由圖可知真實氣體效應(yīng)對俯仰力矩系數(shù)Cm的影響較大[8]。

圖4 真實氣體和完全氣體模型計算馬赫數(shù)等值線比較Fig.4 Comparison of Mach number contours between real gas and ideal gas model

圖5 俯仰力矩系數(shù)計算結(jié)果和OADB數(shù)據(jù)比較(Ma=15)[8]Fig.5 Comparison of pitching moment coefficients for calculation results and OADB data (Ma=15)[8]

高超聲速流動領(lǐng)域的另一個難題是湍流，大渦模擬(LES)、直接數(shù)值模擬(DNS)、雷諾平均Navier-Stokes(RANS)以及RANS/LES混合方法等是目前數(shù)值模擬湍流的主要手段。其中，LES和DNS雖然在理論上較為完備，但極大的計算量難以滿足工程應(yīng)用需求，混合方法近些年雖有所發(fā)展，但要廣泛應(yīng)用于復(fù)雜外形計算仍有難度，而基于RANS的湍流模式理論雖然最初從不可壓流動發(fā)展而來，但目前仍是解決工程問題最有效、最實際的選擇。常用的湍流模型包括B-L代數(shù)模型、Spalart-Allmaras(S-A)一方程模型和剪切應(yīng)力輸運(yùn)(SST)兩方程模型以及多種可壓縮性修正方法[9]。

在湍流模擬中，高超聲速邊界層轉(zhuǎn)捩是其中非常重要但又難度很大且當(dāng)前最為關(guān)注的研究課題。航天飛機(jī)、高速導(dǎo)彈、臨近空間和再入飛行器等，因為全層流、全湍流和轉(zhuǎn)捩在摩阻和熱流上的差異很大，最大峰值會相差數(shù)倍。所以，高超聲速邊界層轉(zhuǎn)捩的建模與模擬研究是飛行器發(fā)展的迫切需求，美國2014—2030年的CFD技術(shù)路線圖將轉(zhuǎn)捩預(yù)測作為物理建模研究中的一項重要工作[4]。

1.2 邊界層轉(zhuǎn)捩

高超聲速邊界層轉(zhuǎn)捩受到來流雷諾數(shù)Re∞、馬赫數(shù)、湍流度、強(qiáng)激波、強(qiáng)逆壓梯度、粗糙度、橫流和第2模式不穩(wěn)定波等影響，準(zhǔn)確建模和預(yù)測十分困難。風(fēng)洞試驗和理論分析手段發(fā)展較早，為數(shù)值模擬建模研究提供了寶貴的參考數(shù)據(jù)和理論基礎(chǔ)，如美國蘭利中心Mach 6靜音風(fēng)洞[10]和波音/空軍Mach 6靜音風(fēng)洞[11]的一系列邊界層穩(wěn)定性試驗，Bi等[12]在炮風(fēng)洞中進(jìn)行的Mach 6圓錐轉(zhuǎn)捩試驗等。圖 6是Holden等在卡爾斯潘大學(xué)巴福羅研究中心(CUBRC)激波風(fēng)洞[13]中研究粗糙度和迎角等對HIFiRE錐體轉(zhuǎn)捩的影響，圖中給出了溫敏漆(TSP)熱流測量結(jié)果。

高超聲速邊界層轉(zhuǎn)捩數(shù)值模擬近年來發(fā)展較快得益于硬件的飛速發(fā)展，使DNS平板和尖錐等簡單外形高超聲速流動轉(zhuǎn)捩成為可能，圖 7為圓錐DNS計算表面熱流等值線和Q等值面[14]。但在工程應(yīng)用領(lǐng)域，RANS湍流/轉(zhuǎn)捩模式還是最為適用的。

圖6 HIFiRE外形轉(zhuǎn)捩[13]Fig.6 Transition on HIFiRE configuration[13]

圖7 表面熱流等值線和Q等值面[14]Fig.7 Surface contours of heat flux and Q isosurfaces[14]

通過類似渦黏性系數(shù)方式?；瘜恿髅}動影響，再結(jié)合湍流模式和間歇因子構(gòu)成完整轉(zhuǎn)捩過程模型，是當(dāng)今轉(zhuǎn)捩模型構(gòu)造的主流思想之一。如McDaniel[15-16]和Xiao[17]等利用該方法對Ma=3.5～8.0的平板、尖錐和裙錐等一系列高超聲速轉(zhuǎn)捩流動問題進(jìn)行數(shù)值模擬研究，其結(jié)果和實驗值接近，圖8是尖錐Bypass轉(zhuǎn)捩的熱流qw分布計算結(jié)果和實驗數(shù)據(jù)的比較，圖中BTU/(ft2·s)為英制熱流單位[17]。Papp和Dash[18-19]通過仔細(xì)調(diào)節(jié)間歇因子方程生成項的系數(shù)獲得了與風(fēng)洞試驗數(shù)據(jù)符合很好的轉(zhuǎn)捩起始位置和轉(zhuǎn)捩區(qū)長度，但所得模型的普適性較差。符松和王亮[20]建立了具有模擬高超聲速邊界層轉(zhuǎn)捩過程機(jī)制的k-ω-γ轉(zhuǎn)捩模型，圖9是高超聲速圓錐表面溫度分布，其中FSTI為自由來流湍流度。宋博和李椿萱[21]將Favré平均層流脈動動能方程和Navier-Stokes方程聯(lián)合構(gòu)造高超聲速轉(zhuǎn)捩模型，對Ma∞=5.91的裙錐的計算表明，該模型能較好地預(yù)測轉(zhuǎn)捩起始位置，但同樣存在參數(shù)選取經(jīng)驗性問題。

圖8 尖錐熱流分布[17]Fig.8 Heat flux distribution along rays of cone[17]

圖9 帶裙帶圓錐的表面溫度分布(Ma∞=5.91, Re∞=9.348×106,T∞=56.2 K, FSTI=0.1%)[20]Fig.9 Distribution of surface temperature on flared cone (Ma∞=5.91, Re∞=9.348×106,T∞=56.2 K, FSTI=0.1%)[20]

近年來，在邊界層轉(zhuǎn)捩模擬RANS方法研究領(lǐng)域，基于間歇因子和轉(zhuǎn)捩起點(diǎn)動量厚度雷諾數(shù)輸運(yùn)方程的γ-Reθ模型嶄露頭角，目前世界上主流的CFD軟件，如Fluent、CFD++、CFL3D和OVERFLOW等都有γ-Reθ模型計算模塊。γ-Reθ模型依據(jù)實驗數(shù)據(jù)和經(jīng)驗公式建立，完全基于局部變量，使得該模型具有很強(qiáng)的適用性和可擴(kuò)展性，在工程實踐中具有很大潛力。γ-Reθ模型最初由Menter等[22]提出，后來又與Langtry合作進(jìn)行改進(jìn)[23-25]。模型中兩個關(guān)鍵函數(shù)——轉(zhuǎn)捩長度函數(shù)和臨界動量厚度雷諾數(shù)函數(shù)，需要通過仔細(xì)的轉(zhuǎn)捩算例進(jìn)行校準(zhǔn)，作為商業(yè)秘密，Ansys CFX軟件直到2009年才公開這兩個函數(shù)[25]。這兩個函數(shù)對數(shù)值模擬平臺的依賴性很高，不同計算格式會對結(jié)果造成較大影響，Malan[26]、Medida[27]和Krause[28]等都在各自軟件平臺上進(jìn)行了專門的參數(shù)校準(zhǔn)。Bensassi等[29]模擬了軸對稱高超聲速圓錐轉(zhuǎn)捩流動，如圖10所示，圖中LCTM(Local Correlation-based Transition Model)為基于局部經(jīng)驗關(guān)系式的轉(zhuǎn)捩模型，Tu為湍流度。國內(nèi)，張曉東[30]和Yan等[31]也對γ-Reθ模型的應(yīng)用開展了研究。

圖10 圓錐轉(zhuǎn)捩(Ma=8)[29]Fig.10 Transition flow on cone (Ma=8)[29]

2 更精確——空間離散算法研究進(jìn)展

數(shù)值算法作為CFD的核心，是影響計算結(jié)果精度的重要因素之一，一直以來都是研究的重點(diǎn)，這其中又以Navier-Stokes控制方程組中對流項的格式構(gòu)造為主要研究對象。高超聲速流動中的激波強(qiáng)間斷和黏性干擾等現(xiàn)象要求格式在精細(xì)刻畫邊界層等黏性流動的同時，又能有效消除間斷附近的非物理振蕩，迎風(fēng)格式對拋物型問題具有天然的激波捕捉能力強(qiáng)、計算穩(wěn)定性好的優(yōu)點(diǎn)，在高超聲速數(shù)值模擬領(lǐng)域得到了廣泛的認(rèn)可和推廣[32]。

2.1 一階和二階迎風(fēng)格式

CFD的發(fā)展與計算機(jī)硬件水平的發(fā)展息息相關(guān)，正因如此，直到19世紀(jì)后期，CFD才真正具備求解實際問題的能力。對于高超聲速數(shù)值模擬來說，下面這些工作具有里程碑意義。

1959年，Godunov[33]提出一階迎風(fēng)格式，開創(chuàng)了基于Riemann精確解的CFD格式構(gòu)造方法；1979年，van Leer[34]提出單調(diào)迎風(fēng)守恒MUSCL方法，配合合適的限制器，這類格式具有很高的間斷捕捉分辨率，幾乎是目前高分辨率格式的通用做法。此后，Steger[35]和van Leer[36]等建立矢通量分裂格式(FVS)，該格式捕捉強(qiáng)間斷(如激波)能力強(qiáng)、可靠性高、格式簡單、計算量小，被廣泛用于高超聲速應(yīng)用領(lǐng)域；Roe[37]和Osher[38]等建立通量差分分裂格式(FDS)，該格式對剪切層等類似于接觸間斷的黏性作用區(qū)具有很高的分辨率，但需要引入熵修正消除紅玉現(xiàn)象。圖11給出了Roe格式雙錐問題計算結(jié)果，能清晰可辨流場中三叉點(diǎn)、交叉激波以及再附點(diǎn)等結(jié)構(gòu)，充分顯示了FDS格式在接觸間斷模擬中的高分辨率[39]特性；而對流迎風(fēng)分裂AUSM系列格式[40]則是結(jié)合以上兩類格式的優(yōu)缺點(diǎn)而構(gòu)造的混合通量分裂格式，包括AUSM+、AUSMPW、AUSMPW+等，兼有間斷分辨率高和計算效率高的優(yōu)點(diǎn)。圖12是AUSM+-up2格式計算結(jié)果[41]。

另外一種重要的格式是1983年Harten[42]提出的總變差減小差分格式(TVD)，由于其捕捉激波無波動且分辨率高，在高超聲速工程應(yīng)用上取得很大成功，后來又衍生出一系列的TVD格式。

圖11 雙錐計算流場結(jié)構(gòu)[39]Fig.11 Flowfield structure for double-cone case[39]

國內(nèi)對CFD格式的研究基本同步于國外先進(jìn)水平，如李松波[43]的耗散守恒格式、張涵信[44]提出無波動無自由參數(shù)耗散性的差分格式(NND)等。其中，NND格式在高超聲速數(shù)值模擬領(lǐng)域有廣泛應(yīng)用[6,45-46]，圖13給出了Ma=8.15時雙橢球馬赫數(shù)等值線及壁面極限流線計算結(jié)果，表明了NND格式具有良好的激波及分離區(qū)流動結(jié)構(gòu)分辨率。

圖12 AUSM+-up2格式計算火箭外形的壓力系數(shù)分布和流線[41]Fig.12 Pressure coefficient and streamlines of rocket configuration computed by AUSM+-up2 scheme[41]

圖13 NND格式計算所得流場Fig.13 Flowfield computed by NND scheme

2.2 高階精度格式

雖然二階精度格式在高超聲速工程應(yīng)用領(lǐng)域獲得了極大成功，然而，湍流、激波/邊界層干擾及剪切層等問題的精細(xì)模擬對格式的高分辨率提出了更高要求。高階精度格式在流場細(xì)節(jié)刻畫方面的魅力驅(qū)使著CFD研究人員在追求更高階精度格式的道路上矢志不渝。

1987年，Harten等[47]提出基本無波動ENO格式，放寬了TVD中強(qiáng)調(diào)消除波動的限制，并得到不斷發(fā)展[48]；1994年，Liu等[49]建立WENO格式，顯著減少流場中虛假波動，并獲得迅速發(fā)展，如Shu[50]和Lin[51]等建立一系列WENO格式，并在高超聲速流動的精細(xì)數(shù)值模擬中有大量應(yīng)用，如湍流、邊界層轉(zhuǎn)捩、聲波感受性和激波/邊界層干擾等問題，圖14為七階WENO高分辨率求解圓錐高超聲速邊界層轉(zhuǎn)捩問題[52]，圖中z為模型流向位置，us為流向瞬時速度，圖中給出了周向-22.5°～-52.5°位置的流向瞬時速度等值面分布。

圖14 流向瞬時速度等值面[52]Fig.14 Contours of instantaneous velocity stream component[52]

圖15 HWCNS計算所得流場Fig.15 Flowfield computed by HWCNS

3 更高效——時間推進(jìn)方法研究進(jìn)展

影響CFD求解效率的因素很多，如數(shù)值算法、軟件數(shù)據(jù)結(jié)構(gòu)以及計算機(jī)運(yùn)行和存儲速度等，其中，數(shù)值算法是提高CFD求解效率的基礎(chǔ)，時間推進(jìn)方法又是目前多數(shù)高超聲速CFD軟件采用的離散方程迭代求解方法。

隱式時間推進(jìn)方法不受CFL(Courant-Friedrichs-Lewy)穩(wěn)定性條件限制，計算效率高，在高超聲速數(shù)值模擬中有廣泛應(yīng)用[63]，隱式時間離散、離散方程線化和線性方程求解[64]是其重要的3個步驟。其中，隱式時間離散有基于Taylor展開的隱式后差格式和隱式(半隱式)Runge-Kutta(R-K)格式，前者格式簡單、易于實現(xiàn)、存儲量小，適用于定常問題求解[65]，后者分辨率高、色散特性好，但存儲量較大[66]。在離散方程線化方面，Rizzetta等[67]發(fā)展了交錯方向隱式(ADI)方法，Whitfield和Taylor[68]建立了數(shù)值求解雅可比矩陣方法，相比于近似雅可比矩陣，線化過程引入誤差小，具有良好的收斂性和魯棒性。線性方程組求解是影響隱式時間推進(jìn)方法整體計算效率的重要方面，在高超聲速數(shù)值模擬中，常用的有松弛迭代、近似因子分解和Krylov方法等。其中，LU-SGS(Lower-Upper Symmetric Gauss-Seidel)方法[69]存儲量小、計算效率高，在高超聲速數(shù)值模擬中應(yīng)用最為廣泛，具有較好的發(fā)展，如Chen等[70]針對有限體積方法建立的BLU-SGS(Block LU-SGS)能夠顯著提高激波問題模擬的收斂速度，Jameson和Caughey[71]結(jié)合LU-SGS和多種加速收斂技術(shù)建立了一種高效隱式時間推進(jìn)方法。

與此同時，隱式時間推進(jìn)方法推廣到高階精度空間離散格式的研究也取得了成功，代表性的工作有：Rizzetta等[67]將ADI方法推廣到高階緊致格式中，具有較高的計算效率和模擬精度；Pareschi和Russo[72]將半隱式R-K方法應(yīng)用于WENO格式；Gottlieb和Mullen[73]將LU-SGS和WENO格式結(jié)合，計算結(jié)果與實驗數(shù)據(jù)吻合較好；Zhong[74]采用半隱式三階R-K求解Navier-Stokes方程，精確模擬了二維Ma=5和Ma=15的拋物線前緣對聲波的感受性；Ekaterinaris[75]針對線性緊致格式，構(gòu)造了四階精度隱式方法；張毅鋒[76]將GMRES算法用于WCNS格式的計算中，顯著提高了計算效率，如圖16所示，采用以LU-SGS、Gauss-Seidel(G-S)點(diǎn)松弛(PR)、線松弛(LR)為預(yù)處理的GMRES方法收斂性明顯優(yōu)于前3種方法，并且作為雙時間步方法的子迭代在聲波感受性問題的DNS模擬中進(jìn)行了非定常計算，圖中DT為時間步長，LUSGS為上下對稱G-S方法。

圖16 最優(yōu)收斂曲線(Ma=10，球頭)[76]Fig.16 Optimal convergence histories (Ma=10, sphere)[76]

由此可知，隱式時間推進(jìn)方法不僅推動了CFD本身的發(fā)展，同時也推動了CFD在復(fù)雜流動工程問題中的應(yīng)用，然而在高階精度數(shù)值模擬方面，仍存在諸多問題需要解決，如誤差和收斂性問題[75]、復(fù)雜網(wǎng)格影響問題[77]以及多方程耦合求解問題[78]等。

4 數(shù)值模擬準(zhǔn)度與精度的正確評價

要求高超聲速數(shù)值模擬的結(jié)果更加精確和準(zhǔn)確，意味著計算結(jié)果的不確定度(Uncertainty)和誤差(Error)更小。其中最為重要的兩個問題是：搞清楚誤差和不確定度的定義和來源，建立誤差和不確定度的分析方法。

4.1 誤差和不確定度的定義和來源

美國航空航天學(xué)會(AIAA)[79]定義誤差為“在建模與模擬的過程中可以被認(rèn)識到的缺陷或問題，而這種缺陷或問題并非是由于認(rèn)知的不足而導(dǎo)致”，不確定度為“在建模過程中由于認(rèn)識不足而導(dǎo)致的潛在的功能與表現(xiàn)上的不足”。另外，美國機(jī)械工程師學(xué)會(ASME)和美國Sandia國家實驗室等都對CFD的誤差和不確定度進(jìn)行了規(guī)范化定義[80]，張涵信院士[81]提出將不確定度解讀為計算值或?qū)嶒炛蹬c真值準(zhǔn)確到前n位，從而給出不確定度表達(dá)式和真值估算的原則。

不考慮計算機(jī)舍入、程序編寫錯誤以及幾何外形差異所帶來的誤差，在高超聲速數(shù)值模擬中，物理建模、時空離散以及來流條件變化是給計算結(jié)果帶來誤差和不確定度的主要來源。其中，物理建模誤差來自建立控制方程和邊界等時的基本假設(shè)、近似和簡化。然而，高超聲速問題的物理模型非常復(fù)雜，使得傳統(tǒng)誤差理論難以適用，通常只能以不同物理模型或參數(shù)導(dǎo)致的數(shù)據(jù)散布度給出名義上的不確定度；時空離散誤差包括時間離散、空間離散以及初始/邊界條件離散帶來的誤差；來流條件變化引入的誤差和不確定度主要來自大氣環(huán)境參數(shù)的不確定性和數(shù)據(jù)庫插值過程帶來的誤差等。圖17是不同大氣密度(包括航天飛機(jī)飛行實測值)與美國1976標(biāo)準(zhǔn)大氣的密度比值(ρ/ρ76)，其不確定度在高度H為48～69 km的范圍內(nèi)約為5%，而在70～80 km的范圍內(nèi)最大達(dá)到16%[82]。

圖17 不同大氣標(biāo)準(zhǔn)的密度比較[82]Fig.17 Comparison of different atmosphere density data[82]

4.2 誤差和不確定度的分析方法

對于高超聲速數(shù)值模擬，在所有誤差和不確定度來源中，只有少部分能進(jìn)行誤差估計，如迭代誤差和空間離散誤差，而大部分只能進(jìn)行不確定度估計，研究方法可分為概率方法(考慮模型輸入與輸出概率結(jié)構(gòu))和非概率方法。

4.2.1 誤差估計方法

CFD計算通常只收斂到可接受的含有一定迭代誤差的結(jié)果。E?a和Hoekstra[83]給出迭代誤差估計的建模方法，假設(shè)數(shù)值解近似為指數(shù)收斂，構(gòu)造幾何級數(shù)誤差模型形式，利用不同迭代步計算數(shù)據(jù)確定待定系數(shù)，最終獲得誤差估計模型。圖18是針對類航天飛機(jī)外形開展的迭代誤差估計，選取適當(dāng)?shù)臄?shù)據(jù)點(diǎn)采集間隔步i，可獲得與真實誤差(Real Error)較為一致的迭代誤差估計值。

而CFD中的空間離散誤差與計算網(wǎng)格密不可分，按所需網(wǎng)格數(shù)量的不同分為兩大類：一類是基于單套網(wǎng)格的方法，如基于離散誤差傳播方程的方法[84-85]；另一類是基于多套網(wǎng)格的方法，基本思想是針對不同尺度下的網(wǎng)格解fk，利用擬合或外插等手段，獲得剔除空間離散誤差的“精確解”fexact，并給出空間離散誤差的估計值。第2類方法的關(guān)鍵是“精確解”的求取，包括最小二乘法、多項式方法、冪指數(shù)律方法以及三次樣條方法等[86]。對于高超聲速數(shù)值模擬而言，基于多項式泰勒展開的Richardson外插法最為常用[87-89]。該方法的基本思想是利用誤差行列展開式，結(jié)合多套網(wǎng)格計算結(jié)果，通過網(wǎng)格收斂性研究，獲得網(wǎng)格外插解，并以此為基礎(chǔ)研究空間離散誤差。該方法要求所有網(wǎng)格解單調(diào)收斂并落于收斂域內(nèi)，且激波等強(qiáng)間斷區(qū)域可能出現(xiàn)不同階誤差項抵消而影響對精度的判斷，但目前來說，在高超聲速數(shù)值模擬方面，該方法仍是預(yù)測數(shù)值空間離散誤差最魯棒的方法。圖19是利用該方法針對T2-97外形開展的離散誤差估計，由圖可知，數(shù)值解隨著網(wǎng)格細(xì)化(h變小)，具有很好的收斂性，x/L=1.139 8位置的Stanton數(shù)數(shù)值預(yù)測精度約為二階[88]。

圖18 迭代誤差估計Fig.18 Estimation of iterative error

圖19 離散誤差估計[88]Fig.19 Estimation of spatial error[88]

4.2.2 不確定度估計方法

(1)

敏感性導(dǎo)數(shù)的求解是該方法的關(guān)鍵，常用的有有限差分、復(fù)變方程、自動微分、敏感性方程以及離散伴隨方法等。圖20是針對類航天飛機(jī)外形開展的來流敏感性不確定度估計研究，其中，壁面溫度Tw=400 K左右時，Tw每變化10 K將引起軸向力系數(shù)CA約0.6%的不確定度。

圖20 壁面溫度變化引起的氣動力系數(shù)不確定度估計值Fig.20 Estimation of aerodynamic coefficient uncertainty induced by variety of wall temperature

5 結(jié) 論

高超聲速流動因其固有的特點(diǎn)(強(qiáng)間斷、強(qiáng)黏性、真實氣體效應(yīng)、稀薄氣體效應(yīng)等)，對相應(yīng)的數(shù)值模擬技術(shù)提出了更高的要求。雖然近年來在網(wǎng)格生成、數(shù)學(xué)物理建模、數(shù)值方法構(gòu)造以及結(jié)果可靠性評估等方面取得了巨大成就，但為了獲得更準(zhǔn)確、更可靠的氣動數(shù)據(jù)，更為精細(xì)地研究高超聲速復(fù)雜物理現(xiàn)象，仍有大量工作要做。

1) 在物理模型方面，目前所發(fā)展的化學(xué)反應(yīng)模型基本能滿足高溫真實氣體效應(yīng)的研究，但對于更高飛行速度范圍，如第二宇宙速度再入問題，目前的物理模型是否適合仍需確認(rèn)。按流域劃分稀薄氣體效應(yīng)研究方法能滿足跨流域計算的需要，但仍需開展具有更高效率的跨流域統(tǒng)一算法研究?；赗ANS的湍流模式理論能解決復(fù)雜布局的流動問題模擬，但這些模型的驗證與確認(rèn)仍是最大的挑戰(zhàn)。針對簡單問題的高超聲速邊界層轉(zhuǎn)捩研究取得了一些進(jìn)展，但對于復(fù)雜流動或復(fù)雜外形仍有待進(jìn)一步研究，尤其需重點(diǎn)關(guān)注γ-Reθ模型能否很好地模擬高超聲速邊界層轉(zhuǎn)捩流動問題。

2) 在數(shù)值計算方法方面，目前常規(guī)的二階迎風(fēng)格式在流動模擬中取得了很大成就，基本能滿足飛行器工程設(shè)計及流動機(jī)理分析的需要，高階格式應(yīng)用也獲得了初步成功，但對于更為精細(xì)的模擬，如臨近空間環(huán)境下飛行器摩擦阻力與表面熱流預(yù)測、精細(xì)的復(fù)雜流動現(xiàn)象研究(如激波/邊界層干擾、分離區(qū)模擬等)，仍需要大力發(fā)展在色散、耗散特性、邊界層分辨率以及分離區(qū)刻畫等方面具有優(yōu)良特性的高階精度格式，尤其是解決這些格式在復(fù)雜問題中的應(yīng)用問題。在提高計算效率方法，不僅要發(fā)展具有更高效率的求解方法，更要與高性能計算機(jī)(HPC)相結(jié)合，充分利用HPC硬件與其先進(jìn)的框架，提高CFD軟件的大規(guī)模并行能力。

3) 在數(shù)值模擬可信度研究方面，針對高超聲速數(shù)值模擬的特點(diǎn)，理清了迭代誤差估計、空間離散誤差估計以及不確定度的非概率及概率方法的原理、流程，初步開展了E?a迭代誤差建模方法、Richardson外插法以及敏感性導(dǎo)數(shù)不確定度估計等方法在高超聲速數(shù)值模擬可信度研究中的適用性驗證。下一步需開展高超聲速CFD模擬的驗證與確認(rèn)研究，關(guān)注確認(rèn)數(shù)據(jù)庫的建設(shè)，重點(diǎn)進(jìn)行誤差與不確定度量化分析。

此外，在復(fù)雜外形幾何處理與網(wǎng)格生成、知識提取(海量數(shù)據(jù)處理、可視化、數(shù)據(jù)庫生成等)等方面仍需開展大量工作?？傮w上，在高超聲速流動數(shù)值模擬方法方面，國內(nèi)研究工作已在國際上占有一席之地，在湍流模擬研究中也有一些可圈可點(diǎn)的工作，但在數(shù)值模擬結(jié)果可靠性評估方法方面，還處于起步階段。

[1] Guo J.Military application value of near space vehicle[J].OME Information, 2010, 27(8): 22-27 (in Chinese). 郭勁. 臨近空間飛行器軍事應(yīng)用價值分析[J].光機(jī)電信息, 2010, 27(8): 22-27.

[2] Knight D.RTO WG 10: test cases for CFD validation of hypersonic flight, AIAA-2002-0433[R].Reston: AIAA, 2002.

[3] Saric W S, Muylaert J, Dujarric C.Hypersonic experimental and computational capability, improvement and validation, AGARD-AR-319[R].Hull (Québec): North Atlantic Treaty Organisation, 1996.

[4] Slotnick J, Khodadoust A, Alonso J.CFD vision 2030 study: a path to revolutionary computational aerosciences, NASA/CR-2014-218178[R].Washington, D.C.: NASA, 2014.

[5] Sarma G S R.Physico-chemical modeling in hypersonic flow simulation[J].Progress in Aerospace Sciences, 2000, 36(3-4): 281-349.

[6] Li Z H, Zhang H X.Study on gas kinetic unified algorithm for flows from rarefied transition to continuum[J].Acta Aerodynamica Sinica, 2003, 21(3): 255-266 (in Chinese). 李志輝, 張涵信. 稀薄流到連續(xù)流的氣體運(yùn)動論統(tǒng)一算法研究[J]. 空氣動力學(xué)學(xué)報, 2003, 21(3): 255-266.

[7] Park C.Nonequilibrium hypersonic aerothermodynamics[M]. New York: John Wiley & Sons Inc., 1989.

[8] Prabhu D K, Papadopoulos P E, Davies C B, et al. Shuttle orbiter contingency abort aerodynamics: real-gas effects and high angles of attack, RTO-EN-AVT-116[R].Rhode St.Genèse, Belgium: von Kármán Institute, 2004.

[9] Smits A J, Martin M P, Girimaji S. Current status of basic research in hypersonic turbulence, AIAA-2009-0151[R]. Reston: AIAA, 2009.

[10] Wilkinson S P.A review of hypersonic boundary layer stability experiments in a quiet Mach 6 wind tunnel, AIAA-1997-1819[R]. Reston: AIAA, 1997.

[11] Abney A, Ward C, Berridge D, et al.Hypersonic boundary-layer transition experiments in the Boeing/AFOSR Mach-6 quiet tunnel, AIAA-2013-0375[R].Reston: AIAA, 2013.

[12] Bi Z X, Zhu N J, Shen Q, et al.Measurements on transition of hypersonic boundary layer over circular cone, AIAA-2007-6728[R].Reston: AIAA, 2007.

[13] Holden M S, Wadhams T P, MacLean M, et al.Review of studies of boundary layer transition in hypersonic flows over axisymmetric and elliptic cones conducted in the CUBRC shock tunnels, AIAA-2009-0782[R]. Reston: AIAA, 2009.

[14] Gronvall J E, Johnson H B, Candler G V.Hypersonic three-dimensional boundary layer transition on a cone at angle of attack, AIAA-2012-2822[R]. Reston: AIAA, 2012.

[15] McDaniel R D, Nantie R P, Hassan H A.Transition onset prediction for high speed flow, AIAA-1999-3792[R].Reston: AIAA, 1999.

[16] McDaniel R D, Hassan H A.Role of bypass transition in conventional hypersonic facilities, AIAA-2001-0209[R].Reston: AIAA, 2001.

[17] Xiao X, Edwards J R, Hassan H A.Transition flow over an elliptic cone at Mach 8, AIAA-2001-0276[R]. Reston: AIAA, 2001.

[18] Papp J L, Dash S M. Extensions of a rapid engineering approach to modeling hypersonic laminar-to-turbulent transitional flows, AIAA-2005-0892[R].Reston: AIAA, 2005.

[19] Papp J L， Dash S M.Modeling hypersonic laminar to turbulent transitional flows for 3D geometries using two-equation onset and intermittency transport models, AIAA-2012-0449[R].Reston: AIAA, 2012.

[20] Fu S, Wang L.Numerical simulation of hypersonic boundary layer transition using RANS[J].Science China: G, 2009, 39(4): 617-626 (in Chinese). 符松, 王亮. 基于雷諾平均方法的高超音速邊界層轉(zhuǎn)捩模擬[J].中國科學(xué): G輯, 2009, 39(4): 617-626.

[21] Song B, Li C H.A Favré averaged transition prediction model for hypersonic flows[J].Science China Technology Science, 2010, 40(8): 879-885 (in Chinese). 宋博, 李椿萱.基于Favré平均的高超聲速可壓縮轉(zhuǎn)捩預(yù)測模型[J].中國科學(xué): 技術(shù)科學(xué), 2010, 40(8): 879-885.

[22] Menter F R, Esch T, Kubacki S.Transition modelling based on local variables[C]∥Proceedings of the 5th International Symposium on Engineering Turbulence Modelling and Measurements.Amsterdam: Elsevier, 2002: 555-564.[23] Menter F R, Langtry R B, Likki S R, et al.A correlation-based transition model using local variables, Part I: model formulation[J].Journal of Turbomachinery, 2004, 128(3): 413-422.

[24] Langtry R B.A correlation-based transition model using local variables for unstructured parallelized CFD codes[D]. Stuttgart: Institut für Thermische Str?mungsmachinen und Maschinenlaboratorium, 2006.

[25] Langtry R B, Menter F R.Correlation-based transition modeling for unstructured parallelized computational fluid dynamic codes[J].AIAA Journal, 2009, 47(12): 2894-2906.

[26] Malan P, Suluksna K, Juntasaro E.Calibrating theγ-Reθtransition model for commercial CFD, AIAA-2009-1142 [R].Reston: AIAA, 2009.

[27] Medida S, Baeder J D.Application of the correlation-basedγ-Reθ ttransition model to the Spalart-Allmaras turbulence model, AIAA-2011-3979[R].Reston: AIAA, 2011.

[28] Krause M, Behr M, Ballmann J.Modeling of transition effects in hypersonic intake flows using a correlation-based intermittency model, AIAA-2008-2598[R].Reston: AIAA, 2008.

[29] Bensassi K, Lani A, Rambaud P.Numerical investigations of local correlation-based transition model in hypersonic flows, AIAA-2012-3151[R].Reston: AIAA, 2012.

[30] Zhang X D, Gao Z H.Numerical discuss to complete empirical correlation in Langtry’s transition model[J].Applied Mathematics and Mechanics, 2010, 31(5): 544-552(in Chinese). 張曉東, 高正紅.關(guān)于補(bǔ)充Langtry的轉(zhuǎn)捩模型經(jīng)驗修正式的數(shù)值探討[J].應(yīng)用數(shù)學(xué)和力學(xué), 2010, 31(5): 544-552.

[31] Yan P, Mao R, Qiang X.Application of PSE analysis method with transitional turbulence model on hypersonic flows, AIAA-2011-3981[R].Reston: AIAA, 2011.

[32] Yan C, Zhang Z, Zhang L X, et al.Characteristic analysis of the upwind scheme[J].Acta Aerodynamica Sinica, 2003, 21(3): 336-341 (in Chinese). 閻超, 張智, 張立新, 等.上風(fēng)格式的若干性能分析[J].空氣動力學(xué)學(xué)報, 2003, 21(3): 336-341.

[33] Godunov S K.A difference method for the numerical calculation of discontinuous solutions of hydrodynamic equations[J].Matematicheskii Sbornik, 1959, 47(89): 271-306.

[34] van Leer B.Towards the ultimate conservation difference scheme V: a second-order sequel to Godunov’s method[J]. Journal of Computational Physics, 1979, 32(1): 101-136.

[35] Steger J L, Warming R F.Flux vector splitting of the inviscid gasdynamics equations with application to finite difference methods[J].Journal of Computational Physics, 1981, 40(2): 263-293.

[36] van Leer B.Flux vector splitting for Euler equations[C]∥Eighth International Conference on Numerical Methods in Fluid Dynamics.Berlin: Springer Heidelberg, 1982: 507-512.

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

[38] Osher S, Solomon F.Upwind difference schemes for hyperbolic conservation laws[J].Mathematics of Computation, 1982, 38(158): 339-374.

[39] Knight D, Longo J, Drikakis D, et al.Assessment of CFD capability for prediction of hypersonic shock interactions[J].Progress in Aerospace Sciences, 2012, 48-49: 8-26.

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

[41] Kitamura K, Shima E.A new pressure flux for AUSM-family schemes for hypersonic heating computations, AIAA-2011-3056[R].Reston: AIAA, 2011.

[42] Harten A.High resolution schemes for hypersonic conservation laws[J].Journal of Computational Physics, 1983, 49(3): 357-393.

[43] Li S B.Dissipation conservative scheme theory[M].Beijing: Higher Education Press, 1997: 1 (in Chinese). 李松波.耗散守恒格式理論[M].北京: 高等教育出版社, 1997: 1.

[44] Zhang H X.None-oscillation, none free coefficients dissipative finite difference scheme[J].Acta Aerodynamica Sinica, 1988, 6(2): 143-165 (in Chinese). 張涵信.無波動、無自由參數(shù)的耗散差分格式[J].空氣動力學(xué)學(xué)報, 1988, 6(2): 143-165.

[45] Chen J Q, Jiang D W, Zhang Y F.The study on the precision of numerical simulation for lateral jets flow and the experiment validation[J].Acta Aerodynamica Sinica, 2010, 28(4): 421-425 (in Chinese). 陳堅強(qiáng), 江定武, 張毅鋒.側(cè)向噴流數(shù)值模擬精度及實驗驗證研究[J].空氣動力學(xué)學(xué)報, 2010, 28(4): 421-425.

[46] Mao M L, Jiang D W, Deng X G.Study of hybrid scheme for the prediction of aerodynamic heat transfer-rate in hypersonic laminar flow[J].Acta Aerodynamica Sinica, 2009, 27(3): 275-280 (in Chinese). 毛枚良, 江定武, 鄧小剛.高超聲速層流氣動熱預(yù)測混合算法研究[J].空氣動力學(xué)學(xué)報, 2009, 27(3): 275-280.

[47] Harten A, Engquist B, Osher S, et al.Uniformly high-order accurate essentially non-oscillatory schemes III[J].Journal of Computational Physics, 1987, 71(2): 231-303.

[48] Stiriba Y.A nonlinear flux split method for hyperbolic conservation laws[J].Journal of Computational Physics, 2002, 176(1): 20-39.

[49] Liu X D, Osher S, Chan T.Weighted essentially non-oscillatory schemes[J].Journal of Computational Physics, 1994, 115(1): 200-212.

[50] Shu C W.Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws, NASA/CR-97-206253[R].Washington, D.C.: NASA, 1997.

[51] Lin S Y, Hu J J.Parametric study of weighted essentially nonoscillatory schemes for computational aeroacoustics[J]. AIAA Journal, 2001, 39(3): 371-379.

[52] Li X L, Fu D X, Ma Y W.Direct numerical simulation of hypersonic boundary layer transition over a blunt cone with a small angle of attack[J].Physics of Fluids, 2010, 22(2): 025105-025105-18.

[53] Fu D X.Fluid dynamics numerical simulation[M].Beijing: National Defense Industry Press, 1993: 235-249 (in Chinese). 傅德薰.流體力學(xué)數(shù)值模擬[M].北京: 國防工業(yè)出版社, 1993: 235-249.

[54] He G H.Third-order ENN scheme and its application to the calculation of complex hypersonic viscous flows[D].Mianyang: China Aerodynamics Research and Development Center, 1994 (in Chinese). 賀國宏. 三階ENN格式及其在高超聲速黏性復(fù)雜流場求解中的應(yīng)用[D].綿陽：中國空氣動力研究與發(fā)展中心, 1994.

[55] Li Q, Zhang H X, Gao S C.Numerical simulations on supersonic shear layer flow[J].Acta Aerodynamica Sinica, 2000, 18(1): 67-77 (in Chinese). 李沁, 張涵信, 高樹椿.關(guān)于超聲速剪切流動的數(shù)值模擬[J]. 空氣動力學(xué)學(xué)報, 2000, 18(1): 67-77.

[56] Deng X G, Zhang H X.Developing high-order accurate nonlinear schemes[J].Journal of Computational Physics, 2000, 165(1): 22-44.

[57] Deng X G, Mao M L, Liu J C.High-order dissipative weighted compact nonlinear schemes for Euler and Navier-Stokes equations, AIAA-2001-2626[R].Reston: AIAA, 2001.

[58] Ren Y X, Liu M, Zhang H X.A characteristic-wise hybrid compact-WENO scheme for solving hyperbolic conservation laws[J].Journal of Computational Physics, 2003, 192(2): 365-386.

[59] Ren Y X.A robust shock-capturing scheme based on rotated Riemann solvers[J].Computers & Fluids, 2003, 32(10): 1379-1403.

[60] Zheng H S, Zhao N, Cheng J.One dimensional high order accurate discontinuous GDQ methods[J].Mathematica Numerica Sinica, 2004, 26(3): 293-302 (in Chinese). 鄭華盛, 趙寧, 成娟.一維高精度離散GDQ方法[J].計算數(shù)學(xué), 2004, 26(3): 293-302.

[61] Zheng H S, Zhao N.A high order accurate TVD difference scheme for hyperbolic conservation laws[J].Chinese Journal of Computational Physics, 2005, 22(1): 13-18 (in Chinese). 鄭華盛, 趙寧.雙曲型守恒律的一種高精度TVD差分格式[J].計算物理, 2005, 22(1): 13-18.

[62] Zheng H S, Zhao N.A high order accurate MmB difference scheme for nonlinear hyperbolic conservation laws[J]. Chinese Journal of Computational Mechanics, 2006, 23(2): 218-222 (in Chinese). 鄭華盛, 趙寧.非線性雙曲型守恒律的高精度MmB差分格式[J].計算力學(xué)學(xué)報, 2006, 23(2): 218-222.

[63] Walpot L.Numerical analysis of the ARD capsule in S4 wind tunnel[C]∥4th European Symposium on Aerothermodynamics for Space Vehicles, Capua.Italy: European Space Agency, 2001.

[64] Dennis J E, Schnabel R B.Numerical methods for unconstrained optimization and nonlinear equations[M].Englewood Cliffs, NJ: Prentice-Hall, 1983: 2-14.

[65] Tysinger T L, Caughey D A.Implicit multigrid algorithm for the Navier-Stokes equations[C]∥29th Aerospace Sciences Meeting.Reno, NV: AIAA, 1991.

[66] Hixon R, Visbal M R.Comparison of high-order implicit time marching schemes for unsteady flow calculations, AIAA-2007-4324[R].Reston: AIAA, 2007.

[67] Rizzetta D P, Visbal M R, Blaisdell G A.A time-implicit high-order compact differencing and filtering scheme for large-eddy simulation[J].International Journal for Numerical Methods in Fluids, 2003, 42(6): 666-693.

[68] Whitfield D L, Taylor L K.Discretized Newton-relaxation solution of high resolution flux-difference split schemes, AIAA-1991-1539[R]. Peston: AIAA, 1991.

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

[70] Chen R F, Wang Z J.Fast, block lower-upper symmetric Gauss-Seidel scheme for Arbitrary grids[J].AIAA Journal, 2000, 38(12): 2238-2245.

[71] Jameson A, Caughey D A.How many steps arerequired to solve the Euler equations of steady, compressibleflow: in search of a fast solution algorithm, AIAA-2001-2673[R].Reston: AIAA, 2001.

[72] Pareschi L, Russo G.Implicit-explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxation[J].Journal of Scientific Computing, 2005, 25(1): 129-155.

[73] Gottlieb S, Mullen J S.An implicit WENO scheme for steady-state computation of scalar hyperbolic equations[J]. Computational Fluid and Solid Mechanics, 2003, 1(2): 1946-1950.

[74] Zhong X. Direct numerical simulation of hypersonic boundary-layer transition over blunt leading edges, Part I: a new numerical method and validation, AIAA-1997-0755[R]. Reston: AIAA, 1997.

[75] Ekaterinaris J A.Implicit, high-resolution, compact schemes for gas dynamics and aeroacoustics[J].Journal of Computational Physics, 1999, 156(2): 272-299.

[76] Zhang Y F. Investigations of convergence acceleration and complex flow numerical simulation for high-order accurate scheme(WCNS)[D]. Mianyang: China Aerodynamics Research and Development Center, 2007 (in Chinese). 張毅鋒.高精度格式(WCNS)加速收斂和復(fù)雜流動數(shù)值模擬的應(yīng)用研究[D]. 綿陽: 中國空氣動力研究與發(fā)展中心, 2007.

[77] Buelow P E O, Venkateswaran S, Merkle C L.Stability and convergence analysis of implicit upwind schemes[J].Computers & Fluids, 2001, 30(7-8): 961-988.

[78] Chisholm T T, Zingg D W.A Jacobian-free Newton-Krylov algorithm for compressible turbulent flows[J]. Journal of Computational Physics, 2009, 228(9): 3490-3507.

[79] AIAA.Guide for the verification and validation of computational fluid dynamics simulations, AIAA-G-077-1998[R].Reston: AIAA, 1998.

[80] Oberkampf W L, Roy C J.Verification and validation in scientific computing[M].New York: Cambridge University Press, 2010: 19-143.

[81] Zhang H X.On the uncertainty about CFD result[J].Acta Aerodynamica Sinica, 2008, 26(1): 47-49 (in Chinese). 張涵信.關(guān)于CFD計算結(jié)果的不確定度問題[J].空氣動力學(xué)學(xué)報, 2008, 26(1): 47-49.

[82] Champion K S W.Middle atmosphere density data and comparison with models[J].Advances in Space Research, 1990, 10(6): 17-26.

[83] E?a L, Hoekstra M.Code verification and verification of calculations with RANS solvers, RTO-MP-AVT-147-P-07[R].New Jersey: North Atlantic Treaty Orgnisation, 2010.

[84] Zhang X D, Pelletier D, Trépanier J Y, et al.Numerical assessment of error estimators for Euler equations[J].AIAA Journal, 2001, 39(9): 1706-1715.

[85] Qin Y.A discrete transport equation for error estimation in CFD, AIAA-2002-0906R].Reston: AIAA, 2002.

[86] Celik I, Li J.Assessment of numerical uncertainty for the calculations of turbulent flow over a backward facing step[J].International Journal for Numerical Methods in Fluids, 2005, 49(9): 1015-1031.

[87] Roy C J.Grid convergence error analysis for mixed-order numerical schemes, AIAA-2001-2606[R].Reston: AIAA, 2001.

[88] Chen J Q, Zhang Y R.Verification and validation in CFD based on the Richardson extrapolation method[J].Acta Aerodynamica Sinica, 2012, 30(2): 176-183 (in Chinese). 陳堅強(qiáng), 張益榮.基于Richardson插值法的CFD驗證和確認(rèn)方法的研究[J].空氣動力學(xué)學(xué)報, 2012, 30(2): 176-183.

[89] Oberkampf W L, Trucano T G.Validation methodology in computational fluid dynamics, AIAA-2002-2549[R].Reston: AIAA, 2002.

[90] Perez R A.Uncertainty analysis of computational fluid dynamics via polynomial chaos[D].Blacksburg: Viginia Polytechnic Institute and State University, 2008.

[91] Walters R W, Huyse L.Uncertainty analysis for fluid mechaics with applications, NASA/CR-2002-211449[R].Washington, D.C.: NASA, 2002.

Tel: 0816-2472011

E-mail: tangzhigong@sina.com

陳堅強(qiáng) 男， 博士， 研究員， 博士生導(dǎo)師。主要研究方向：高超聲速數(shù)值模擬研究。

Tel: 0816-2463011

E-mail: jq-chen@263.net

*Corresponding author. Tel.： 0816-2463011 E-mail: jq-chen@263.net

More fidelity, more accurate, more efficient—progress on numerical simulations for hypersonic flow

TANG Zhigong, ZHANG Yirong, CHEN Jianqiang*, MAO Meiliang, ZHANG Yifeng, LIU Huayong

ChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China

The progress of numerical simulations capability, which is related to fidelity, accuracy and efficiency, for prediction of hypersonic flow during the last several decades are reviewed. For physical modeling, the current status of modeling and simulation about high-temperature gas effects, rarefied gas effects and turbulence effects in hypersonic flow simulations is reviewed. And great emphasis is placed on hypersonic boundary layer transition modeling based on Reynolds-averaged Navier-Stokes (RANS) equations. For spatial discretization schemes, a review of development and application for usual 2nd order upwind schemes and high-order schemes are given. For time marching schemes, development and application for implicit time marching schemes is briefly summarized. For uncertainty quantification, a review of its concept, source and analysis method is firstly addressed, and then the primitive results of iterative error estimation, spatial discretization error analyses using Richardson extrapolation and uncertainty quantification using sensitivity derivatives are also presented. At last, some possible future directions on numerical simulations for hypersonic flow are discussed.

hypersonic flow; numerical simulation; physical modeling; numerical method; error; uncertainty

2014-07-25； Revised： 2014-08-26； Accepted： 2014-09-15； Published online： 2014-09-16 16:02

National Nature Science Foundation of China (11372342)

2014-07-25; 退修日期： 2014-08-26; 錄用日期： 2014-09-15; 網(wǎng)絡(luò)出版時間： 2014-09-16 16：02

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

國家自然科學(xué)基金 (11372342)

Tang Z G, Zhang Y R, Chen J Q, et al. More fidelity, more accurate, more efficient—progress on numerical simulations for hypersonic flow[J]. Acta Aeronautica et Astronautica Sinica， 2015， 36(1)： 120-134. 唐志共， 張益榮， 陳堅強(qiáng)， 等. 更準(zhǔn)確、更精確、更高效——高超聲速流動數(shù)值模擬研究進(jìn)展[J]. 航空學(xué)報， 2015， 36(1)： 120-134.

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

10.7527/S1000-6893.2014.0231

V211.3

A

1000-6893(2015)01-0120-15

唐志共 男， 博士， 研究員， 博士生導(dǎo)師。主要研究方向：超聲速與高超聲速空氣動力學(xué)研究。

*通訊作者.Tel.： 0816-2463011 E-mail: jq-chen@263.net

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