可壓縮流動模擬的多重網(wǎng)格-擾動域推進(jìn)方法

胡姝瑤，蔣崇文，李椿萱

北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院 國家計算流體力學(xué)實(shí)驗(yàn)室，北京 100191

多重網(wǎng)格方法是一種在粗、細(xì)網(wǎng)格間循環(huán)求解控制方程的加速收斂技術(shù)，是飛行器繞流數(shù)值模擬加速中廣泛使用的手段之一［1］。多重網(wǎng)格法的加速原理有2 方面：一是粗網(wǎng)格以更大的時間步長推進(jìn)求解，從而帶動細(xì)網(wǎng)格更快達(dá)到收斂態(tài)；二是粗、細(xì)網(wǎng)格間的循環(huán)求解有助于衰減阻礙細(xì)網(wǎng)格收斂的高頻誤差，從而使細(xì)網(wǎng)格的整體誤差能更快地衰減［2］。

多重網(wǎng)格方法起源于橢圓型方程的數(shù)值計算，并已被證明對于橢圓型問題是一種最優(yōu)化的數(shù)值方法［3］。Ni［4］最先將多重網(wǎng)格方法推廣至雙曲型的非定常Euler 方程，建立了定?？蓧嚎s流動的顯式多重網(wǎng)格求解器。Jameson 等［5］將隱式時間格式引入多重網(wǎng)格方法，進(jìn)一步提高了Euler 方程的求解效率。Yadlin 等［6］將隱式多重網(wǎng)格方法推廣至多塊網(wǎng)格，驗(yàn)證了多重網(wǎng)格內(nèi)邊界的處理方式。隨后，該方法又被推廣至黏性［7］、高超聲速［8］、化學(xué)反應(yīng)［9］等復(fù)雜流動的模擬，應(yīng)用于螺旋槳滑流［1］、外掛物投放［10］、壓氣機(jī)內(nèi)流［11］、噴流干擾［12］等工程問題的研究中。

當(dāng)前，多重網(wǎng)格方法［13］已成功應(yīng)用于各類流動的研究中，可實(shí)現(xiàn)數(shù)十倍的收斂速率增長。不過，無論是在細(xì)網(wǎng)格還是粗網(wǎng)格上，現(xiàn)有多重網(wǎng)格方法在每個迭代步內(nèi)均必須對所有網(wǎng)格單元更新求解。這一全局更新的求解方式雖便于實(shí)現(xiàn)，卻未考慮流場擾動傳播與數(shù)值迭代收斂的特征，可能導(dǎo)致大量無效計算，限制計算效率。

針對傳統(tǒng)全局更新求解方式存在無效計算的問題，本文作者前期工作［14-17］提出了擾動域推進(jìn)方法，以提升可壓縮流動數(shù)值模擬的計算效率。擾動域推進(jìn)方法的加速原理是采用隨擾動傳播而增大、隨求解收斂而縮小的動態(tài)計算域，實(shí)現(xiàn)了僅在未收斂的受擾區(qū)域求解，N-S 方程、Euler 方程動態(tài)分區(qū)耦合求解的局部更新求解方式。通過有效避免傳統(tǒng)方法中的無效計算，擾動域推進(jìn)方法可同時顯著降低可壓縮流動數(shù)值模擬的單步計算量和存儲需求。對于典型亞聲速流動計算時間節(jié)省可達(dá)54.0%，對于典型超聲速流動時間節(jié)省可達(dá)74.5%［17］。

本文將在擾動域推進(jìn)方法的計算框架下融入多重網(wǎng)格方法，著重建立粗網(wǎng)格上動態(tài)計算域的演化邏輯，設(shè)計對結(jié)構(gòu)網(wǎng)格更為普適的新型粗網(wǎng)格生成策略，提出多重網(wǎng)格-擾動域推進(jìn)方法，并通過典型算例驗(yàn)證多重網(wǎng)格、擾動域推進(jìn)2 類加速技術(shù)的耦合加速效果。

1 控制方程及數(shù)值方法

1.1 控制方程

擾動域推進(jìn)方法采用了有黏、無黏流動的動態(tài)分區(qū)計算，其中模擬無黏、層流和湍流的控制方程分別為Euler 方程、N-S 方程和RANS 方程。若忽略體積力與熱源，3 類方程可統(tǒng)一描述為

式中：t為時間；Ω表示控制體，?Ω為其邊界；dΩ、dS分別為體積分及面積分的微元。在笛卡爾坐標(biāo)系（x1，x2，x3）下，守恒變量W、對流通量Fc、黏性通量Fv以及源項(xiàng)Q可表示為

式中：ρ、P、T、E和H分別表示密度、壓強(qiáng)、溫度、比總能和比總焓；λ為熱傳導(dǎo)系數(shù)；μ為動力黏性系數(shù)；ui為速度分量；τij為分子剪應(yīng)力張量；δij為Kronecker 函數(shù)（i，j=1，2，3）；下標(biāo)“T”表示湍流模型。

對于無黏流動，去除式（1）中包含μ以及湍流相關(guān)的項(xiàng)，包括Fv、Q、WT、Fc，T，式（1）即可退化為Euler 方程。對于層流流動，去除湍流相關(guān)項(xiàng)，μ設(shè)為滿足Sutherland 公式的分子黏性系數(shù)μL，式（1）即為N-S 方程。對于湍流流動，保留所有項(xiàng)，μ設(shè)為分子黏性與湍流渦黏性系數(shù)之和，即μ=μL+μT，式（1）即為滿足Boussinesq 渦黏性假設(shè)、省略了Reynolds 平均與Favre 平均符號的RANS 方程。

1.2 定?？蓧嚎s流動的有限體積方法

本文方法將基于結(jié)構(gòu)貼體網(wǎng)格的有限體積法實(shí)現(xiàn)。假設(shè)式（1）中，W、Q在一個網(wǎng)格單元內(nèi)恒定，通量的曲面積分可由網(wǎng)格單元面的通量加和近似，則式（1）可離散為

式中：殘差R可表示為

其中：|Ω|表示網(wǎng)格單元的體積；Nf為網(wǎng)格單元的面數(shù)；nm和ΔSm分別為網(wǎng)格單元第m個面的單位外法向和面積（1≤m≤Nf）。

