基于非結(jié)構(gòu)／混合網(wǎng)格的高階精度DG／FV混合方法研究進(jìn)展

張來(lái)平，李明，劉偉，赫新，張涵信

（1．空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川 綿陽(yáng) 621000；

2．中國(guó)空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所，四川 綿陽(yáng) 621000）

基于非結(jié)構(gòu)／混合網(wǎng)格的高階精度DG／FV混合方法研究進(jìn)展

張來(lái)平1，2，李明2，劉偉1，2，赫新2，張涵信2

（1．空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川 綿陽(yáng) 621000；

2．中國(guó)空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所，四川 綿陽(yáng) 621000）

DG／FV混合方法因其具有緊致、易于推廣獲得高階格式及相比同階精度DG方法計(jì)算量、存儲(chǔ)量小等優(yōu)點(diǎn)，自提出以來(lái)已成功應(yīng)用于一維、二維標(biāo)量方程和Euler／N-S方程的求解。綜述了DG／FV混合方法的研究進(jìn)展，重點(diǎn)介紹了DG／FV混合方法的空間重構(gòu)算法、針對(duì)RANS方程的求解方法、隱式時(shí)間離散格式、數(shù)值色散耗散及穩(wěn)定性分析、計(jì)算量理論分析，并給出了系列粘性流算例的計(jì)算結(jié)果，包括用于驗(yàn)證混合方法數(shù)值精度的庫(kù)埃特流，以及方腔流、亞聲速剪切層、低速平板湍流、NACA0012翼型湍流繞流等。數(shù)值計(jì)算結(jié)果表明DG／FV混合方法達(dá)到了設(shè)計(jì)的精度階，且相比同階DG方法計(jì)算量減少約40%，而隱式方法能大幅提高定常流的收斂歷程，較顯式Runge-Kutta的收斂速度提高1～2個(gè)量級(jí)。

非結(jié)構(gòu)／混合網(wǎng)格；間斷Galerkin方法；有限體積方法；DG／FV混合方法；RANS方程

0 引 言

隨著計(jì)算機(jī)技術(shù)和計(jì)算流體力學(xué)（CFD）的發(fā)展，近幾十年來(lái)網(wǎng)格生成技術(shù)及數(shù)值計(jì)算方法取得了飛速的進(jìn)步。在網(wǎng)格生成技術(shù)中，非結(jié)構(gòu)網(wǎng)格具有適合復(fù)雜外形、方便網(wǎng)格自適應(yīng)等突出優(yōu)點(diǎn)，進(jìn)一步發(fā)展的混合網(wǎng)格技術(shù)部分克服了非結(jié)構(gòu)網(wǎng)格需要大量計(jì)算資源的不足，代表了網(wǎng)格技術(shù)的發(fā)展趨勢(shì)［1］。而對(duì)數(shù)值計(jì)算方法來(lái)說(shuō)，目前已知的絕大多數(shù)CED商業(yè)軟件和in-house工業(yè)應(yīng)用軟件均以二階精度的有限體積方法（FVM）為基礎(chǔ)，盡管其已成功解決了大量的工程實(shí)際問(wèn)題，但在CFD領(lǐng)域內(nèi)仍有許多問(wèn)題需要高階精度（三階或以上精度）方法才能較好地解決。這主要是因?yàn)榈碗A方法具有較大的數(shù)值耗散與色散，對(duì)一些非常復(fù)雜的流動(dòng)現(xiàn)象，如旋渦主導(dǎo)的流動(dòng)、分離、湍流等問(wèn)題，其通常難以給出精細(xì)的流場(chǎng)結(jié)構(gòu)。尤其對(duì)于湍流的大渦模擬（LES）、直接數(shù)值模擬（DNS）等問(wèn)題，低階方法固有的數(shù)值耗散可能大到掩蓋這些真實(shí)的物理粘性，采用高階方法已是CFD界的共識(shí)［2］。在計(jì)算氣動(dòng)聲學(xué)、計(jì)算電磁學(xué)等領(lǐng)域，需要對(duì)長(zhǎng)時(shí)間歷程的波傳播進(jìn)行高精度的數(shù)值模擬，低階方法在網(wǎng)格規(guī)模受計(jì)算機(jī)資源限制時(shí)往往無(wú)法準(zhǔn)確模擬流場(chǎng)的特性，帶來(lái)無(wú)法接受的計(jì)算結(jié)果，此時(shí)亦有必要使用高階方法。此外分析表明，對(duì)于誤差水平要求不高的問(wèn)題，低階方法可以用較小的計(jì)算代價(jià)來(lái)解決；而對(duì)誤差水平要求很高的問(wèn)題，使用高階方法計(jì)算代價(jià)更小，效率更高［2］。

近二十年來(lái)，基于非結(jié)構(gòu)／混合網(wǎng)格的高階精度計(jì)算方法發(fā)展迅速，CFD工作者已提出大量的計(jì)算方法，主要包括k-exact FVM［3-4］、間斷Galerkin方法（DGM）［5-6］、有限譜體積（SV）方法［7］、有限譜差分（SD）方法［8］、及將DG、SV、SD等方法統(tǒng)一在一個(gè)框架之內(nèi)的CPR（Correction Procedure via Reconstruction）方法［9］等等。更多相關(guān)內(nèi)容可以參考Ekaterinaris［10］、Wang［2］和張來(lái)平等［11］的綜述文章。

