高超聲速湍流直接數(shù)值模擬技術(shù)

李新亮

中國科學(xué)院力學(xué)研究所 高溫氣體動力學(xué)國家重點實驗室， 北京 100190

高超聲速湍流直接數(shù)值模擬技術(shù)

李新亮*

中國科學(xué)院力學(xué)研究所 高溫氣體動力學(xué)國家重點實驗室， 北京 100190

概述了近年來國內(nèi)外高超聲速湍流直接數(shù)值模擬(DNS)技術(shù)方面最新的研究進展，主要集中在高精度、高魯棒性數(shù)值方法方面，同時也介紹了近年來典型的高超聲速湍流DNS算例。在數(shù)值方法方面，主要介紹了高精度激波捕捉格式以及保持計算穩(wěn)定的數(shù)值技術(shù)，重點是WENO格式及高階保單調(diào)格式的最新進展。在高超聲速湍流DNS算例方面，介紹了壓縮性影響、壁溫影響、真實氣體效應(yīng)以及高超聲速轉(zhuǎn)捩等方面的DNS研究。此外，還簡要介紹了作者開發(fā)的可壓縮高精度計算流體力學(xué)軟件OpenCFD。

高超聲速； 湍流； 直接數(shù)值模擬； 高精度數(shù)值方法； OpenCFD

受航空航天等領(lǐng)域的需求牽引，高超聲速湍流(通常指來流馬赫數(shù)超過5的湍流)成為目前流體力學(xué)領(lǐng)域的熱點。高超聲速飛行器的內(nèi)、外流動以湍流及轉(zhuǎn)捩為主要特征。層流轉(zhuǎn)捩到湍流后，摩阻及熱流增加為原來的數(shù)倍。因而，無論氣動力還是熱防護設(shè)計，都需要對轉(zhuǎn)捩與湍流進行精細預(yù)測。此外，燃料混合、氣動光學(xué)以及氣動聲學(xué)都與湍流及轉(zhuǎn)捩密切相關(guān)，這就對湍流及轉(zhuǎn)捩的機理、模型及控制提出了嚴格的要求。直接數(shù)值模擬(DNS)是研究湍流與轉(zhuǎn)捩重要的數(shù)值計算手段。該方法不采用任何湍流及轉(zhuǎn)捩模型，而是使用足夠密集的網(wǎng)格直接離散求解控制流動的Navier-Stokes方程，從而得到湍流及轉(zhuǎn)捩時空演化的全部細節(jié)。DNS無湍流模型誤差，是目前最為準確可靠的湍流計算手段，也是校驗及改進湍流及轉(zhuǎn)捩模型的有效方法。同時，DNS可以給出湍流全部尺度的信息，非常有利于進行湍流和轉(zhuǎn)捩的機理及控制研究。DNS的不足之處是計算量非常巨大，目前尚無法直接應(yīng)用于飛行器全機尺度的工程問題計算。但這并不妨礙DNS作為研究湍流機理、改進湍流模型以及探索湍流控制的有力工具。隨著計算機性能的快速發(fā)展，DNS將在湍流研究中發(fā)揮越來越重要的作用[1]。

相對于亞、跨及超聲速湍流，高超聲速湍流DNS開展得比較晚，直到近幾年才逐漸開展起來。這主要是受制于計算方法及計算量的限制。湍流DNS對數(shù)值方法要求苛刻，要求數(shù)值方法必須是低耗散的，否則就會將湍流小尺度耗散掉，從而影響小尺度的分辨率。為了捕捉激波，必須用足夠的耗散抑制數(shù)值振蕩。低耗散與保持計算穩(wěn)定這一對矛盾很難平衡，尤其是對于高超聲速湍流[2]。此外，由于高超聲速湍流通常具有很高的雷諾數(shù)，其DNS計算量非常大，也給計算造成了一定的困難。近幾年來，隨著高精度激波捕捉方法的快速發(fā)展以及計算機性能的提升，高超聲速湍流DNS逐漸發(fā)展起來。

1 面向高超聲速湍流的高精度方法

高超聲速湍流的高分辨率計算(DNS及大渦模擬)對數(shù)值方法要求非?？量獭Ｓ绕涫荄NS，要求數(shù)值方法必須是高分辨及低耗散的，而且必須有非常強的激波捕捉能力及計算穩(wěn)定性。本節(jié)主要介紹面向高超聲速湍流高分辨率計算的高精度激波捕捉方法以及必要的保持計算穩(wěn)定性技術(shù)。

1.1 高精度激波捕捉格式

由于DNS需要分辨湍流的全部流動細節(jié)，因而要求數(shù)值方法具有高分辨率、低耗散及低色散等特征。譜方法具有理論上的無窮階精度，對湍流DNS非常有利，但譜方法無法處理含間斷的流場(否則會出現(xiàn)非物理振蕩)，而對于高超聲速湍流，激波、接觸間斷等間斷隨處可見，這就限制了譜方法的應(yīng)用。差分方法易于實現(xiàn)高精度、高分辨率及激波捕捉，易于實現(xiàn)湍流DNS。與有限體積方法相比，差分法不易處理復(fù)雜網(wǎng)格，因而復(fù)雜工程計算中更傾向于有限體積方法。但由于目前湍流DNS多采用較為簡單的幾何外形，因而采用高精度差分方法更為便捷。

差分方法有中心差分及迎風(fēng)差分2種類型。中心差分具有更低的耗散，但數(shù)值穩(wěn)定性較差，需要添加人工黏性或數(shù)值濾波方法保持計算穩(wěn)定。為了克服非線性效應(yīng)產(chǎn)生的混淆誤差，使用中心差分時需要對對流項進行特殊處理[2]，另外，在有激波的情況下，中心差分格式需要與迎風(fēng)型的激波捕捉格式(如WENO(Weighted Essentially Non-Oscillatory)格式等)混合使用，在激波區(qū)切換到激波捕捉格式。相對于中心格式，迎風(fēng)格式(通常為激波捕捉格式)耗散更大，但魯棒性更強。對于像高超聲速湍流這樣復(fù)雜的流場，更適合采用迎風(fēng)型格式。目前常用的迎風(fēng)型高階激波捕捉格式有WENO格式[3-4]及保單調(diào) (MP) 格式[5]等。此外，中國學(xué)者開發(fā)的加權(quán)緊致非線性(WCNS)[6]以及群速度控制 (GVC) 格式[7-8]等也是較好的選擇。