對于定常可壓縮流動，數(shù)值模擬旨在通過迭代求解式（3），確定出滿足dW/dt=0的W。其中，近似式（3）等號右端項(xiàng)與左端項(xiàng)的步驟分別稱為殘差估計和時間積分。殘差估計以上一迭代步的守恒量W（n?1）為輸入，根據(jù)式（4），利用空間離散格式近似當(dāng)前迭代步的殘差R（n）。時間積分則以R（n）為輸入，利用時間離散格式估計當(dāng)前迭代步的守恒量更新量ΔW（n），并更新守恒量，即W（n）=W（n?1）+ΔW（n）。在每一迭代步通過所有單元的||ΔW（n）||2判斷收斂情況，若所有單元的最大值小于給定的收斂閾值εc，即||ΔW（n）||2，max<εc，則停止計算。

傳統(tǒng)方法中，殘差估計與時間積分要在預(yù)設(shè)的計算域內(nèi)對所有單元執(zhí)行，兩者的耗時約占總計算時間的99%［14］，十分耗時。擾動域推進(jìn)方法正是旨在減少這2 個步驟的計算量。為避免混淆，下文中將同時求解所有單元的方法稱為傳統(tǒng)全局更新方法，其所使用的靜態(tài)計算域稱為預(yù)設(shè)計算域。

1.3 多重網(wǎng)格方法

多重網(wǎng)格方法需在一層細(xì)網(wǎng)格和若干層粗網(wǎng)格上，按指定規(guī)律循環(huán)求解式（3）。對于結(jié)構(gòu)網(wǎng)格，常用的粗網(wǎng)格生成方法是將細(xì)網(wǎng)格沿網(wǎng)格方向稀疏一倍［2］。常用的多重網(wǎng)格循環(huán)有V 型、W 型等，其均由基礎(chǔ)循環(huán)組成［18］?；A(chǔ)循環(huán)包括以下3 步：

1）將細(xì)網(wǎng)格的守恒量與殘差插值到粗網(wǎng)格。

在細(xì)網(wǎng)格求解更新后，將粗網(wǎng)格所包含細(xì)網(wǎng)格的體積加權(quán)平均守恒量作為粗網(wǎng)格的守恒量初值，即

式中：Nc為粗網(wǎng)格所含細(xì)網(wǎng)格的個數(shù)；下標(biāo)h、2h分別表示細(xì)網(wǎng)格和其稀疏一倍所得粗網(wǎng)格。粗網(wǎng)格的體積|Ω|2h可表示為

2）在粗網(wǎng)格上迭代求解控制方程。

由于粗網(wǎng)格的數(shù)值精度不會影響細(xì)網(wǎng)格，為提高效率，多重網(wǎng)格法通常只采用一階迎風(fēng)格式求解粗網(wǎng)格上的對流通量，也不在粗網(wǎng)格上求解湍流模型方程。為將細(xì)網(wǎng)格精度傳給粗網(wǎng)格，以源項(xiàng)形式在控制方程中引入考慮粗、細(xì)網(wǎng)格間精度差異的強(qiáng)迫函數(shù)（QF）2h，可定義為

粗網(wǎng)格所采用的時間推進(jìn)格式與細(xì)網(wǎng)格保持一致。粗網(wǎng)格上的迭代步數(shù)可任意給定，文獻(xiàn)［13］指出當(dāng)以LU-SGS 格式推進(jìn)求解時在粗網(wǎng)格上迭代2 步可使效率最佳。

3）將粗網(wǎng)格的守恒量修正量插值回細(xì)網(wǎng)格。

定義粗網(wǎng)格第n次迭代所得守恒量與粗網(wǎng)格的守恒量初值之差為粗網(wǎng)格的守恒量修正量，即

則細(xì)網(wǎng)格的新解可表示為

當(dāng)粗網(wǎng)格多于1 層且采用迎風(fēng)格式求解對流通量時，文獻(xiàn)［18］推薦粗網(wǎng)格通過式（11）獲得新解后，先迭代一步再將其修正量插值給更密網(wǎng)格以獲得最佳的多重網(wǎng)格性能。

1.4 擾動域推進(jìn)方法

本文作者的前期工作表明，定常可壓縮流動的求解過程具有4 點(diǎn)特征［15］：①流場的數(shù)值擾動始于定?？刂品匠涛礉M足之處；②擾動向周圍逐漸傳播，尚未受擾區(qū)域仍保持初始值；③整個預(yù)設(shè)計算域并非同時收斂，位于上游及遠(yuǎn)離壁面的區(qū)域會先收斂；④無論流動是否分離，黏性效應(yīng)僅在有限區(qū)域起主導(dǎo)作用。因此，全局更新方法在未受擾動和已收斂區(qū)域內(nèi)的求解對數(shù)值迭代不僅沒有積極貢獻(xiàn)，反而可能引入更多數(shù)值誤差，延遲迭代收斂。

為利用這4 點(diǎn)特征提高計算效率，本文提出了采用動態(tài)計算域的擾動域推進(jìn)方法。該方法的特征如下：①定義了對流、黏性2 類動態(tài)計算域，且黏性動態(tài)域?對流動態(tài)域?預(yù)設(shè)計算域；②初始動態(tài)計算域包括壁面邊界毗鄰單元或不滿足來流條件的單元；③動態(tài)計算域隨數(shù)值擾動傳播而增大，隨求解收斂而縮小；④殘差估計與時間積分僅在對流動態(tài)域中執(zhí)行，并僅在黏性動態(tài)域中考慮式（4）中的黏性項(xiàng)。

圖1 將擾動域推進(jìn)方法與全局更新方法的求解流程進(jìn)行了對比。圖中白色框?yàn)? 種方法的共有步驟，黃色框則為擾動域推進(jìn)方法的特有步驟。對于全局更新方法，白色框步驟將全在預(yù)設(shè)計算域中執(zhí)行。對于擾動域推進(jìn)方法，紅色虛線框中的步驟則在對流動態(tài)域中執(zhí)行，而綠色虛線框的步驟則僅在黏性動態(tài)域中執(zhí)行。

圖1 擾動域推進(jìn)方法與全局更新方法的流程對比Fig.1 Comparison of flow charts between Disturbance Region Update Method（DRUM）and Global-Update Method（GUM）