高階k-exact FV方法和以DG方法為代表的DG／SV／SD／CPR等方法都各具優(yōu)點(diǎn)，但仍有可以改進(jìn)的空間。如高階FVM需要擴(kuò)充模板來(lái)提高重構(gòu)精度，在非結(jié)構(gòu)網(wǎng)格上模板的搜尋擴(kuò)展很不方便并且方法不緊致，而DGM的計(jì)算量和存儲(chǔ)量非常大。由此構(gòu)造結(jié)合FV方法和DG方法優(yōu)點(diǎn)的混合方法是一種較好的選擇，其已受到許多學(xué)者的關(guān)注，并提出了多種混合方案。這些混合方法的核心思想是使用本單元和相鄰單元的多個(gè)自由度重構(gòu)一個(gè)更高階的分布，使用這個(gè)更高階的分布來(lái)得到本單元的較低階自由度信息的高階更新。

Cockburn等［12］為了進(jìn)一步提高DG方法精度最早提出對(duì)解使用重構(gòu)算法的思想，后來(lái)Ryan等［13］做了進(jìn)一步的發(fā)展。上述工作僅僅是在最后輸出解時(shí)使用重構(gòu)算法，因而可以認(rèn)為是DG方法的后處理技術(shù)。Dumbser和Munz首次提出從DG方法計(jì)算開(kāi)始時(shí)刻即使用線性重構(gòu)算子，并構(gòu)造出一類重構(gòu)的DG方法，命名為PNPM方法［14-15］；基于類似思想，Luo等提出了RDG［16-17］（Reconstructed DG）方法，構(gòu)造了3階RDG（P1P2）格式；張來(lái)平等構(gòu)造出一類DG／FV混合方法［18-22］；在原始CPR格式的基礎(chǔ)上，王志堅(jiān)等借鑒PNPM方法和DG／FV混合方法的思想，構(gòu)造了一系列PNPM-CPR格式［23-24］。這幾種混合方法均顯示出了很好的性能，具有每個(gè)自由度比FV方法和DG方法更高效、方法緊致、隱式時(shí)間離散內(nèi)存需求更低等很多優(yōu)點(diǎn)［11］。

張來(lái)平等在構(gòu)造DG／FV混合方法［18-22］時(shí)提出了“靜態(tài)重構(gòu)”和“動(dòng)態(tài)重構(gòu)”的概念，并基于靜動(dòng)態(tài)“混合重構(gòu)”的思想構(gòu)造了一類三階以上精度的DG／FV混合格式，其已應(yīng)用于三角形／四邊形混合網(wǎng)格下的標(biāo)量方程和Euler／N-S方程的數(shù)值模擬，大量的數(shù)值算例表明該方法達(dá)到了設(shè)計(jì)精度，且相比同階精度DG方法具有更高的計(jì)算效率［25］。

本文對(duì)DG／FV混合方法研究進(jìn)展進(jìn)行了簡(jiǎn)要綜述，重點(diǎn)介紹了“混合重構(gòu)”的基本思想、針對(duì)RANS方程的求解方法、隱式時(shí)間離散格式、色散耗散特性及穩(wěn)定性條件、計(jì)算量理論分析等。在此基礎(chǔ)上，給出了系列典型算例的層流和湍流計(jì)算結(jié)果，并與相應(yīng)的精確解和文獻(xiàn)結(jié)果進(jìn)行了比較，以驗(yàn)證DG／FV方法對(duì)粘性流動(dòng)問(wèn)題模擬的適用性。

1 DG／FV混合方法

1.1 控制方程

在CED中控制方程一般為Navier-Stokes方程，它的守恒形式可簡(jiǎn)寫(xiě)為：

其中U為守恒變量；Fc為無(wú)粘（對(duì)流）通量；Fv為粘性通量，它們的具體形式可參考文獻(xiàn)［26］。

對(duì)于湍流計(jì)算我們使用雷諾平均N-S（RANS）方程，本文使用Spalart-Allmaras（S-A）一方程湍流模型。S-A模型的湍流渦粘性表達(dá)式為νt＝fv1，其中為模型因變量，其滿足的微分方程為：

上式中的各種系數(shù)可參考文獻(xiàn)［26］。本文中對(duì)式（1）和式（2）使用解耦的方式求解，并且目前僅初步使用低階方式求解式（2）。

1.2 靜態(tài)重構(gòu)和動(dòng)態(tài)重構(gòu)

本質(zhì)上說(shuō)，數(shù)值格式的構(gòu)造過(guò)程就是離散和重構(gòu)的過(guò)程。離散即利用網(wǎng)格技術(shù)將計(jì)算域分解為離散網(wǎng)格單元；重構(gòu)可以看成把解耦的信息進(jìn)行耦合，把丟失的信息“還原”。不同的“還原”方式導(dǎo)致了不同的格式，因此重構(gòu)技術(shù)對(duì)格式的構(gòu)造起到至關(guān)重要的作用。在DG方法中，解一般認(rèn)為是跨單元間斷的多項(xiàng)式。解多項(xiàng)式系數(shù)的約束關(guān)系是時(shí)間相關(guān)的，其隨“時(shí)間”同步推進(jìn)計(jì)算，因此我們稱之為“動(dòng)態(tài)重構(gòu)”（或“時(shí)間相關(guān)重構(gòu)”）。解多項(xiàng)式的信息由初邊值條件和控制方程經(jīng)過(guò)Galerkin有限元方法“提煉”而來(lái)。而在k-exact FV方法中，解只有單元平均值隨時(shí)間推進(jìn)更新，且更新時(shí)使用一個(gè)重構(gòu)的解高階多項(xiàng)式分布來(lái)計(jì)算數(shù)值通量。解高階分布的系數(shù)由相鄰單元的單元平均值插值而來(lái)，不與單元平均值一起推進(jìn)計(jì)算，可以看成是時(shí)間上的一種“后處理”過(guò)程。重構(gòu)多項(xiàng)式的高階信息來(lái)源于鄰近單元，因此我們稱之為“靜態(tài)重構(gòu)”（或“網(wǎng)格相關(guān)重構(gòu)”）。