1.1.1 WENO格式

Liu等[3]對高精度的ENO (Essentially Non-Oscillatory) 格式進行改進，利用加權(quán)平均的思想提高了方法的計算效率，構(gòu)造了WENO格式。隨后，Jiang 和Shu[4]對該方法進行了進一步改進，構(gòu)造了新的光滑度量因子，進一步提高了計算效率。Jiang 和Shu改進的WENO格式[4]是目前應(yīng)用最為廣泛的WENO格式 (以下簡稱為WENO-JS格式). 近年來，人們對WENO-JS格式進行了改進和優(yōu)化，推出了一系列改進的WENO格式。WENO格式的最大特點是保持高精度的同時還具有非常好的魯棒性，因而在目前可壓縮尤其是超、高超聲速復(fù)雜流動中得到了較為廣泛的應(yīng)用。WENO格式的不足之處在于其數(shù)值耗散仍然有些偏大，尤其是對于湍流小尺度耗散過大[1]，此外相對于線性格式，WENO格式的計算量仍然有些大。

針對WENO-JS格式的不足，近年來發(fā)展出了各種改進版本。Borges等[9]對WENO-JS格式的光滑度量因子進行了改進，提出了WENO-Z格式。該格式僅需在WENO-JS基礎(chǔ)上改動幾行代碼即可實現(xiàn)，非常便于WENO-JS格式程序的移植。與WENO-JS格式相比，WENO-Z格式在極值點(包括高階極值點)附近的計算精度得到了改善。作者也對WENO-Z格式和WENO-JS格式進行了對比測試，發(fā)現(xiàn)對于較為簡單的問題(如一、二維的模型問題)WENO-Z格式有較好的表現(xiàn)，但對于三維復(fù)雜問題(如湍流DNS),WENO-Z格式與WENO-JS格式計算結(jié)果相近，而魯棒性測試則顯示W(wǎng)ENO-JS格式的魯棒性略優(yōu)于WENO-Z格式。

Martin等[10]針對WENO格式進行了系數(shù)優(yōu)化，提出了分辨率優(yōu)化的WENO(WENO-SYMBO/SYMOO)格式。與經(jīng)典的WENO-JS格式不同，該格式基于對稱型網(wǎng)格基架(Stencil)，具有更小的數(shù)值耗散。其中WENO-SYMOO格式在光滑區(qū)逼近8階精度的中心差分格式，而WENO-SYMBO格式在光滑區(qū)逼近4階精度的帶寬優(yōu)化格式。數(shù)值測試顯示，2個格式差異不大。由于WENO-SYMBO/SYMOO的一個子格式使用了純單側(cè)差分格式，為了保持計算的穩(wěn)定性，需要對該子格式的權(quán)重進行限制， Martin等[10]強制該子格式的權(quán)重不超過其他子格式。 后來，Wu 和Martin[11]對WENO-SYMBO格式進行了進一步改進，首先通過變差函數(shù)判斷當(dāng)?shù)氐墓饣?。?dāng)總變差極大值與極小值的比值小于5時，認為流場光滑，關(guān)閉WENO格式的加權(quán)機制，從而使用線性權(quán)重。這時WENO格式退化為理想情況下的線性格式，從而進一步減小數(shù)值耗散及計算量。作者對WENO-SYMBO格式與WENO-JS格式進行了測試對比，發(fā)現(xiàn)WENO-SYMBO格式具有更小的數(shù)值耗散及更高的小尺度分辨率，但其魯棒性不如WENO-JS格式。

Sun等[12]對WENO格式的線性格式進行了優(yōu)化。他們發(fā)現(xiàn)，對于含有2個自由參數(shù)的線性格式，可以通過參數(shù)重組，使得其中一個參數(shù)控制格式的色散誤差，而另外一個參數(shù)控制格式的耗散誤差。在此基礎(chǔ)上，可以對格式的色散、耗散特性進行獨立優(yōu)化，從而大幅降低了優(yōu)化過程的復(fù)雜程度。在此基礎(chǔ)上他們對WENO格式的基礎(chǔ)格式(線性格式)進行了優(yōu)化，使得該格式具有最小的色散誤差(即色散分辨率最高)，而耗散誤差則可進行調(diào)整，因而該格式被稱為“色散最小、耗散可控的WENO(MDCD-WENO)格式。他們給出了低耗散、中耗散以及高耗散的3組MDCD-WENO格式，用戶可以根據(jù)計算的問題進行選擇。通常情況下，耗散越高則穩(wěn)定性越好。

WENO混合方法(Hybrid-WENO)也是目前常用的一種方法，該方法首先對當(dāng)?shù)氐墓饣赃M行判斷，如果足夠光滑，則切換到低耗散、低計算量的線性差分格式，反之，則切換到WENO格式。其中光滑區(qū)的線性格式可以是中心格式也可以是迎風(fēng)格式，可以是緊致格式也可以是非緊致格式。如果光滑區(qū)采用中心格式，則需要對其中的對流項進行特殊處理，以抑制混淆誤差，同時，還需要進行濾波或添加人工黏性以保持計算的穩(wěn)定。Hhybrid-WENO需要在WENO格式及線性格式之間進行切換，而硬性的切換會破壞差分的光滑性，從而造成一定的誤差。Ren 等[13]采用2種格式(WENO及線性格式)加權(quán)的方法來代替硬性切換的方法，其中權(quán)重根據(jù)當(dāng)?shù)氐墓饣潭仍?和1之間過渡，因而格式可保持整體上的光滑性。

1.1.2 高階保單調(diào)差分格式

通常的做法：如果fi+1/2超出了該區(qū)間的上限或者下限，則強制fi+1/2取該區(qū)間的上限或者下限值[5]，即

(1)

圖1為Shu-Osher問題的密度分布圖(200網(wǎng)格點)，可以看出，相同網(wǎng)格點上，6階精度OMP格式(OMP6)比7階精度的WENO格式(WENO7)具有更好的小尺度分辨率[14]。