擾動域推進(jìn)方法的特有步驟可分為3 類：①迭代求解前，建立動態(tài)計算域；②在每一個迭代步中，更新動態(tài)計算域，依次執(zhí)行增大對流動態(tài)域、縮小對流動態(tài)域、增大黏性動態(tài)域和縮小黏性動態(tài)域4 個子步；③根據(jù)對流動態(tài)域，調(diào)整存儲空間。前2 類特有步驟旨在確保擾動域推進(jìn)方法能夠消除全局更新方法中的無效計算，第3 類則使得擾動域推進(jìn)方法可在提高計算效率的同時減少存儲空間占用。為便于閱讀，下文將簡要介紹動態(tài)計算域建立與更新算法的原理，算法的數(shù)學(xué)表達(dá)如表1 所示。算法的具體推導(dǎo)、經(jīng)驗(yàn)參數(shù)取值的討論詳見文獻(xiàn)［14-17］。

表1 擾動域推進(jìn)方法動態(tài)計算域更新算法Table 1 Update algorithms of dynamic computational domains in DRUM

1.4.1 建立動態(tài)計算域

根據(jù)流場初始化方式的不同，對流、黏性動態(tài)計算域的建立分為2 種情況：①當(dāng)根據(jù)來流條件初始化時，對流、黏性2 類動態(tài)域依據(jù)壁面邊界建立，分別取緊鄰壁面的10 層和1 層單元作為初始單元；②當(dāng)根據(jù)給定流場初始化時，對流動態(tài)域包含與來流條件不符的單元，黏性動態(tài)域包含對流動態(tài)域中的黏性主導(dǎo)單元。

1.4.2 增大對流動態(tài)域

遍歷對流動態(tài)域的邊界單元，對每一邊界單元執(zhí)行如下2 個子步：①判斷是否為受擾單元；②若是受擾單元，將其可能受到擾動的相鄰單元加入對流動態(tài)域；亞聲速流動中，可能受到擾動的是所有相鄰單元；而超聲速流動中，則是位于下游的單元。

1.4.3 縮小對流動態(tài)域

遍歷對流動態(tài)域的邊界單元，對每一邊界單元判斷如下4 個條件：①求解已收斂且周圍不存在新增受擾單元；②位于對流動態(tài)域的最上游；③對于亞聲速流動，不再影響對流動態(tài)域中其他單元的求解；④不再受對流動態(tài)域中其他單元的影響。若上述4 個條件均滿足，則將該單元從對流動態(tài)域中移除，該單元的緊鄰單元將成為對流動態(tài)域新的邊界單元，繼續(xù)判斷。由于黏性動態(tài)域?對流動態(tài)域，因此若該單元同時位于黏性動態(tài)域中，將其也從黏性動態(tài)域中移除。

對于超聲速流動，由于流動速度明顯大于聲速，數(shù)值求解過程具有“上游先于下游收斂”的規(guī)律。因此，對流動態(tài)域的更新須保持“從上游向下游縮小”，即嚴(yán)格滿足條件②。而對于亞聲速可壓縮流動，由于流動速度與聲速量級接近并且預(yù)設(shè)計算域較大，數(shù)值求解過程呈現(xiàn)“遠(yuǎn)場先于近壁區(qū)收斂”的規(guī)律。因此，條件②被適當(dāng)放寬以提高計算效率，但增加了條件③以避免降低計算精度。

1.4.4 增大黏性動態(tài)域

遍歷黏性動態(tài)域的邊界單元，對每一邊界單元執(zhí)行如下2 個子步：①判斷是否受黏性效應(yīng)主導(dǎo)；②若受黏性效應(yīng)主導(dǎo)，則將其所有緊鄰單元加入黏性動態(tài)域。

1.4.5 縮小黏性動態(tài)域

遍歷黏性動態(tài)域的邊界單元，對每一邊界單元判斷如下2 個條件：①位于黏性動態(tài)域的最上游；②不再受黏性效應(yīng)主導(dǎo)。若上述2 個條件均滿足，則將該單元從黏性動態(tài)域中移除，該單元的緊鄰單元將成為黏性動態(tài)域新的邊界單元，繼續(xù)判斷。

2 多重網(wǎng)格-擾動域推進(jìn)方法

為整合多重網(wǎng)格、擾動域推進(jìn)2 類方法的加速優(yōu)勢，本文提出多重網(wǎng)格-擾動域推進(jìn)方法。其求解流程如圖2 所示，有別于擾動域推進(jìn)方法的步驟由綠色框標(biāo)出。由圖可知，無論是在細(xì)網(wǎng)格還是粗網(wǎng)格上，擾動域推進(jìn)方法依舊在對流動態(tài)域內(nèi)求解控制方程，僅在黏性動態(tài)域內(nèi)考慮黏性效應(yīng)。動態(tài)計算域更新算法的運(yùn)算僅在細(xì)網(wǎng)格上執(zhí)行，而粗網(wǎng)格的動態(tài)計算域則根據(jù)粗、細(xì)網(wǎng)格間的對應(yīng)關(guān)系生成。本節(jié)將對多重網(wǎng)格-擾動域推進(jìn)方法的特有步驟進(jìn)行介紹，包括粗網(wǎng)格初始化、建立動態(tài)計算域、增大對流動態(tài)域、調(diào)整存儲空間、動態(tài)計算域中粗、細(xì)網(wǎng)格間的插值等。未提及的步驟與擾動域推進(jìn)方法完全相同。

圖2 多重網(wǎng)格-擾動域推進(jìn)方法流程圖Fig.2 Flow chart of DRUM with multigrid

2.1 粗網(wǎng)格初始化

在讀入細(xì)網(wǎng)格的網(wǎng)格格點(diǎn)坐標(biāo)、邊界條件和計算參數(shù)后，需首先對粗、細(xì)網(wǎng)格的幾何參數(shù)與流場參數(shù)進(jìn)行初始化，為后續(xù)控制方程的迭代計算提供參數(shù)與初值。粗網(wǎng)格初始化的關(guān)鍵有2 點(diǎn)：一是計算粗網(wǎng)格幾何參數(shù)與流場參數(shù)的方法；二是粗網(wǎng)格與細(xì)網(wǎng)格間的對應(yīng)關(guān)系。