在本文DG／FV混合方法中，對(duì)每個(gè)單元構(gòu)造其模板，并使用模板中每個(gè)單元上物理量初始較低的PDG階多項(xiàng)式分布重構(gòu)出更高階的PFV階多項(xiàng)式分布（PFV＞PDG），此過(guò)程稱為“靜態(tài)重構(gòu)”；使用此PFV階高階分布來(lái)計(jì)算數(shù)值通量和數(shù)值積分，從而對(duì)單元中物理量的PDG階多項(xiàng)式分布系數(shù)使用常規(guī)的DG方法存儲(chǔ)和時(shí)間推進(jìn)（此稱為“動(dòng)態(tài)重構(gòu)”），得到下一時(shí)刻的PDG階分布系數(shù)。理論分析和數(shù)值計(jì)算表明，此靜動(dòng)態(tài)“混合重構(gòu)”算法能夠明顯減少計(jì)算量和存儲(chǔ)量。上述靜態(tài)重構(gòu)具體策略可以采用類似k-exact的最小二乘重構(gòu)［14-15，25］、強(qiáng)插值重構(gòu)［16-17］、Gauss-Green公式重構(gòu)［18-22］等算法。

1.3 混合重構(gòu)及DG／FV空間離散

引入輔助變量Z，將式（1）重寫(xiě)為如下一階方程組：

單元e中，變量Uh和重構(gòu)變量Wh可以分別寫(xiě)成如下的形式：

其中cl和是多項(xiàng)式系數(shù)，bl（ξ）取參考單元內(nèi)的正交多項(xiàng)式基函數(shù)，ξ是參考單元的局部坐標(biāo)。式（4）中變量Uh即為需要計(jì)算和存儲(chǔ)的自由度，它的時(shí)間推進(jìn)過(guò)程即為“動(dòng)態(tài)重構(gòu)”。另外可以看到，重構(gòu)變量Wh可以分為兩部分，它的低階部分為了保證守恒性直接取為Uh，而高階部分即為需要通過(guò)網(wǎng)格模板“靜態(tài)重構(gòu)”的部分。

我們數(shù)值求解式（3）即是尋求Uh∈，滿足式（3）的弱形式。設(shè)由某種重構(gòu)算法得到重構(gòu)變量Wh∈（具體參考文獻(xiàn)［14-22］），并令輔助變量Zh∈）d。將式（3）乘以檢驗(yàn)函數(shù)bj，然后在單元Ωe上使用分部積分，可得其弱形式為：

上式中為了消除Wh、Fc和Fv在邊界?Ωe上的二義性，分別使用了其數(shù)值通量、Hc和Hv。和分別表示相鄰單元的重構(gòu)變量和輔助變量，n為邊界單位外法向。無(wú)粘數(shù)值通量Hc可以使用常規(guī)的歐拉黎曼求解器得到［18-22］，本文使用Roe格式［27］。與粘性相關(guān)的數(shù)值通量、Hv如何計(jì)算是DG／FV混合方法求解N-S方程的關(guān)鍵內(nèi)容。本文使用由Bassi和Rebay提出的BR格式［28-29］。由于BR1格式需要求解和存儲(chǔ)輔助變量Zh，其計(jì)算量和存儲(chǔ)量大，且BR1格式需要用到鄰居的鄰居單元信息，其不緊致，所以此處使用其改進(jìn)后的BR2格式。具體如下，由式（5）再經(jīng)過(guò)一次分部積分，得：

依此定義提升算子Rh＝Zh－▽W(xué)h，則在每個(gè)單元中Rh可由式（7）方便的算出（由于使用正交基函數(shù)，不用單元體積分），而▽W(xué)h可以直接得到，因此可得Zh＝▽W(xué)h＋Rh。進(jìn)一步定義修正提升算子rf為：

其中f是單元e的一個(gè)邊界面。

粘性數(shù)值通量Hv的計(jì)算公式為：

其中ηf為穩(wěn)定性因子，一般取3。

式（9）定義的關(guān)于面f的修正提升算子rf僅依賴本單元和面f對(duì)應(yīng)的相鄰單元，因此式（10）中面f的粘性通量只與此兩個(gè)單元相關(guān)，從而具有緊致性。

最終的半離散方程可以寫(xiě)為：

1.4 時(shí)間離散

隨著空間精度階的提高，使用顯式時(shí)間離散時(shí)DGM及DG／FV方法對(duì)應(yīng)的穩(wěn)定性條件將越來(lái)越嚴(yán)格，時(shí)間推進(jìn)步長(zhǎng)受到明顯的限制，而隱式時(shí)間離散方法的時(shí)間推進(jìn)步長(zhǎng)受計(jì)算精度階和計(jì)算網(wǎng)格尺度的限制較小，因此本文使用之前發(fā)展的基于Newton／Gauss-Seidel迭代的時(shí)間隱式離散方法［22］。以如下一階Euler后差時(shí)間離散為例：

定義關(guān)于Un＋1的非定常殘差Run（Un＋1）為：