圖1 Shu-Osher問題的密度分布[14]Fig.1 Density distribution of Shu-Osher problem[14]

在OMP格式的基礎(chǔ)上，Li和Fu[15]進一步優(yōu)化了格式的耗散，構(gòu)造了自適應(yīng)耗散的保單調(diào)(OMP-AD)格式。該格式的耗散可根據(jù)函數(shù)的光滑程度自行調(diào)節(jié)，在光滑區(qū)，格式具有非常低的數(shù)值耗散，而在間斷或大梯度區(qū)，格式具有較高的數(shù)值耗散以保持計算的穩(wěn)定。

OMP-AD格式為

(2)

(3)

式中：fj+k(k=-2～2)為j點周圍物理量的值。

1.1.3 WCNS格式

Deng和Zhang[6]利用緊致及加權(quán)思想，構(gòu)造出一類高階加權(quán)緊致非線性(WCNS)格式 。非線性譜分析顯示，該格式具有非常好的色散誤差特性及較低的數(shù)值耗散水平[16]。目前，該格式除了被用于簡單幾何外形的高分辨率計算，還被推廣到復(fù)雜工程問題的計算中[17]。該格式的對稱特征可以使得其易于在復(fù)雜網(wǎng)格上實現(xiàn)幾何守恒律[18]，從而盡量避免因為網(wǎng)格畸變而造成的附加誤差。

1.1.4 群速度控制格式

Ma和Fu[7-8]提出一種利用群速度控制以抑制間斷附近數(shù)值振蕩的激波捕捉格式。該格式在間斷兩側(cè)分別使用具有不同色散誤差特性的差分格式(“快格式”和“慢格式”)，以避免數(shù)值振蕩向間斷兩側(cè)擴展。由于群速度控制的主要思想是利用數(shù)值色散來抑制振蕩，因而該格式僅需使用較小的數(shù)值黏性即可抑制數(shù)值振蕩，因而格式的總體耗散較低。之后，He等[19]利用該思想提出了高分辨率的加權(quán)群速度控制(WGVC)格式，并通過與WENO格式混合提升了格式的魯棒性。

圖2為R-T不穩(wěn)定性問題中的密度分布圖(計算網(wǎng)格為121×481)[19]?？梢钥闯?，WGVC-WENO混合格式比同階精度的WENO格式具有更高的尺度分辨率。

圖2 R-T不穩(wěn)定性問題的密度分布[19]Fig.2 Density distribution of R-T instability problem[19]

2 穩(wěn)定計算技術(shù)

對于超聲速、高超聲速湍流的高分辨率數(shù)值模擬，如何保持計算穩(wěn)定性是非常重要的。計算發(fā)現(xiàn)，當(dāng)馬赫數(shù)很高時，用高精度格式計算很容易出現(xiàn)負壓力或負溫度，從而造成計算失敗?？蓧嚎s流動計算時大多從守恒方程出發(fā)，利用質(zhì)量及動量方程計算出速度從而得到動能，利用總能量方程計算出總能量，由總能量減去動能得到內(nèi)能，從而計算出溫度及壓力。而當(dāng)馬赫數(shù)很高時，動能非常接近總能量，受數(shù)值誤差的影響，計算出的動能有可能超過總能量，從而計算出的內(nèi)能為負，造成計算失敗。尤其是在物理量劇烈波動的區(qū)域，非常容易出現(xiàn)負壓力和負溫度。在壁湍流的DNS中，轉(zhuǎn)捩區(qū)物理量脈動強烈，更容易出現(xiàn)問題。另外在一些膨脹區(qū)，本身溫度就很低，當(dāng)出現(xiàn)數(shù)值波動時，也容易出現(xiàn)負值問題。因而，高馬赫數(shù)湍流的數(shù)值模擬必須采用某些方法以保持計算穩(wěn)定。

2.1 對流項的處理

(4)

(5)

(6)

式中：ρ為流場密度；ui為流動速度；φ為流場中的物理量，可以是速度，也可以是總能量或總焓等物理量。

如果采用迎風(fēng)格式離散，由于格式中含有數(shù)值耗散，對混淆誤差有較強的抑制作用，因而無需采用上述格式。

2.2 濾波

濾波技術(shù)可以濾掉流場中的高波數(shù)成分，從而可以保持計算的穩(wěn)定。通常每隔一定時間步，對流場進行一次濾波，濾波公式為

(7)

α=0時，為顯式(非緊致型)濾波；α≠0時，為隱式(緊致型)濾波。具體系數(shù)可見Pirozzoli[2]的綜述文章。濾波格式中的系數(shù)決定了濾波尺度，濾波尺度越大，穩(wěn)定性越強，但帶來的數(shù)值耗散也越嚴重。濾波不足會導(dǎo)致計算不穩(wěn)定，而濾波過度則小尺度流場被過度耗散，從而使計算分辨率下降。恰當(dāng)使用濾波非常關(guān)鍵，濾波的尺度、區(qū)域以及使用頻率都需要精確控制。目前還沒有通用的濾波準則，如何恰當(dāng)使用濾波大多還要靠經(jīng)驗。

2.3 保正性格式

在計算流體力學(xué)中，保正性就是保證密度、壓力、溫度在計算過程中不會出現(xiàn)負值。采用保正性格式可以保證計算穩(wěn)定，但不幸的是，保正性格式與單調(diào)格式一樣，很難推廣到高階精度。真正意義的保正性格式目前只有一階精度(如一階精度的Godnov類格式[20])，某些限制下(如嚴格的CFL(Courant-Friedrichs-Lewy)條件限制)，某些二階精度格式也能做到保正性，而更高階精度的保正性格式很難構(gòu)造。混合技術(shù)[21](光滑區(qū)采用高階格式，不穩(wěn)定的區(qū)域切換為保正性格式)是目前解決該問題的常用方法。

2.4 間斷識別器