數(shù)值實(shí)驗(yàn)表明，由于存在強(qiáng)迫函數(shù)，無論是依據(jù)粗網(wǎng)格的頂點(diǎn)坐標(biāo)精確計算其體積、面向量、格心等幾何參數(shù)，還是通過細(xì)網(wǎng)格的幾何參數(shù)近似，均對多重網(wǎng)格的加速效果沒有影響。本文采用通過細(xì)網(wǎng)格估計粗網(wǎng)格幾何參數(shù)的方式，即采用式（6）計算粗網(wǎng)格單元體積，格心坐標(biāo)、面向量均采用細(xì)網(wǎng)格單元值的加權(quán)平均值。

在迭代求解中，由于僅在對流動態(tài)域中進(jìn)行多重網(wǎng)格循環(huán)，因此多重網(wǎng)格-擾動域推進(jìn)法在流場初始化時就必須采用式（5）為粗網(wǎng)格的所有單元賦守恒量初值。

對于一個結(jié)構(gòu)網(wǎng)格，文獻(xiàn)［2］采用沿一網(wǎng)格方向?qū)⒁粚o鄰單元合二為一的方式生成粗網(wǎng)格，即對于一個單元數(shù)為2a×2b×2c的網(wǎng)格塊，第1 層粗網(wǎng)格的單元數(shù)為2a?1×2b?1×2c?1。若目標(biāo)粗網(wǎng)格層數(shù)是L，則網(wǎng)格單元數(shù)需滿足a≥L，b≥L，c≥L；若不滿足，則會出現(xiàn)單個無法合并單元，影響多重網(wǎng)格方法的使用和效果。

為消除多重網(wǎng)格使用對單元數(shù)的限制，本文提出一種新型粗網(wǎng)格生成策略。令m（0）表示讀入的網(wǎng)格在某一網(wǎng)格方向上的單元數(shù)，則第i層粗網(wǎng)格在網(wǎng)格方向的單元數(shù)可表示為

式中：mod（?）代表取模。第i層網(wǎng)格的單元I（i）對應(yīng)于第i+1 層網(wǎng)格單元的標(biāo)號為

圖3 以L=2、m（0）=11 的情況為例，示意了本文所提粗網(wǎng)格生成策略的特點(diǎn)?？紤]到細(xì)網(wǎng)格是網(wǎng)格尺度最小的網(wǎng)格層，相比其他層粗網(wǎng)格，將細(xì)網(wǎng)格單元合二為一或合三為一的尺度差異最小，更有助于生成分布均勻的粗網(wǎng)格。因此，當(dāng)細(xì)網(wǎng)格沿一網(wǎng)格方向的單元數(shù)無法整除2L時，本文方法僅在第1 層粗網(wǎng)格出現(xiàn)合并3 個緊鄰單元作為1 個粗網(wǎng)格的情況，如圖3 所示。

圖3 粗網(wǎng)格生成策略示意圖Fig.3 Schematic of generation of coarse grids

2.2 建立粗、細(xì)網(wǎng)格上的動態(tài)計算域

在開始控制方程的迭代求解前，還需建立動態(tài)計算域以確定第1 步參與計算的范圍。細(xì)網(wǎng)格上動態(tài)計算域的建立方法與1.4.1 節(jié)所述的擾動域推進(jìn)方法完全相同。粗網(wǎng)格的動態(tài)計算域根據(jù)細(xì)網(wǎng)格的生成。若某一細(xì)網(wǎng)格單元包含于某類動態(tài)域中，則對應(yīng)于該細(xì)網(wǎng)格單元的粗網(wǎng)格單元也被加入粗網(wǎng)格相應(yīng)的動態(tài)域中。在迭代求解中，當(dāng)細(xì)網(wǎng)格的動態(tài)計算域更新完畢，粗網(wǎng)格的動態(tài)計算域也采用同一方式更新。

2.3 增大粗、細(xì)網(wǎng)格上的對流動態(tài)域

在完成控制方程一個迭代步的求解后，需根據(jù)更新后的流場對動態(tài)計算域進(jìn)行更新。這其中，由于對流動態(tài)域的范圍直接影響求解的收斂性，因此需首先執(zhí)行增大對流動態(tài)域，將可能需要計算的單元加入對流動態(tài)域。在多重網(wǎng)格-擾動域推進(jìn)方法中，增大對流動態(tài)域的判斷算法僅在細(xì)網(wǎng)格執(zhí)行。

細(xì)網(wǎng)格上對流動態(tài)域的增大算法與1.4.2 節(jié)所述的擾動域推進(jìn)方法完全相同，對每個對流動態(tài)域邊界單元均執(zhí)行2 個子步：①判斷是否為受擾單元；②若是受擾單元，將其可能受到擾動的相鄰單元加入對流動態(tài)域。

對于子步②，擾動域推進(jìn)方法忽略了時間格式對擾動傳播速度的影響，至多允許一層緊鄰單元加入對流動態(tài)域。不過，與多重網(wǎng)格法結(jié)合后，數(shù)值擾動在一個時間步內(nèi)可通過最粗層網(wǎng)格傳向更遠(yuǎn)的地方。為避免因?qū)α鲃討B(tài)域增長緩慢而影響收斂速率，在多重網(wǎng)格-擾動域推進(jìn)方法中，細(xì)網(wǎng)格的對流動態(tài)域不再每次僅新增受擾單元的一層緊鄰單元，而是依據(jù)多重網(wǎng)格的層數(shù)調(diào)整對流動態(tài)域的新增單元層數(shù)。

考慮到粗網(wǎng)格是通過稀疏細(xì)網(wǎng)格得到的，若要求最粗一層網(wǎng)格每次僅新增受擾單元的緊鄰一層單元，則在細(xì)網(wǎng)格上對流動態(tài)域的新增單元層數(shù)La，c可估計為

式中：當(dāng)不采用多重網(wǎng)格時，L取0，則式（15）對應(yīng)于擾動域推進(jìn)方法的情況。