若能求解Run（Un＋1）＝0則可得到滿足一階精度的Un＋1，但上式關(guān)于Un＋1是非線性的，不能直接求解。對(duì)其使用Newton迭代方法求解，可得：

其中l(wèi)是Newton迭代指標(biāo)，l＝0時(shí)Un＋1，l取Un，l→∞時(shí)Un＋1，l＋1＝Un＋1，l＝Un＋1。

式（14）是關(guān)于ΔUn＋1，l的大型線性方程組，但注意到RHS只與目標(biāo)單元以及若干鄰近的單元相關(guān)，從而式（14）左端矩陣是一個(gè)塊稀疏矩陣。本文采用Gauss-Seidel迭代來(lái)求解式（14），具體方法見(jiàn)文獻(xiàn)［21-22］。

以上求解過(guò)程可能還需要對(duì)曲邊界進(jìn)行修正，其它如邊界條件和限制策略［30］等限于篇幅在此不再贅述。

2 數(shù)值特性分析

2.1 數(shù)值色散耗散特性

流體力學(xué)方程是復(fù)雜的非線性方程組，一般來(lái)說(shuō)此方程組的數(shù)學(xué)性質(zhì)，如解的存在性、唯一性、數(shù)學(xué)提法的適定性等等都還是正在研究中的問(wèn)題，且很難找到一般情況下方程的解析解或精確解。所以這里我們僅針對(duì)一維線性對(duì)流方程分析DG／FV混合方法的色散、耗散特性，得到幾種典型精度階的DG／FV格式的數(shù)值特性，并和典型精度階DG格式的數(shù)值特性進(jìn)行比較。

圖1給出了幾種典型格式的數(shù)值色散關(guān)系曲線，此處及下文中各格式括號(hào)中的數(shù)字表示其設(shè)計(jì)精度階。從圖1可以看出，在高波數(shù)區(qū)，DG／FV（3）的色散特性接近DGM（2），在低波數(shù)區(qū)，DG／FV（3）的色散特性介于DGM（2）和DGM（3）之間，表現(xiàn)出令人滿意的色散特性。圖2給出了幾種典型格式的數(shù)值耗散關(guān)系曲線，可以看到其數(shù)值耗散均小于零，說(shuō)明格式是穩(wěn)定的。DG／FV（2）格式的耗散較大，DG／FV（3）的耗散特性同樣介于DGM（2）和DGM（3）之間，既表現(xiàn)出低波數(shù)區(qū)有盡量小的耗散，在高波數(shù)區(qū)又具有較大的耗散，耗散特性比較理想。

圖1 幾種典型格式的數(shù)值色散關(guān)系Fig.1 Dispersive behavior of some schemes

圖2 幾種典型格式的數(shù)值耗散關(guān)系Fig.2 Dissipative behavior of some schemes

圖3和圖4分別為DG／FV（3）格式的非物理模態(tài)的數(shù)值色散和耗散曲線。從圖3可以看出，非物理模態(tài)波的傳播方向與物理模態(tài)波相反，不過(guò)從圖4可以看出非物理模態(tài)波數(shù)值耗散很大，在短距離、短時(shí)間內(nèi)就會(huì)被耗散掉，這與Hu和Atkins［31］研究所得DGM的數(shù)值色散耗散特性是類似的。

圖3 DG／FV（3）格式的非物理模態(tài)數(shù)值色散關(guān)系Fig.3 Dispersive of DG／FV（3）in spurious mode

2.2 穩(wěn)定性分析

使用Eourier方法對(duì)模型方程進(jìn)行穩(wěn)定性分析，得到了一些典型階數(shù)DG和DG／FV格式的穩(wěn)定性條件（見(jiàn)表1），其中時(shí)間離散使用1-5階的顯式RK方法。由表1可以看出DG／FV格式的CEL數(shù)與低一階DG格式的CEL數(shù)相當(dāng)，比同階DG格式的CEL數(shù)有大幅提高，說(shuō)明了此DG／FV混合格式在時(shí)間推進(jìn)步長(zhǎng)方面具有顯著優(yōu)勢(shì)。

圖4 DG／FV（3）格式的非物理模態(tài)的數(shù)值耗散關(guān)系Fig.4 Dissipative of DG／FV（3）in spurious mode

表1 幾種DG與DG／FV格式穩(wěn)定性條件Table 1 CFL number of several DG and DG／FV schemes

2.3 計(jì)算量理論分析

文獻(xiàn)［32］在對(duì)數(shù)值通量的計(jì)算量作了合理的假設(shè)后給出了幾種精度DG格式的計(jì)算量的理論分析結(jié)果。本文參考其方法，理論分析了常用的3-4階精度DG和DG／FV混合格式的計(jì)算量，其中二維標(biāo)量方程三角形網(wǎng)格下所得結(jié)果見(jiàn)表2。從中我們可以看出：1）無(wú)論是DG還是DG／FV格式，單元高斯積分的計(jì)算量都占總計(jì)算量很大比重（40%～50%），這說(shuō)明了減少單元積分點(diǎn)數(shù)目的重要性；2）邊界高斯積分的計(jì)算量所占比重隨著方法精度階數(shù)提高逐漸減?。?）DG／FV格式中重構(gòu)運(yùn)算計(jì)算量占總計(jì)算量約15%，數(shù)值計(jì)算也表明重構(gòu)計(jì)算量比重較小；4）DG／FV格式計(jì)算量只有同階DG格式計(jì)算量的50%左右。