由于高精度與計算穩(wěn)定性很難兼顧，很多數(shù)值計算采用了混合方法，在光滑區(qū)使用低耗散的高精度格式，在容易出問題的區(qū)域(間斷、大梯度區(qū)或劇烈波動區(qū))使用高穩(wěn)定性的差分格式(如WENO、保單調(diào)、保正性格式，或較大的人工黏性等)。這種情況下，合適的間斷識別器是十分必要的。Wu和Martin[11]采用當(dāng)?shù)啬０迳系目傋儾頣Vk作為激波識別器，當(dāng)其幅值大于一定數(shù)值時，就判斷當(dāng)?shù)貫殚g斷區(qū)，從而切換到捕捉激波的WENO-SYMBO格式。Darian 等[22]利用量階分析，構(gòu)造了一種激波識別器，用于實現(xiàn)2階濾波與高階濾波之間的切換。Shen和Zha[23]提出了一種基于WENO-Z格式[9]的關(guān)鍵參數(shù)τ5的間斷識別器為

τ5≤min(β0,β1,β2)

(8)

式中：βk(k=0，1，2)為各模板上的光滑因子[9]。是否滿足式(8)可作為區(qū)分光滑區(qū)與間斷區(qū)的判據(jù)。在此基礎(chǔ)上，Shen和Zha[23]提出了緊致格式與WENO格式的混合格式。

2.5 “壞點”的處理

對于具有強間斷，且有很強剛性源項的流動(例如高速化學(xué)反應(yīng)的流動)，流場的波動會被源項放大。例如，當(dāng)化學(xué)反應(yīng)對溫度非常敏感時，溫度場很小的波動就會造成化學(xué)反應(yīng)速率的劇烈變化，從而大幅放大振蕩。這時，僅采用激波捕捉格式可能無法抑制數(shù)值振蕩，從而使得計算發(fā)散。而簡單的濾波也很難處理這種情況。針對這種情況，Kotov等[24]構(gòu)造了一種處理方法。該方法首先采用光滑性(或單調(diào)性)判據(jù)，識別出流場中的“壞點(TroubledCell)”。對于這些“壞點”，使用其左側(cè)或右側(cè)光滑區(qū)的值重構(gòu)出這些點上的物理量，例如，使用左側(cè)或右側(cè)光滑點ENO插值構(gòu)造這些“壞點”的值。再通過單調(diào)判據(jù)選擇采用左側(cè)還是右側(cè)的重構(gòu)值。這樣，經(jīng)過重構(gòu)之后，物理量重新恢復(fù)了光滑性或單調(diào)性，從而避免了振蕩發(fā)生。

3 高超聲速湍流直接數(shù)值模擬算例

與亞、跨及超聲速湍流相比，高超聲速湍流更為復(fù)雜。高溫真實氣體效應(yīng)、壁溫效應(yīng)、熵層效應(yīng)都會對流動特性產(chǎn)生明顯的影響。此外，當(dāng)馬赫數(shù)足夠高時，Mack第2模態(tài)(甚至更高模態(tài))的不穩(wěn)定擾動波將被激發(fā)，導(dǎo)致高超聲速流動的轉(zhuǎn)捩機制更為復(fù)雜，轉(zhuǎn)捩預(yù)測的難度更大。由于計算復(fù)雜性，高超聲速湍流DNS研究直到最近幾年才開展起來。雖然起步晚，但由于航空航天等領(lǐng)域的強勁需求，近幾年發(fā)展非常迅速。

3.1 壓縮性效應(yīng)

由于壁面影響，相對于同樣來流條件下的自由湍流，壁湍流的壓縮性效應(yīng)更弱一些。通常認為，當(dāng)來流馬赫數(shù)不超過5時，壁湍流適用于Morkovin理論[25]。這時壓縮性效應(yīng)不是很強，而且壓縮性主要體現(xiàn)在對平均密度等平均量的影響上，對湍流脈動量影響較小。這種情況下，可以通過某些參數(shù)變換(如VanDirst變換)，將不可壓縮湍流的統(tǒng)計規(guī)律(如平均速度的對數(shù)律)推廣到可壓縮湍流。而對于馬赫數(shù)超過6的高超聲速壁湍流，壓縮性效應(yīng)則很難被完全忽略，尤其是面對低壁溫時的情況。

文獻[26]～文獻[29]對高超聲速湍流邊界層進行了一系列DNS研究， 包括壁溫效應(yīng)、馬赫數(shù)效應(yīng)以及高焓效應(yīng)。文獻[29] 對來流馬赫數(shù)為0.3～12.0，來流雷諾數(shù)(以邊界層厚度定義)約為1 500的平板邊界層湍流進行了DNS，其壁溫接近恢復(fù)溫度或接近絕熱壁溫。DNS結(jié)果顯示，對于平均速度、平均溫度以及湍流度等主要統(tǒng)計特征而言，壓縮性效應(yīng)并不明顯：VanDrist變換后的平均速度剖面與不可壓縮湍流的對數(shù)律吻合很好，平均溫度也符合Walz方程(Crocco關(guān)系)。而歸一化的湍流度與低速平板邊界層的分布之間差別也不大。如果采用半局部坐標度量[30]，不同馬赫數(shù)湍流的湍能生成、耗散以及輸運的分布也基本吻合。其中，半局部坐標采用壁面處摩擦應(yīng)力、當(dāng)?shù)孛芏纫约爱?dāng)?shù)仞ば远x的量為

(9)

這些研究表明，在較高壁溫(如接近絕熱壁)的情況下，Morkovin理論甚至可以部分拓展到馬赫數(shù)為12這樣的高超聲速邊界層。

Liang和Li[31]對馬赫數(shù)為3、6和8的平板邊界層湍流進行了DNS研究，其計算參數(shù)中包含了強冷壁面的工況(壁面溫度低于恢復(fù)溫度的0.5倍)。計算結(jié)果表明，在強冷壁情況下，平均量明顯偏離Morkovin理論的預(yù)測。這說明降低壁溫可明顯增強壓縮性效應(yīng)。此外，文獻[31]還對強雷諾比擬進行了修正，給出了適用于高馬赫數(shù)及強冷壁面的雷諾比擬。

Li等[32]對馬赫數(shù)為6的平板湍流進行了DNS，發(fā)現(xiàn)高馬赫數(shù)平板湍流的擬序結(jié)構(gòu)與低馬赫數(shù)平板有很大差別，高馬赫數(shù)下發(fā)卡渦很少出現(xiàn)，擬序結(jié)構(gòu)以準流向渦為主。