為提高動態(tài)計算域更新運(yùn)算的效率，式（15）并未考慮第1 層粗網(wǎng)格可能出現(xiàn)將細(xì)網(wǎng)格單元合三為一的情況。數(shù)值實(shí)驗(yàn)表明，這一簡化不會對多重網(wǎng)格-擾動域推進(jìn)方法的收斂性產(chǎn)生影響。一方面，對流動態(tài)域的增大十分迅速，迭代求解不會因計算域在一個迭代步中不足而立即發(fā)散；另一方面，初始對流動態(tài)域留有余量。初始對流動態(tài)域中包含了10 層緊鄰壁面的細(xì)網(wǎng)格單元或所有與來流條件不符的單元，這2 種情況中的對流動態(tài)域都大于會導(dǎo)致求解發(fā)散的最小求解范圍。

從加速效果的角度來看，式（15）估計的是在下一個迭代步中細(xì)網(wǎng)格上可能受到擾動的單元，其也可能會將過多單元加入對流動態(tài)域。不過，由于對流動態(tài)域縮小步驟會在下一個迭代步及時將未受擾動單元從對流動態(tài)域中移除，因此即使式（15）導(dǎo)致對流動態(tài)域偏大，也不會對多重網(wǎng)格-擾動域推進(jìn)方法的計算效率產(chǎn)生明顯影響。

2.4 調(diào)整存儲空間

擾動域推進(jìn)方法將數(shù)據(jù)分為固有信息和迭代相關(guān)信息2 類。固有信息是指在輸出時必須包含的幾何信息和流場信息，如網(wǎng)格坐標(biāo)、守恒量等。迭代相關(guān)信息是指與迭代求解方法相關(guān)的信息，如守恒量更新量、當(dāng)?shù)貢r間步長、上一迭代步的守恒量等［14］。對于固有信息，仍與傳統(tǒng)全局更新法一樣，采用靜態(tài)數(shù)組存儲。對于迭代相關(guān)信息，則采用動態(tài)數(shù)據(jù)結(jié)構(gòu)僅存儲對流動態(tài)域中單元的信息，以節(jié)省內(nèi)存需求。

多重網(wǎng)格-擾動域推進(jìn)方法在粗網(wǎng)格上也沿用這種數(shù)據(jù)存儲方式。在每個迭代步完成動態(tài)計算域的更新后，便根據(jù)對流動態(tài)域的范圍調(diào)整粗、細(xì)網(wǎng)格上迭代相關(guān)信息的存儲空間，從而同時節(jié)省粗、細(xì)網(wǎng)格上的內(nèi)存需求。

2.5 動態(tài)計算域中粗、細(xì)網(wǎng)格間的插值

各層網(wǎng)格間依靠守恒量、殘差的相互插值建立聯(lián)系。當(dāng)從第i層網(wǎng)格進(jìn)入第i+1 層時，需采用式（5）為第i+1 層單元的守恒量賦初值，采用式（7）確定第i+1 層單元的插值殘差。對于守恒量初值，無需區(qū)分所對應(yīng)第i層單元是否位于對流動態(tài)域中，仍以所有對應(yīng)單元的守恒量體積加權(quán)平均值作為第i+1 層單元的初值。對于插值殘差，未包含于第i層對流動態(tài)域中的細(xì)網(wǎng)格單元的殘差取0。當(dāng)從第i+1 層網(wǎng)格進(jìn)入第i層時，需采用式（10）計算第i+1 層單元的守恒量修正量。此外，未包含于第i+1 層對流動態(tài)域中的單元的δW2h取0。

3 算例驗(yàn)證

本文的多重網(wǎng)格-擾動域推進(jìn)方法在Win‐dows 10 系統(tǒng)下以Fortran 2003 語言實(shí)現(xiàn)。粗、細(xì)網(wǎng)格上的時間離散均采用LU-SGS 隱式格式［19］。細(xì)網(wǎng)格上的對流離散格式為熵修正的Roe 格式［20］，重構(gòu)格式為結(jié)合van Albada 限制器的MUSCL 插值［2］。細(xì)網(wǎng)格上的湍流模型選用壓縮修正的S-A 模型［21］。粗網(wǎng)格上，對流離散采用一階迎風(fēng)格式，并且不在粗網(wǎng)格上求解湍流模型方程。使用1～4 層粗網(wǎng)格：1 層粗網(wǎng)格情況即為執(zhí)行一次基礎(chǔ)循環(huán)，2～4 層粗網(wǎng)格時采用W 型循環(huán)。計算迭代的收斂閾值εc取10?7，動態(tài)計算域更新算法中參數(shù)的取值如表1 所示。多重網(wǎng)格方法在亞聲速流動中具有更優(yōu)秀的加速效果，因此選取2 個典型亞聲速問題進(jìn)行驗(yàn)證。

3.1 跨聲速無黏NACA0012 翼型

本算例考慮了以馬赫數(shù)0.85、1°迎角流經(jīng)NACA0012 翼型的跨聲速無黏流動。網(wǎng)格采用C 型拓?fù)洌膺吔绱笥?0 倍的翼型弦長以消除有限計算域的影響，網(wǎng)格量為468×120。本算例為無黏流動問題，多重網(wǎng)格-擾動域推進(jìn)方法中僅存在對流動態(tài)域。圖4 展示了3 層粗網(wǎng)格加速時多重網(wǎng)格-擾動域推進(jìn)方法中對流動態(tài)域的演化過程。圖中，ηc表示對流動態(tài)域與預(yù)設(shè)計算域的網(wǎng)格量之比；Nmax表示達(dá)到收斂態(tài)的迭代總步數(shù)。

圖4 多重網(wǎng)格-擾動域推進(jìn)方法的跨聲速無黏NACA0012 翼型求解過程（3 層粗網(wǎng)格）Fig.4 Solutions of transonic inviscid NACA0012 airfoil by DRUM with three coarse grids

由圖4 的第1 行可知，流場以來流條件初始化，因此對流動態(tài)域依據(jù)壁面邊界建立。在求解初期，擾動逐漸從壁面邊界向遠(yuǎn)場邊界傳播，對流動態(tài)域也隨之增大，如第10 步情況所示。采用3 層粗網(wǎng)格后，對流動態(tài)域每次新增8 層相鄰單元。因此，僅用10 步就將約40%的預(yù)設(shè)計算域單元納入對流動態(tài)域中。亞聲速流動中，擾動可隨聲波傳向四周；相應(yīng)地，對流動態(tài)域會隨擾動擴(kuò)展至整個預(yù)設(shè)計算域。