表2 幾種DG與DG／FV格式計(jì)算量Table 2 Costs of several DG and DG／FV schemes

表3給出了本文使用的幾種DG和DG／FV格式所用的單元高斯積分和邊界積分節(jié)點(diǎn)數(shù)。從中可以看出，相比同階精度的DG格式，DG／FV混合格式所用積分點(diǎn)數(shù)大幅減少。這正是DG／FV格式較同階精度DG格式計(jì)算量小的主要原因。

表4給出了二維標(biāo)量方程DG和DG／FV格式計(jì)算效率的數(shù)值計(jì)算結(jié)果。與三階DGM（3）相比，同階精度的DG／FV（3）混合方法計(jì)算量節(jié)省了約40%～50%，存儲(chǔ)量也節(jié)省了約30%～40%。

表3 幾種DG與DG／FV格式積分點(diǎn)數(shù)目Table 3 Quadrature points of DG and DG／FV schemes

表4 DG與DG／FV計(jì)算效率、存儲(chǔ)量比較（二維標(biāo)量方程）Table 4 Efficiency and storage of DG and DG／FV schemes

3 數(shù)值計(jì)算結(jié)果與分析

3.1 庫(kù)埃特流動(dòng)

為了驗(yàn)證DG／FV混合方法的數(shù)值精度階，首先計(jì)算了庫(kù)埃特（Couette）流動(dòng)問(wèn)題。此問(wèn)題描述的是在兩個(gè)平板y＝0和y＝H中，由上平板勻速運(yùn)動(dòng)所導(dǎo)致的槽道流動(dòng)。取上平板速度U＝1，溫度T1＝0．85；下平板溫度T0＝0．8；粘性系數(shù)μ＝0．01為常數(shù)；槽道高度H＝2。計(jì)算區(qū)域取為［0，4］×［0，2］，左右邊界使用周期邊界條件，上下邊界使用等溫壁條件。計(jì)算使用四套網(wǎng)格，粗網(wǎng)格為120個(gè)三角形單元，其余網(wǎng)格由粗網(wǎng)格逐次加倍而得。此問(wèn)題的精確解如下［8］：

表5給出了計(jì)算所得密度與其精確解的誤差的L2模及數(shù)值精度階。可以看出各階DG／FV混合格式的數(shù)值精度階數(shù)均達(dá)到了設(shè)計(jì)精度，同時(shí)DG／FV格式的計(jì)算誤差位于低一階DGM和同階DGM結(jié)果之間，且隨著網(wǎng)格加密更接近于同階DGM結(jié)果。

表5 Couette流動(dòng)問(wèn)題密度L2模誤差Table 5 Error of Couette flow problem in L2norm ofdensity

3.2 方腔流動(dòng)

方腔流動(dòng)是一個(gè)層流計(jì)算的經(jīng)典算例，本文采用的是邊長(zhǎng)為1的正方形空腔，計(jì)算使用混合網(wǎng)格，其中包含1440個(gè)四邊形單元和1110個(gè)三角形單元，如圖5所示。計(jì)算條件是Ma＝0．1，Re＝1×104。

圖5 計(jì)算使用混合網(wǎng)格，單元數(shù)為2550Fig.5 Hybrid grid with 2550 cells

計(jì)算所得流線和文獻(xiàn)［33］結(jié)果如圖6所示，其中文獻(xiàn)結(jié)果采用二階精度有限差分方法和257×257的均勻結(jié)構(gòu)網(wǎng)格。從流線圖看，在如此稀疏的網(wǎng)格上，本文4階DG／FV（4）格式和DGM（4）格式計(jì)算結(jié)果與文獻(xiàn)都符合的很好，方腔內(nèi)一個(gè)主渦在中心位置，在三個(gè)角附近分布了三個(gè)較小的二次渦。本文計(jì)算比較準(zhǔn)確地分辨出了三個(gè)二次渦以及右下角的小渦。圖7給出了計(jì)算域內(nèi)x＝0.5直線位置的速度u剖面分布和y＝0.5直線上的速度v剖面分布，并與文獻(xiàn)結(jié)果進(jìn)行了對(duì)比?？梢钥闯?，本文DG／FV（4）格式和DGM（4）格式計(jì)算結(jié)果均與文獻(xiàn)結(jié)果符合的很好。從圖7還可以看出，4階DG／FV（4）格式計(jì)算結(jié)果比3階DGM（3）格式計(jì)算結(jié)果有很大改進(jìn)。

圖6 幾種格式計(jì)算所得流線與文獻(xiàn)結(jié)果比較Fig.6 Streamlines of computation compared with result in ref.［33］

圖7 幾中格式計(jì)算所得速度剖面Fig.7 Velocity profile of several schemes 7

圖8給出了幾種典型格式的顯式和隱式時(shí)間離散的收斂歷程比較，可以看出各階格式采用隱式時(shí)間離散時(shí)計(jì)算CPU時(shí)間都大大減少，收斂性能提高約1～2個(gè)量級(jí)。另外可以看出，盡管迭代步數(shù)幾乎相同，相比同階精度DGM，DG／FV混合格式計(jì)算CPU時(shí)間顯著減少，計(jì)算效率明顯提高。

圖8 幾種典型格式殘差收斂曲線Fig.8 Convergence history of several schemes

3.3 亞聲速剪切層流動(dòng)