3.2 高焓及高溫真實氣體效應(yīng)湍流的DNS

對于高焓高超聲速流動，氣動加熱會使邊界層溫度升高，從而產(chǎn)生高溫真實氣體效應(yīng)。當(dāng)溫度足夠高時，空氣的比熱將不再保持常數(shù)，這時完全氣體假設(shè)將不再適用。當(dāng)溫度繼續(xù)升高時，空氣中的O2和N2將發(fā)生解離及化學(xué)反應(yīng)，從而形成復(fù)雜的化學(xué)非平衡流動。 高溫真實氣體湍流非常復(fù)雜，以往該流動的數(shù)值模擬多以雷諾平均模擬(RANS)為主，高分辨率模擬，(DNS和大渦模擬(LES))直到近幾年才開展起來。

Duan和Martin[27]以來流馬赫數(shù)為21,高度為30km的來流的楔形體邊界層(半楔角分別為8°和35°)為背景流動，進行了邊界層某一局部區(qū)域的DNS，對應(yīng)的局部雷諾數(shù)分別為3.5×105(35°半楔角工況)和4.5×106(8°半楔角工況)。計算結(jié)果顯示，流場表現(xiàn)出了明顯的高溫真實氣體效應(yīng)，近壁區(qū)的比熱比γ在1.28～1.30之間，明顯偏離了完全氣體的值(γ=1.4)。其中35°半楔角工況下，邊界層氧原子的質(zhì)量比分接近20%，表明該條件下大部分O2已經(jīng)解離。根據(jù)Duan和Martin的研究結(jié)果，高焓高超聲速湍流邊界層的統(tǒng)計特性與低速情況有很大差別，壓縮性效應(yīng)更為明顯，流動與Morkovin假設(shè)偏差很大。

Chen和Li[33]對來流馬赫數(shù)為6和10的槽道湍流進行了DNS，計算考慮了高溫情況下比熱的變化以及O2及N2的化學(xué)反應(yīng)等高溫真實氣體效應(yīng)，并與同樣來流情況下完全氣體流動進行了對比。計算結(jié)果顯示，與完全氣體相比，真實氣體的平均溫度及湍動能明顯降低，說明高溫真實氣體效應(yīng)對湍流有一定抑制作用。

3.3 高超聲速流動的轉(zhuǎn)捩

當(dāng)馬赫數(shù)足夠高時，Mack第2模態(tài)擾動波將被激發(fā)甚至主導(dǎo)轉(zhuǎn)捩。由于存在多種不穩(wěn)定模態(tài)，因而高超聲速邊界層的轉(zhuǎn)捩過程更為復(fù)雜。文獻[34]～文獻[38]利用數(shù)值模擬研究了高超聲速鈍錐以及超聲速平板邊界層對來流擾動的感受性過程。對于該流動，自由來流的擾動首先要穿過頭激波，擾動與激波相互作用后，分解為熵波、渦波和聲波。其中熵波和渦波的傳播速度為當(dāng)?shù)芈暳魉賣,聲波的傳播速度為u+c與u-c, 其中c為當(dāng)?shù)芈曀?。傳播速度為u+c的聲波稱為“快聲波”，傳播速度為u-c的聲波為“慢聲波”。這些波被邊界層吸收的機制是不同的。

不穩(wěn)定TS(Tollmien-Schlichting)波的相速度通常會低于當(dāng)?shù)氐膶α魉俣?，很多計算結(jié)果顯示，相速度大約為邊界層外緣速度的0.9倍左右(不穩(wěn)定TS波相速度處于0.8倍的邊界層外緣速度到0.95倍的邊界外緣速度之間)。因而以u-c速度傳播的慢聲波有可能落入不穩(wěn)定TS波的相速度區(qū)間，從而可能與某些不穩(wěn)定TS波發(fā)生共振，激發(fā)起不穩(wěn)定的TS波。因而，超/高超聲邊界層對于慢聲波有較好的感受性。 但這并不意味著其他擾動波無法被邊界層吸收。根據(jù)Zhong 和Ma的研究[34], 快聲波也可通過特殊的機制激發(fā)出邊界層內(nèi)不穩(wěn)定波。 由于相速度差異，快聲波無法直接與不穩(wěn)定TS波形成共振，但可以與某些穩(wěn)定TS波發(fā)生共振。有些穩(wěn)定TS波在前緣附近的傳播速度可以達到流動速度的1.2倍左右, 因而可以和快聲波形成共振。通過共振，快聲波將能量傳遞給邊界層內(nèi)的穩(wěn)定TS波。 該穩(wěn)定TS波向下游傳播過程中振幅會逐漸衰減，但某些情況下，其衰減率并不是很高，因而可以將能量攜帶到較遠的下游。

隨著該波的傳播，由于邊界層厚度及平均剖面的變化，該穩(wěn)定擾動波的相速度逐漸降低，從而進入不穩(wěn)定TS波的相速度區(qū)間。在下游的某一位置，該穩(wěn)定擾動波與某一不穩(wěn)定擾動波達到了相同的相速度，進而形成共振，激發(fā)出不穩(wěn)定TS波。因而，即使存在與不穩(wěn)定TS波的相速度差異，快聲波也可通過某一穩(wěn)定TS波作為載體，激發(fā)出不穩(wěn)定的TS波，從而導(dǎo)致流動失穩(wěn)。

Li等[8,39-40]運用DNS，研究頭半徑為1 mm,半錐角為5°的小頭鈍錐的轉(zhuǎn)捩特性。其中來流馬赫數(shù)為6，以頭半徑定義的來流雷諾數(shù)為10 000, 攻角為1°。為了觸發(fā)轉(zhuǎn)捩，在距頭部90～100 mm處的壁面上添加1%幅值的隨機吹吸擾動(模擬壁面上的粗糙單元)。

計算結(jié)果顯示，在擾動區(qū)下游，Mack第1模態(tài)及第2模態(tài)的擾動波均被激發(fā)出來，但第2模態(tài)的擾動波增長更快，并主導(dǎo)了下游的轉(zhuǎn)捩過程。計算還發(fā)現(xiàn)，錐身上的轉(zhuǎn)捩線呈現(xiàn)非單調(diào)特征，如圖3所示，其中Cf為摩擦阻力系數(shù)。從背風(fēng)面算起周角20°子午面附近的區(qū)域轉(zhuǎn)捩大幅推遲，其機理尚有待進一步分析。