隨著時間推進(jìn)，遠(yuǎn)場附近單元先達(dá)到定常態(tài)，對流動態(tài)域的縮小起始于最上游遠(yuǎn)場區(qū)域，如17%Nmax的情況所示。隨后，對流動態(tài)域逐漸向壁面縮小，如38%Nmax的情況所示。激波附近存在的流動特性大梯度會導(dǎo)致求解更難以達(dá)到定常態(tài)。因此，在73%Nmax～96%Nmax的情況中，對流動態(tài)域的縮小被翼型上、下兩道激波所阻擋。在計算結(jié)束時，即100%Nmax的情況，對流動態(tài)域消失，翼型上、下表面的兩道激波均被清晰捕捉。

圖4 的下半子圖對比了不同網(wǎng)格層上對流動態(tài)域的區(qū)別。由于粗網(wǎng)格上的動態(tài)域是依據(jù)細(xì)網(wǎng)格上的生成，因此粗網(wǎng)格上對流動態(tài)域的形狀與細(xì)網(wǎng)格上的保持一致。但隨著網(wǎng)格粗化，對流動態(tài)域所包含的范圍會有所擴(kuò)大。這是由于，即使某一粗網(wǎng)格單元所對應(yīng)的多個細(xì)網(wǎng)格單元中，僅有一個包含在細(xì)網(wǎng)格上的動態(tài)域中，該粗網(wǎng)格單元仍會被納入該層粗網(wǎng)格的動態(tài)域。

本算例中，翼型的上、下表面均會產(chǎn)生激波。從圖4 第1 行中100%Nmax情況的局部放大圖可以看出，翼型上、下表面的兩道激波均能被清晰捕捉；受1°迎角的影響，上表面激波更靠近后緣，強(qiáng)度也較高。圖5 將采用3 層粗網(wǎng)格的多重網(wǎng)格-擾動域推進(jìn)方法所得壁面壓力分布與傳統(tǒng)多重網(wǎng)格方法的、文獻(xiàn)［22-23］的結(jié)果進(jìn)行了對比。對比表明，本文方法與傳統(tǒng)多重網(wǎng)格方法的結(jié)果完全重合，本文方法提供的數(shù)值結(jié)果與文獻(xiàn)［22-23］的數(shù)值、實(shí)驗(yàn)結(jié)果也能很好吻合。

圖5 NACA0012 翼型壓力系數(shù)分布對比Fig.5 Comparison of pressure coefficient distributions on NACA0012 airfoil

表2 將本文采用3 層粗網(wǎng)格所得翼型表面氣動力系數(shù)與文獻(xiàn)［22］的結(jié)果進(jìn)行了對比。其中，CL、CD分別表示升力系數(shù)和阻力系數(shù)；ΔCL/CL、ΔCD/CD分別為本文所得升力系數(shù)和阻力系數(shù)與文獻(xiàn)［22］結(jié)果的相對偏差。由對比可知，本文數(shù)值結(jié)果與文獻(xiàn)［22］的氣動力偏差均在1.92%以內(nèi)。因此，本文所使用的傳統(tǒng)多重網(wǎng)格方法和多重網(wǎng)格-擾動域推進(jìn)方法能夠提供正確的無黏流動數(shù)值解。

表2 NACA0012 翼型的氣動力對比Table 2 Comparison of forces on NACA0012 airfoil

圖6 對比了不同粗網(wǎng)格層數(shù)下對流動態(tài)域的網(wǎng)格量變化。由圖可知，所有采用多重網(wǎng)格情況的ηc增長率均會大于未采用多重網(wǎng)格的情況，且粗網(wǎng)格層數(shù)越多，ηc增長率越快。這一現(xiàn)象有2 方面原因：其一，多重網(wǎng)格-擾動域推進(jìn)方法要求對流動態(tài)域的新增單元層數(shù)與粗網(wǎng)格層數(shù)成正比；其二，粗、細(xì)網(wǎng)格間的循環(huán)迭代能加速擾動的傳播，也會促使對流動態(tài)域增長更快。

圖6 NACA0012 翼型問題對流動態(tài)域網(wǎng)格量變化曲線Fig.6 Evolution of number of cells in advective dy‐namic computational domain for NACA0012 airfoil

圖7（a）和圖7（b）分別展示了本算例中傳統(tǒng)全局更新方法和擾動域推進(jìn)方法采用0～4 層粗網(wǎng)格求解時的收斂曲線。由圖7（a）中不同粗網(wǎng)格層數(shù)的收斂曲線對比可知，多重網(wǎng)格方法有助于消除求解后期收斂曲線的振蕩，從而大幅降低總迭代步數(shù)。觀察圖7（b）中的收斂曲線可知，由于動態(tài)計算域的縮小，擾動域推進(jìn)方法在求解后期也可起到與多重網(wǎng)格相同的效果，因此也同樣有大幅降低總迭代步數(shù)的效果。

圖7 NACA0012 翼型問題收斂曲線Fig.7 Convergence histories of NACA0012 airfoil

表3 對比了本算例中采用不同層粗網(wǎng)格擾動域推進(jìn)方法的加速效果與求解精度。其中，擾動域推進(jìn)時間節(jié)省是指相同粗網(wǎng)格層數(shù)下多重網(wǎng)格-擾動域推進(jìn)方法相對于傳統(tǒng)全局更新方法的計算總耗時節(jié)省量；求解精度是指相同粗網(wǎng)格層數(shù)下多重網(wǎng)格-擾動域推進(jìn)方法相對于傳統(tǒng)全局更新方法的升力、阻力系數(shù)最大偏差。

表3 NACA0012 翼型的擾動域推進(jìn)效果Table 3 Acceleration effects stemming of DRUM for NACA0012 airfoil