對(duì)于剪切層算例，計(jì)算條件為Ma1＝0．5和Ma2＝0．25，Re＝ρ1U1δ／μ＝500。計(jì)算域?yàn)椋?，800］×［-100，100］，使用25，272個(gè)單元的混合網(wǎng)格。計(jì)算初始條件和邊界條件參考Colonius等的文獻(xiàn)［34］，計(jì)算最終時(shí)間為t＝1357．28＝68Tf。圖9給出了兩種格式計(jì)算所得渦量等值線云圖并與文獻(xiàn)［34］結(jié)果進(jìn)行了比較?？梢钥闯?，本文4階DG／FV（4）和5階DG／FV（5）計(jì)算結(jié)果很一致并和文獻(xiàn)結(jié)果類似。

圖9 剪切層算例渦量等值線云圖與文獻(xiàn)結(jié)果比較Fig.9 Vorticity contours at time t＝68 Tf

3.4 低速平板湍流流動(dòng)

對(duì)于低速平板湍流算例，計(jì)算條件為來(lái)流Ma∞＝0．2，Re＝1．03×107，計(jì)算網(wǎng)格為76×61的結(jié)構(gòu)網(wǎng)格，x方向平板端點(diǎn)附近第一層網(wǎng)格尺度為2×10-3，y方向第一層網(wǎng)格尺度為2×10-6。湍流模型為前述的SA模型。圖10給出了湍流平板某個(gè)站位下速度剖面分布的三階DG／FV（3）格式計(jì)算結(jié)果，并和二階FVM結(jié)果進(jìn)行了比較，其中二階FVM計(jì)算使用網(wǎng)格單元數(shù)量4倍于三階格式計(jì)算所用網(wǎng)格。圖10可以看出，較粗網(wǎng)格下的三階計(jì)算結(jié)果與密網(wǎng)格下的二階FV計(jì)算結(jié)果相符良好。

3.5 NACA0012翼型湍流流動(dòng)

圖10 湍流平板計(jì)算速度剖面分布Fig.10 Velocity profile of turbulent flat plate

對(duì)于NACA0012翼型繞流算例，計(jì)算條件為來(lái)流Ma∞＝0.7，α＝1.49°，Re＝9×106。湍流模型仍是前述的SA模型。計(jì)算使用網(wǎng)格分布為前緣0.0015，尾緣0.003，法向4.75×10-6，遠(yuǎn)場(chǎng)15倍弦長(zhǎng)，網(wǎng)格量為265×132，如圖11所示。計(jì)算所得升力、阻力系數(shù)同文獻(xiàn)［35］結(jié)果比較見(jiàn)表6，從表中可以看出，3階DG／FV（3）計(jì)算結(jié)果介于DGM（2）和DGM（3）之間，和文獻(xiàn)結(jié)果也比較符合。

表6 NACA0012翼型升阻力系數(shù)Table 6 Lift coefficient anddrag coefficient of NACA0012

圖11 NACA0012翼型計(jì)算所用網(wǎng)格Fig.11 Grid of NACA 0012 airfoil

圖12 NACA0012翼型壓力等值線分布（DG／FV（3））Fig.12 Pressure contours of NACA0012 airfoil（DG／FV（3））

圖13 NACA0012翼型計(jì)算所得湍流粘性系數(shù)分布（DG／FV（3））Fig.13 Turbulent viscosity contours of NACA0012 airfoil（DG／FV（3））

4 結(jié) 論

本文對(duì)基于非結(jié)構(gòu)／混合網(wǎng)格的高階精度DG／FV混合格式進(jìn)行了簡(jiǎn)要的綜述，重點(diǎn)介紹了“混合重構(gòu)”基本思想、與BR2方法結(jié)合的RANS方程求解方法、Newton／Gauss-Seidel耦合隱式計(jì)算格式，理論分析了其色散耗散特性及穩(wěn)定性條件，并給出DG／FV混合方法與DG方法計(jì)算量的定性理論分析和數(shù)值結(jié)果，表明相比同階精度的DG方法，DG／FV混合方法計(jì)算量減少約40%。隨后，利用多個(gè)典型算例考察了DG／FV的計(jì)算精度和計(jì)算效率，并與相應(yīng)的DG方法進(jìn)行了對(duì)比。對(duì)有解析解的Couette流動(dòng)問(wèn)題計(jì)算表明，本文方法達(dá)到了設(shè)計(jì)精度；方腔流動(dòng)、剪切層算例和低速平板湍流、NACA0012翼型繞流等經(jīng)典算例結(jié)果表明，隨著精度的提高，在較粗的計(jì)算網(wǎng)格上亦能得到高精度的計(jì)算結(jié)果，其計(jì)算結(jié)果的精度與同階的DG方法相當(dāng)，但是計(jì)算效率大幅提高；而隱式方法進(jìn)一步提高了本文方法的計(jì)算效率。由此可見(jiàn)高階精度DG／FV混合格式將具有良好的工程應(yīng)用前景。

下一步我們將針對(duì)三維復(fù)雜外形的湍流數(shù)值模擬進(jìn)行研究，重點(diǎn)研究曲邊界修正方法及高精度邊界條件的實(shí)現(xiàn)、針對(duì)高速流動(dòng)的限制器和間斷偵測(cè)方法或人工粘性方法、基于GMRES的隱式計(jì)算方法、三維大規(guī)模分區(qū)并行計(jì)算技術(shù)、湍流模型的緊耦合高精度計(jì)算方法、基于動(dòng)網(wǎng)格的高精度非定常計(jì)算方法及幾何守恒問(wèn)題、與DES和LES方法的有機(jī)結(jié)合等，期待早日實(shí)現(xiàn)在復(fù)雜工程問(wèn)題中的成功應(yīng)用。