圖3 錐身上的摩擦阻力系數(shù)分布及轉(zhuǎn)捩線[40]Fig.3 Skin fraction coefficient and transition line on surface of blunt cone[40]

針對高超聲速鈍錐轉(zhuǎn)捩的預(yù)測問題，Su和Zhou[41]基于擾動波的發(fā)展規(guī)律，對傳統(tǒng)的eN方法進行了修正。傳統(tǒng)的eN方法假設(shè)擾動波按照線性穩(wěn)定性理論進行指數(shù)增長，從進入非穩(wěn)定狀態(tài)開始，振幅增長到eN倍時就發(fā)生了轉(zhuǎn)捩(N通常為8～11之間，由實驗或經(jīng)驗確定)。而修正的eN方法假設(shè)擾動波振幅增長到固定值(例如主流的1%或2%振幅)就發(fā)生轉(zhuǎn)捩，該方法既考慮了擾動波的增長過程，又考慮了前期的衰減過程。該方法的預(yù)測結(jié)果與Li等[40]的數(shù)值結(jié)果吻合較好。

最近，Sandham等[42]進行了馬赫數(shù)為6的來流條件下平板湍流以及平板湍流-入射激波干擾的DNS及實驗研究，并將DNS與實驗進行了交叉比對。通過DNS及實驗，探討了雷諾數(shù)效應(yīng)、擾動強度、壁面冷卻以及入射激波對轉(zhuǎn)捩的影響。研究結(jié)果顯示，在該條件下，Mack第2模態(tài)不穩(wěn)定波主導(dǎo)了轉(zhuǎn)捩過程。

4 面向可壓縮湍流的高精度計算軟件

近年來計算流體力學(xué)(CFD)軟件得到了迅速發(fā)展，目前比較常用的商用CFD軟件有FLUENT, STAR-CD、CFX、CFD-ACE以及PHONEX等。此外，國外很多科研機構(gòu)開發(fā)的CFD軟件也得到了廣泛應(yīng)用，例如美國國家航空航天局(NASA)的CFL3D, 德國宇航中心(DLR)的TAU等。這些CFD軟件在工程應(yīng)用中發(fā)揮了重要作用。但目前流行的CFD軟件大多以工程應(yīng)用為目標，強調(diào)程序及方法的普適性及魯棒性。這些軟件多以2階精度有限體積方法為主要算法，在計算精度及分辨率方面有一定不足，因而很難被用于像DNS這樣的高分辨率計算。

面向可壓縮湍流高分辨率數(shù)值模擬方面的需求，作者開發(fā)了一套高精度CFD軟件— OpenCFD。該軟件包括2個求解器：高精度差分求解器(Open CFD code for Scientific Computing，OpenCFD-SC)[43]以及多塊結(jié)構(gòu)網(wǎng)格有限體積求解器 (Open CFD code for Engineering Computing，OpenCFD-EC)[44]。其中OpenCFD-SC主要用于進行可壓縮湍流的高分辨率數(shù)值計算(DNS及大渦模擬)。

圖4為OpenCFD-SC求解器的主要框架。該求解器基于高精度有限差分法，其核心為差分庫，其中既包含了作者課題組開發(fā)的高精度差分格式(如加權(quán)群速度控制格式、優(yōu)化保單調(diào)差分格式等)，又包含了像WENO系列格式(包括WENO-JS及目前流行的多種改進WENO格式)等流行的差分格式，還包括線性迎風(fēng)和中心差分格式等常規(guī)高精度格式。該求解器的特點是精度高、耗散小，具有很高的數(shù)值分辨率。作者近年來的DNS算例(包括超聲速及高超聲速湍流DNS)大多采用該求解器計算[8]。OpenCFD-SC的另一個特點是具有很強的并行可擴展性。

圖5為OpenCFD-SC在國產(chǎn)某眾核并行計算機上的并行加速測試情況，本次測試最多使用了983 040個CPU核心(包括主核和從核)?？梢钥闯觯?dāng)使用524 288個CPU核心時(包括主核及從核)，并行效率仍可達59%(以32 768個CPU核心為基準)，顯示了該軟件具有非常強的大規(guī)模擴展性。

圖4 OpenCFD-SC的主要框架Fig.4 Main frame of OpenCFD-SC

圖5 OpenCFD 軟件的大規(guī)模并行測試加速比情況Fig.5 Speed-up ratio of OpenCFD in a large-scale parallel test

5 結(jié) 論

回顧了高超聲速湍流數(shù)值模擬方面的進展，包括數(shù)值方法以及直接數(shù)值模擬(DNS)研究的進展。同時，也簡要回顧了目前高超聲速湍流邊界層DNS方面的研究進展。

1)在數(shù)值格式方面，WENO類格式仍是目前高超聲速湍流數(shù)值模擬的主要方法，為了避免WENO格式耗散過大的不足，目前出現(xiàn)了很多優(yōu)化/改進版的WENO格式，如WENO-Z,WENO-SYMBO等。此外WENO混合格式也是目前的一種主要方法。與WENO格式相比，高階保單調(diào)(MP)格式具有格式簡單、計算量小、耗散低等優(yōu)點，但其魯棒性有待進一步改進。此外，中國學(xué)者構(gòu)造的高階加權(quán)緊致非線性(WCNS)格式,群速度控制(GVC)格式等也是很好的選擇。

2)當(dāng)馬赫數(shù)足夠高時，為了保持計算穩(wěn)定，防止出現(xiàn)負溫度或壓力，很有必要采用某些技術(shù)以保證計算穩(wěn)定。這時，可以采用可抑制混淆誤差的Skew-Symmetry形式的對流項、采用保正性格式以及數(shù)值濾波來保持計算穩(wěn)定。同時，為了盡量減小數(shù)值耗散，以激波識別器為基礎(chǔ)，采用分區(qū)混合計算(光滑區(qū)采用低耗散方法，間斷區(qū)或易出現(xiàn)問題的區(qū)域采用高魯棒性方法)也是很好的選擇。