得益于同時降低細(xì)網(wǎng)格與粗網(wǎng)格上的計算量并減少達(dá)到收斂所需的迭代步數(shù)，多重網(wǎng)格-擾動域推進(jìn)方法在傳統(tǒng)多重網(wǎng)格方法加速的基礎(chǔ)上還可再提升25.0%～39.9%的計算效率。對比不同粗網(wǎng)格層數(shù)的情況可知，隨著粗網(wǎng)格層數(shù)的增加，多重網(wǎng)格-擾動域推進(jìn)方法的時間節(jié)省量會逐漸降低，但求解精度也逐漸提升。這一現(xiàn)象的主要原因在于，多重網(wǎng)格和擾動域推進(jìn)2 種方法都有在求解后期消除殘差振蕩的效果，如圖7所示；而當(dāng)粗網(wǎng)格層數(shù)>1 時，傳統(tǒng)多重網(wǎng)格方法已可完全消除殘差振蕩現(xiàn)象；此時，擾動域推進(jìn)方法的迭代步數(shù)節(jié)省量顯著減小，其加速效果主要來源于通過對流動態(tài)域的縮小同時降低細(xì)網(wǎng)格與粗網(wǎng)格上的計算量。因此，擾動域推進(jìn)時間節(jié)省量隨著粗網(wǎng)格層數(shù)的增加而降低。

此外，表3 的效率與精度結(jié)果對應(yīng)于動態(tài)計算域更新中閾值恒定的情況。多重網(wǎng)格-擾動域推進(jìn)方法的計算效率還可通過縮小對流刪除閾值εd、增大上游容差角θd等閾值來提升。

表4 給出了多重網(wǎng)格-擾動域推進(jìn)方法相比于未使用多重網(wǎng)格的傳統(tǒng)全局更新方法的時間節(jié)省量。本算例中，采用3 層粗網(wǎng)格的情況獲得了最優(yōu)的加速效果，可節(jié)省90.8%的計算時間。此外，本算例網(wǎng)格在I方向的單元數(shù)為468，無法整除23。因此，采用3 層粗網(wǎng)格情況獲得最優(yōu)加速效果也證明了本文所提出的新型粗網(wǎng)格生成策略切實(shí)可行，可有效提高多重網(wǎng)格-擾動域推進(jìn)方法的適用性。

表4 NACA0012 翼型的總加速效果Table 4 Total acceleration effects of DRUM for NACA0012 airfoil

綜上所述，本算例驗(yàn)證了多重網(wǎng)格-擾動域推進(jìn)方法對亞聲速無黏流動模擬的適用性。得益于同時降低細(xì)、粗網(wǎng)格的計算量并減少總迭代步數(shù)，多重網(wǎng)格-擾動域推進(jìn)方法相比于傳統(tǒng)多重網(wǎng)格方法可展現(xiàn)出更為顯著的加速效果。

3.2 跨聲速湍流RAE2822 翼型

本算例考慮了跨聲速RAE2822 翼型風(fēng)洞試驗(yàn)［24］的第6 種狀態(tài)，即來流馬赫數(shù)為0.729，雷諾數(shù)為6.5×106，迎角為2.31°。使用NASA 網(wǎng)站所發(fā)布的網(wǎng)格進(jìn)行數(shù)值模擬，網(wǎng)格量為368×64。本算例為湍流流動，多重網(wǎng)格-擾動域推進(jìn)方法中包含對流、黏性2 類動態(tài)計算域。圖8 展示了采用3 層粗網(wǎng)格時求解跨聲速湍流RAE2822 翼型中動態(tài)域的演化過程；其中，ηv表示黏性動態(tài)域與預(yù)設(shè)計算域的網(wǎng)格量之比。

圖8 多重網(wǎng)格-擾動域推進(jìn)方法的跨聲速湍流RAE2822 翼型求解過程（3 層粗網(wǎng)格）Fig.8 Solutions of transonic turbulent RAE2822 airfoil by DRUM with three coarse grids

由圖8 的第1 行可知，初始流場設(shè)定為具有來流狀態(tài)的均勻流場，因此對流、黏性2 類動態(tài)域均依據(jù)壁面邊界建立。隨著擾動向周圍流動傳播，包圍著翼型的這2 類動態(tài)計算域也逐漸擴(kuò)大。采用多重網(wǎng)格方法后，對流動態(tài)域隨壁面擾動迅速向遠(yuǎn)處擴(kuò)張。在第5 步時，對流動態(tài)域已經(jīng)包含54.5%的預(yù)設(shè)計算域單元。由4%Nmax的情況可知，對流動態(tài)域會增長到整個預(yù)設(shè)計算域，而最大黏性動態(tài)域只停留在近壁區(qū)。由于流場解的變化，當(dāng)對流動態(tài)域保持在最大范圍時，黏性動態(tài)域會先開始緩慢縮小。當(dāng)流場求解開始收斂時，對流動態(tài)域會從最上游遠(yuǎn)場區(qū)域開始逐漸向壁面收縮，如33%Nmax的情況所示，并帶動黏性動態(tài)域一起縮小，如57%Nmax的情況所示。求解結(jié)束時，2 類動態(tài)計算域消失，流場達(dá)到定常態(tài)，上壁面的激波被很好地捕捉，如100%Nmax的情況所示。

圖8 的下半子圖對比了不同網(wǎng)格層上對流、粘性動態(tài)域的區(qū)別。由于網(wǎng)格的差異，2 類動態(tài)域的范圍都會隨著網(wǎng)格的粗化而有所變大。由57%Nmax的情況可知，由于上翼面激波附近的流動變化劇烈，求解較難收斂，因此，上翼面動態(tài)域的縮小被激波所阻擋，而下翼面的動態(tài)域縮小明顯快于上翼面。

表5 和圖9 將本文采用3 層粗網(wǎng)格的數(shù)值結(jié)果與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了對比，表5 中CN表示法向力系數(shù)，ΔCN/CN表示法向力系數(shù)偏差。本文所計算的傳統(tǒng)多重網(wǎng)格方法和多重網(wǎng)格-擾動域推進(jìn)方法所得的氣動力均能與實(shí)驗(yàn)結(jié)果很好地吻合，氣動力偏差低于1.6%。圖9 中，傳統(tǒng)多重網(wǎng)格方法和多重網(wǎng)格-擾動域推進(jìn)方法預(yù)測的翼型表面壓力分布完全重合，兩者與實(shí)驗(yàn)數(shù)據(jù)也吻合良好。因此，本文所使用的傳統(tǒng)多重網(wǎng)格方法和多重網(wǎng)格-擾動域推進(jìn)方法能夠提供正確的湍流流動數(shù)值解。

表5 RAE2822 翼型的氣動力對比Table 5 Comparison of forces on RAE2822 airfoil