［1］BAKER T J．Mesh generation：art or science［J］．Progress inAerospace Sciences，2005，41：29-63．

［2］WANG Z J．High-order methods for the Euler and Navier-Stokes equations on unstructured grids［J］．Progress in Aerospace Sciences，2007，43：1-41．

［3］BARTH T J，EREDERICKSON P O．Higher order solution of the Euler equations on unstructured grids using quadratic reconstruction［R］．AIAA Paper 90-0013．

［4］OLLIVIER-GOOCH C，NEJAT A，MICHALAK K．Obtaining and verifying high-order unstructured finite volume solutions to the Euler equations［J］．AIAA Journal，2009，47（9）：2105-2120．

［5］COCKBURN B，SHU C W．The Runge-Kuttadiscontinuous Galerkin method for conservation laws V：multidimensional systems［J］．J．Comput．Phys．，1998，141：199-224．

［6］HE L X，ZH ANG L P，ZHANG H X．Discontinuous Galerkin finite element method on 3D arbitrary elements［J］．ACTA Aerod ynamica Sinica，2007，25（2）：157-162．（in Chinese）

賀立新，張來(lái)平，張涵信．任意單元間斷Galerkin有限元計(jì)算方法研究［J］．空氣動(dòng)力學(xué)學(xué)報(bào)，2007，25（2）：157-162．

［7］SUN Y，WANG Z J，LIU Y．Spectral（finite）volume method for conservation laws on unstructured grids VI：extension to viscous flow［J］．J．Comput．Phys．，2006，215：41-58．

［8］SUN Y Z，WANG Z J，LIU Y．High-order multidomain spectraldifference method for the Navier-Stokes Equations on unstructured hexahedral grids［J］．Com mun．Comput．Phys．，2007，2（2）：310-333．

［9］WANG Z J，GAO H，HAGA T．A unifyingdiscontinuous formulation for hybrid meshes［M］．In：Adaptive High-Order Methods in Computational Eluid Dynamics．Edited by WANG Z J．Singapore：World Scientific Publishing，2011．

［10］EKATERINARIS J A．High-order accurate，low numericaldiffusion methods for aerodynamics［J］．Prog．Aero．Sci．，2005，41：192-300．

［11］ZHANG L P，HE L X，LIU W，et al．Reviews of high-order methods on unstructured and hybrid grid［J］．Advances in Mechanics，2013，43（2）：202-236．（in Chinese）

張來(lái)平，賀立新，劉偉，等．基于非結(jié)構(gòu)／混合網(wǎng)格的高階精度格式研究進(jìn)展［J］．力學(xué)進(jìn)展，2013，43（2）：202-236．

［12］COCKBURN B，LUSKIN M，SHU C W，SULI E．Enhanced accuracy by post-processing for finite element methods for hyperbolic equations［J］．Mathematics of Computation，2003，72：577-606．

［13］RYAN J K，SHU C W，ATKINS H L．Extension of a post-processing technique for thediscontinuous Galerkin method for hyperbolic equations with applications to an aeroacoustic problem［J］．SIAM Journal on Scientific Computing，2005，26：821-843．

［14］DUMBSER M，BALSARA D S，TORO E E．A unified framework for the construction of one-step finite volume anddiscontinuous Galerkin schemes on unstructured meshes［J］．J．Comput．Phys．，2008，227：8209-8253．

［15］DUMBSER M．Arbitrary high orderPNPMschemes on unstructured meshes for the compressible Navier-Stokes equations［J］．Computers and Eluids，2010，39：60-76．

［16］LUO H，LUO L P，NOURGALIEV R，et al．A reconstructeddiscontinuous Galerkin method for the compressible Navier-Stokes equations on arbitrary grids［J］．J．Comput．Phys．，2010，229：6961-6978．

［17］LUO H，LUO L P，ALI A，et al．A parallel，reconstructeddiscontinuous Galerkin method for the compressible flows on arbitrary grids［J］．Com mun．Comput．Phys．，2011，9（2）：363-389．

［18］ZHANG L P，LIU W，HE L X，et al．A class of hybrid DG／FV methods for conservation laws I：basic formulation and onedimensional systems［J］．J．Comput．Phys．，2012，231：1081-1103．

［19］ZHANG L P，LIU W，HE L X，et al．A class of hybrid DG／FV methods for conservation laws II：two-dimensional cases［J］．J．Comput．Phys．，2012，231：1104-1120．

［20］ZHANG L P，LIU W，HE L X，et al．A class of hybrid DG／FV methods for conservation laws III：two-dimensional Euler equations［J］．Commun．Comput．Phys．，2012，12（1）：284-314．

［21］ZHANG L P，LIU W，LI M，et al．A class of DG／FV hybrid schemes for conservation law IV：2D viscous flows and implicit algorithm［J］，Computers and Eluids，2014，97：110-125．

［22］ZHANG L P，LI M，LIU W，et al．An implicit algorithm of high-order DG method and hybrid DG／FV schemes based on Newton／Gauss-Seidel iteration for 2D inviscid flows on arbitrary grids［J］．Com munication in Computational Physics，2014（in press）．

［23］WANG Z J，SHI L，EU S，et al．APNPM-CPR framework for hyperbolic conservation laws［R］．AIAA Paper 2011-3227．

［24］SHI L，WANG Z J，EU S，et al．APNPM-CPR method for Navier-Stokes equations［R］．AIAA Paper 2012-460．