3)由于計算難度很大，高超聲速湍流邊界層的DNS研究近幾年才剛剛開展，但進展很快。這些DNS在壓縮性效應(yīng)、壁溫效應(yīng)、真實氣體效應(yīng)以及高超聲速轉(zhuǎn)捩方面進行了有益探索。研究發(fā)現(xiàn)，相對于自由湍流、壁湍流的壓縮性效應(yīng)較弱。對于高壁溫(接近恢復(fù)溫度)的壁湍流、Morkovin假設(shè)甚至可以部分推廣到來流馬赫數(shù)為12這樣的高超聲速流動。而對于低壁溫的情況，壓縮性效應(yīng)則明顯增強。隨著馬赫數(shù)升高，Mack第2模態(tài)擾動波被激發(fā)，引起高超聲速壁湍流轉(zhuǎn)捩更加復(fù)雜。

近年來航天領(lǐng)域的強勁需求牽引人們對高超聲速湍流進行更加深入地探索。作為研究湍流的重要工具，DNS在高超聲速湍流研究方面將會發(fā)揮更加重要的作用。相信近年內(nèi)在高精度、高魯棒性數(shù)值方法以及高超聲速湍流DNS方面會有更快的發(fā)展。

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

[2] Pirozzoli S. Numerical methods for high-speed flows[J]. Annual Review of Fluid Mechanics, 2011,43:163-194.

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

[4] Jiang G S, Shu C W. Efficient implementation of weighted ENO schemes[J]. Journal of Computational Physics, 1996, 126(1): 202-228.

[5] Suresh A, Huynh H T. Accurate monotonicity-preserving schemes with Runge-Kutta time stepping[J]. Journal of Computational Physics, 1997, 136(1): 83-99.

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

[7] Ma Y W, Fu D X. Fourth order accurate compact scheme with group velocity control (GVC)[J].Science in China Series A: Mathematics Physics & Astronomy, 2001, 44(9): 1197-1204.

[8] Fu D X, Ma Y W, Li X L, et al. Direct numerical simulation of compressible turbulence[M]. Beijing: Science Press, 2011: 144-164 (in Chinese). 傅德薰， 馬延文， 李新亮， 等. 可壓縮湍流直接數(shù)值模擬 [M]. 北京： 科學(xué)出版社， 2011: 144-164.

[9] Borges R, Carmona M, Costa B, et al. An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws[J]. Journal of Computational Physics, 2008, 227(6): 3191-3211.

[10] Martin M P, Taylor E M, Wu M, et al. A bandwidth-optimized WENO scheme for the effective direct numerical simulation of compressible turbulence[J]. Journal of Computational Physics, 2006, 220(1): 270-289.

[11] Wu M, Martin M P. Direct numerical simulation of supersonic turbulent boundary layer over a compression ramp [J]. AIAA Journal, 2007, 45(4): 879-889.

[12] Sun Z S, Ren Y X, Larricq C, et al. A class of finite difference schemes with low dispersion and controllable dissipation for DNS of compressible turbulence[J]. Journal of Computational Physics, 2011, 230(12): 4616-4635.

[13] 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.

[14] Li X L, Leng Y, He Z W. Optimized sixth-order monotonicity-preserving scheme by nonlinear spectral analysis [J]. International Journal for Numerical Methods in Fluids, 2013, 73(6): 560-577.

[15] Li X L, Fu D X. Optimized MP scheme with adaptive dissipation and DNS of supersonic turbulent flows in DLR scramjet intake[C]∥Eighth International Conference on Computational Fluid Dynamics. Chengdu: ICCFD, 2014.

[16] Tu G H, Deng X G, Mao M L. Spectral property comparison of fifth-order nonlinear WCNS and WENO difference schemes[J]. Acta Aerodynamics Sinica, 2012, 30(6): 709-712 (in Chinese). 涂國華， 鄧小剛， 毛枚良. 5階非線性WCNS和WENO差分格式頻譜特性比較[J]. 空氣動力學(xué)報， 2012, 30(6)：709-712.

[17] Deng X G, Mao M L, Tu G H, et al. High-order and high accurate CFD methods and their applications for complex grid problems[J]. Communications in Computational Physics, 2011, 11(4): 1081-1102.

[18] Deng X G, Mao M L, Tu G H, et al. Geometric conservation law and applications to high-order finite difference schemes with stationary grids[J]. Journal of Computational Physics, 2011, 230(4): 1100-1115.

[19] He Z W, Li X L, Liang X. Nonlinear spectral-like schemes for hybrid schemes[J]. Science China: Physics, Mechanics & Astronomy, 2014, 57(4): 753-763.

[20] Toro E F. Riemann solvers and numerical methods for fluid dynamics: a practical introduction[M]. Berlin: Springer-Verlag, 2009: 174-184.

[21] Hu X Y，Adams N A，Shu C W. Positivity-preserving method for high-order conservative schemes solving compressible Euler equations[J]. Journal of Computational Physics, 2013, 242: 169-180.

[22] Darian H M, Esfahanian V, Hejranfar K. A shock-detecting sensor for filtering of high-order compact finite difference schemes[J]. Journal of Computational Physics, 2011, 230(3): 494-514.

[23] Shen Y Q, Zha G C. Generalized finite compact difference scheme for shock/complex flowfield interaction[J]. Journal of Computational Physics, 2011, 230(12): 4419-4436.

[24] Kotov D V, Yee H C, Sjogreen B, et al. Performance of four high-order shock-capturing schemes for stiff source terms with discontinuities: preliminary results[R]. Center for Turbulence Research Annual Research Briefs, 2011: 393-403.

[25] Morkovin M V. Effects on compressibility on turbulent flows[M]∥Favre A J. Mecanique de la turbulence. Paris: CNRS, 1962: 367-380.

[26] Martin M P. Direct numerical simulation of hypersonic turbulent boundary layers, Part 1. Initialization and comparison with experiments[J]. Journal of Fluid Mechanics, 2007,570: 347-364.

[27] Duan L, Martin M P. Direct numerical simulation of hypersonic turbulent boundary layers, Part 4. Effect of high enthalpy[J]. Journal of Fluid Mechanics, 2011, 684: 25-59.