圖9 RAE2822 翼型壓力系數(shù)分布對比Fig.9 Comparison of pressure coefficient distributions on RAE2822 airfoil

圖10 對比了RAE2822 翼型算例中多重網(wǎng)格對對流、黏性2 類動態(tài)域網(wǎng)格量變化的影響。多重網(wǎng)格方法影響對流動態(tài)域的規(guī)律與無黏情況一致。對于黏性動態(tài)域，多重網(wǎng)格-擾動域推進(jìn)方法中黏性動態(tài)域仍保持每次只新增緊鄰的一層單元，因此粗網(wǎng)格層數(shù)的不同并未對黏性動態(tài)域的增大速率產(chǎn)生影響。不過，粗網(wǎng)格層數(shù)為2 層、3 層情況的最大黏性動態(tài)域會明顯小于未使用多重網(wǎng)格的情況。這一現(xiàn)象主要源自黏性動態(tài)域的范圍與流動參數(shù)密切相關(guān)，而多重網(wǎng)格方法可加速流場的收斂；即在相同迭代步時，采用粗網(wǎng)格情況對應(yīng)的流場解更接近定常解。因此，采用多重網(wǎng)格與否會對最大黏性動態(tài)域的范圍產(chǎn)生影響。

圖10 RAE2822 翼型動態(tài)域網(wǎng)格量變化曲線Fig.10 Evolution of number of cells in dynamic compu‐tational domains for RAE2822 airfoil

圖11 展示了RAE2822 算例中考慮的所有情況的收斂曲線。所有情況均達(dá)到了收斂條件。對于湍流問題，多重網(wǎng)格-擾動域推進(jìn)方法也同樣具有顯著的加速效果。加速效果源自3 方面：其一，對流動態(tài)域的自適應(yīng)增大、縮小可同時降低細(xì)、粗網(wǎng)格的計算量；其二，黏性動態(tài)域可有效甄別黏性效應(yīng)主導(dǎo)區(qū)域，非黏性效應(yīng)主導(dǎo)區(qū)域未計算控制方程黏性項(xiàng)；其三，達(dá)到相同收斂條件所需的總迭代步數(shù)降低。

圖11 RAE2822 翼型收斂曲線Fig.11 Convergence histories of RAE2822 airfoil

表6 展示了本算例中采用不同層粗網(wǎng)格擾動域推進(jìn)方法的加速效果與求解精度。多重網(wǎng)格-擾動域推進(jìn)方法在傳統(tǒng)多重網(wǎng)格方法加速的基礎(chǔ)上還可再提升53.9%～57.9%的計算效率。與無黏NACA0012 情況相同，多重網(wǎng)格-擾動域推進(jìn)方法的求解精度會隨粗網(wǎng)格層數(shù)的增加而提升。與無粘NACA0012 情況不同的是，本算例中采用3 層粗網(wǎng)格情況所獲得的擾動域推進(jìn)時間節(jié)省量大于采用2 層粗網(wǎng)格或未采用多重網(wǎng)格的情況。由圖11 可知，這主要是由于采用3 層粗網(wǎng)格的傳統(tǒng)全局更新方法的加速效果不佳而導(dǎo)致的。采用3 層粗網(wǎng)格時，傳統(tǒng)全局更新方法的收斂曲線出現(xiàn)了明顯的振蕩，導(dǎo)致其相比于采用2 層粗網(wǎng)格的情況并未提升效率。

表6 RAE2822 翼型的擾動域推進(jìn)時間效果Table 6 Acceleration effects stemming of DRUM for RAE2822 airfoil

表7 給出了本算例中多重網(wǎng)格-擾動域推進(jìn)方法相比于未使用多重網(wǎng)格的傳統(tǒng)全局更新方法的時間節(jié)省量。對比可知，采用3 層粗網(wǎng)格的多重網(wǎng)格-擾動域推進(jìn)方法獲得了最優(yōu)的加速效果，可節(jié)省90.4%的計算時間。

表7 RAE2822 翼型的總加速效果Table 7 Total acceleration effects of DRUM for RAE2822 airfoil

綜上所述，本算例驗(yàn)證了多重網(wǎng)格-擾動域推進(jìn)方法對亞聲速湍流流動模擬的適用性。得益于同時降低細(xì)、粗網(wǎng)格的計算量、無黏-黏性流動自適應(yīng)分區(qū)計算和減少總迭代步數(shù)，多重網(wǎng)格-擾動域推進(jìn)方法可在傳統(tǒng)多重網(wǎng)格方法的加速效果上，進(jìn)一步顯著提升計算效率。

4 結(jié)論

在作者前期提出的擾動域推進(jìn)方法計算框架中，融入多重網(wǎng)格加速技術(shù)，提出了多重網(wǎng)格-擾動域推進(jìn)方法。著重建立了粗網(wǎng)格上動態(tài)計算域的演化邏輯，設(shè)計了對結(jié)構(gòu)網(wǎng)格更為普適的新型粗網(wǎng)格生成策略，并對典型算例進(jìn)行了驗(yàn)證演示，得到以下結(jié)論：

1）多重網(wǎng)格-擾動域推進(jìn)方法通過采用隨數(shù)值擾動傳播而增大、隨求解收斂而縮小的動態(tài)計算域，對于可壓縮流動，尤其是亞聲速流動，具有顯著的加速效果。典型算例驗(yàn)證結(jié)果表明，多重網(wǎng)格-擾動域推進(jìn)方法可在傳統(tǒng)多重網(wǎng)格方法加速的基礎(chǔ)上再將計算效率提升57.9%。

2）多重網(wǎng)格-擾動域推進(jìn)方法的加速原因包括3 方面：其一，對流動態(tài)域的自適應(yīng)增大、縮小可同時降低細(xì)、粗網(wǎng)格的計算量；其二，黏性動態(tài)域可有效甄別黏性效應(yīng)主導(dǎo)區(qū)域，非黏性效應(yīng)主導(dǎo)區(qū)域無需計算流動控制方程的黏性項(xiàng)；其三，達(dá)到相同收斂條件所需的總迭代步數(shù)降低。

3）提出的新型粗網(wǎng)格生成策略，可消除結(jié)構(gòu)網(wǎng)格生成粗網(wǎng)格對單元數(shù)目的限制，提高了多重網(wǎng)格-擾動域推進(jìn)方法的適用性。