［25］LI M，LIU W，ZHANG L P，et al．A class of high order hybrid DG／FV methods for 2D viscous flows［J］．Transactions of Nanjing University of Aeronautics and Astronautics，2013，30（S）：68-73．

［26］PETERSON W，DAVID W Z．Three-dimensional aerodynamic computations on unstructured grids using a Newton-Krylov approach［J］．Computers and Eluids，2008，37（2）：107-120．

［27］ROE P L．Approximate Riemann solvers，parameter vectors，anddifference schemes［J］．J．Comput．Phys．，1981，43：357-372．

［28］BASSI E，REBAY S．A high-order accuratediscontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations［J］．J．Comput．Phys．，1997，131：267-293．

［29］BASSI E，REBAY S，MARIOTTI G，et al．A high-order accuratediscontinuous finite element method for inviscid and viscous turbomachinery flows［A］．In：2nd European Conference on Turbomachinery Eluid Dynamics and Thermodynamics［C］．Decuypere R，Dibelius G，editors．Technologisch Instituut，Antwerpen，Belgium，1997：99-108．

［30］ZHANG L P，LIU W，HE L X，et al．A shockdetection method and applications in DGM for hyperbolic conservation laws on unstructured grids［J］．ACTA Aerod ynamica Sinica，2011，29（4）：401-406．（in Chinese）

張來(lái)平，劉偉，賀立新，等．一種新的間斷偵測(cè)器及其在DGM中的應(yīng)用［J］．空氣動(dòng)力學(xué)學(xué)報(bào)，2011，29（4）：401-406．

［31］HU E，ATKINS H．Eigensolution analysis of thediscontinuous Galerkin method with nonuniform grids，I：one spacedimension［J］．J．Comput．Phys．，2002，182（2）：516-545．

［32］SUN Y Z，WANG Z J．Evaluation ofdiscontinuous Galerkin and spectral volume methods for scalar and system conservation laws on unstructured grids［J］．Int．J．Numer．Meth．Eluids，2004，45：819-838．

［33］GHIA U，GHIA K N，SHIN C T．HighResolutions for incompressible flow using the Navier-Stokes equations and a multigrid method［J］．J．Comput．Phys．，1982，48：387-411．

［34］COLONIUS T，LELE S K，MOIN P．Sound generation in a mixing layer［J］．J．Eluid Mech．，1997，330：375-409．

［35］MAKSYMIUK C M，PULLIAM T H．Viscous transonic airfoil workshop results using ARC2D［R］．AIAA Paper 87-0415．

Recent development of high order DG／FV hybrid methods

ZHANG Laiping1，2，LI Ming2，LIU Wei1，2，HE Xin1，2，ZHANG Hanxin2
（1．State Key Laboratory of Aerodynamics，China Aerod ynamics Research and Development Center，Mianyang，Sichuan 621000，China；2．Computational Aerodynamics Institute，China Aerod ynamics Research and Development Center，Mianyang，Sichuan 621000，China）

A concept of‘static reconstruction’and‘dynamic reconstruction’had been introduced for higher-order（third-order and higher）numerical methods in our previous work．Based on this concept，a class of DG／FV hybrid methods had beendeveloped for the scalar equations and Euler／NS equations on triangular and Cartesian／triangular hybrid grids．In this paper，the recent progress of the DG／FV hybrid methods was presented．The basic idea of‘hybrid reconstruction’，the procedure of solving NS equations with BR2 approach，and the implicit algorithm were reviewed briefly．And then thedissipative anddispersive property，as well as the stability，of the DG／FV hybrid schemes were analyzed．In order to show the high efficiency in the term of CPU time of the present DG／FV hybrid schemes，the computational costs werediscussed and compared with the corresponding DG methods．The numerical accuracy was validated by some typical test cases of viscous flow，including the Couette flow，laminar flow in a square，compressible mixing layer problem，turbulent flows by RANS equations with S-A turbulent model over a flat plate and over NACA0012 airfoil．The accuracy study shows that the hybrid DG／FV method achieves thedesired order of accuracy，and they can capture the flow structure accurately．Qualitative analysis and numerical applicationsdemonstrate that they can reduce the CPU time greatly（up to 40%）comparing with the traditional DG method with the same order of accuracy．Meanwhile，the implicit algorithm can accelerate the convergence history obviously，one to two orders faster than the explicit Runge-Kutta method．

unstructured／hybrid grid；discontinuous Galerkin method；finite volume method；DG／FV hybrid method；RANS equations

V211.3

Adoi：10.7638／kqdlxxb-2014.0123

0258-1825（2014）06-0717-10

2014-10-16；

2014-12-02

國(guó)家自然科學(xué)基金（91130029，11402290）

張來(lái)平（1968-），研究員，博士生導(dǎo)師，主要從事非結(jié)構(gòu)／混合網(wǎng)格生成技術(shù)、基于非結(jié)構(gòu)／混合網(wǎng)格的計(jì)算格式、非定常流動(dòng)機(jī)理等方面的研究與應(yīng)用工作．E-mail:zhanglp＿cardc＠126.com

張來(lái)平，李明，劉偉，等．基于非結(jié)構(gòu)／混合網(wǎng)格的高階精度DG／FV混合方法研究進(jìn)展［J］．空氣動(dòng)力學(xué)學(xué)報(bào)，2014，32（6）：717-726．

10．7638／kqdlxxb-2014．0123 ZHANG L P，LI M，LIU W，et al．Recentdevelopment of high order DG／FV hybrid methods［J］．ACTA Aerodynamica Sinica，2014，32（6）：717-726．