[28] Duan L, Beekman I, Martin M P. Direct numerical simulation of hypersonic turbulent boundary layers, Part 2. Effect of wall temperature[J]. Journal of Fluid Mechanics, 2010, 655: 419-445.

[29] Duan L, Beekman I, Martin M P. Direct numerical simulation of hypersonic turbulent boundary layers, Part 3. Effect of Mach number[J]. Journal of Fluid Mechanics, 2011, 672: 245-267.

[30] Huang P G, Coleman G N, Bradshaw P. Compressible turbulent channel flows: DNS results and modeling[J]. Journal of Fluid Mechanics, 1995, 305: 185-218.

[31] Liang X, Li X L. DNS of a spatially evolving hypersonic turbulent boundary layer at Mach 8[J]. Science China: Physics Mechanics and Astronomy, 2013, 56(7):1408-1418.

[32] Li X L, Fu D X, Ma Y W. Direct numerical simulation of a spatially evolving supersonic turbulent boundary layer atMa=6[J]. Chinese Physics Letters, 2006, 23(6): 1519-1522.

[33] Chen X P， Li X L. Direct numerical simulation of chemical non-equilibrium turbulent flow[J]. Chinese Physical Letters, 2013, 30(6): 064702.

[34] Zhong X, Ma Y. Boundary-layer receptivity of Mach 7.99 flow over a blunt cone to free-stream acoustic waves[J]. Journal of Fluid Mechanics, 2006, 556: 55-103.

[35] Ma Y, Zhong X. Receptivity of a supersonic boundary layer over a flat plate. Part 1. Wave structures and interactions[J]. Journal of Fluid Mechanics, 2003,488:31-78.

[36] Ma Y, Zhong X. Receptivity of a supersonic boundary layer over a flat plate， Part 2. Receptivity to freestream sound[J]. Journal of Fluid Mechanics, 2003,488:79-121.

[37] Ma Y, Zhong X. Receptivity of a supersonic boundary layer over a flat plate， Part 3. Effects of different types of free-stream disturbances[J]. Journal of Fluid Mechanics, 2005,532: 63-109.

[38] Zhong X, Wang X. Direct numerical simulation on the receptivity, instability and transition of hypersonic boundary layers[J]. Annual Review of Fluid Mechanics, 2012,44:527-561.

[39] Li X L, Fu D X, Ma Y W. Direct numerical simulation of hypersonic boundary-layer transition over a blunt cone [J]. AIAA Journal, 2008,46(11): 2899-2913.

[40] 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.

[41] Su C H, Zhou H. Transition prediction of a hypersonic boundary layer over a cone at small angle of attack-with the improvement of eNmethod[J], Science in China G, Mechanics and Astronomy, 2009, 52 (1): 115-123.

[42] Sandham N D, Schulein E A, Wagner E A, et al. Transitional shock-wave/boundary-layer interactions in hypersonic flow[J]. Journal of Fluid Mechanics, 2014, 752:349-382.

[43] Li X L. OpenCFD-SC user’s manual[EB/OL]. (2011-4-12) [2014-9-20]. http:∥www.cfluid.com/bbs/viewthread.php?tid=89129&extra=page%3D1 (in Chinese). 李新亮. OpenCFD-SC用戶手冊[EB/OL]. (2011-4-12) [2014-9-20]. http:∥www.cfluid.com/bbs/viewthread.php?tid=89129&extra=page%3D1.

[44] Li X L. OpenCFD-EC User’s manual[EB/OL]. (2011-4-12) [2014-9-20]. http:∥www.cfluid.com/bbs/viewthread.php?tid=89129&extra=page%3D1 (in Chinese). 李新亮.OpenCFD-EC理論手冊[EB/OL]. (2011-4-12)[2014-9-20]. http:∥www.cfluid.com/bbs/viewthread.php?tid=91376&extra=page%3D1

Tel: 010-82543801

E-mail: lixl@imech.ac.cn

*Corresponding author. Tel.： 010-82543801 E-mail: lixl@imech.ac.cn

Direct numerical simulation techniques for hypersonic turbulent flows

LI Xinliang*

StateKeyLaboratoryofHigh-temperatureGasDynamics,InstituteofMechanics,ChineseAcademyofSciences,Beijing100190,China

The recent developments of high resolution schemes, especially, high-order and high-robustness shock-capture schemes, and direct numerical simulation (DNS) cases for hypersonic turbulent flows are reviewed in this paper. The numerical methods include the high-resolution shock-capture methods and the technique to stabilize computation for hypersonic flows, as well as, the developments of WENO and monotonicity preserving schemes. The DNS studies include the effects of compressibility, wall temperature and high-temperature real gas on the turbulent flows, and the studies of hypersonic transition flows are also reviewed briefly. Furthermore, an OpenCFD code developed by the author which is compressible and high-resolution, is addressed briefly

hypersonic; turbulence; direct numerical simulation; high-resolution numerical method; OpenCFD

2014-07-25； Revised： 2014-09-16； Accepted： 2014-09-20； Published online： 2014-10-31 17:05

s： National Natural Science Foundation of China (1372330, 11472278, 11472010, 91441103); National High Technology Research and Development Program (2012AA01A304); Chinese Academy of Sciences Innovation Programs(KJCX2-EW-J01，XXH12503-02-02-04)

2014-07-25; 退修日期： 2014-09-16; 錄用日期： 2014-09-20; 網(wǎng)絡(luò)出版時間： 2014-10-31 17：05

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

國家自然科學(xué)基金(1372330, 11472278, 11472010, 91441103);國家“863”計劃(2012AA01A304);中國科學(xué)院知識創(chuàng)新工程(KJCX2-EW-J01, XXH12503-02-02-04)

Li X L. Direct numerical simulation techniques for hypersonic turbulent flows[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(1): 147-158.李新亮． 高超聲速湍流直接數(shù)值模擬技術(shù)[J]． 航空學(xué)報, 2015, 36(1): 147-158.

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

10.7527/S1000-6893.2014.0233

V211.3; O357.5

A

1000-6893(2015)01-0147-12

李新亮 男， 博士， 研究員。主要研究方向： 計算流體力學(xué)， 湍流。

*通訊作者.Tel.： 010-82543801 E-mail: lixl@imech.ac.cn

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