胡姝瑤，蔣崇文，李椿萱

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

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

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

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

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

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

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

1.1 控制方程

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

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

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

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

1.2 定?？蓧嚎s流動(dòng)的有限體積方法

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

式中：殘差R可表示為

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

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

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

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

多重網(wǎng)格方法需在一層細(xì)網(wǎng)格和若干層粗網(wǎng)格上，按指定規(guī)律循環(huán)求解式（3）。對(duì)于結(jié)構(gòu)網(wǎng)格，常用的粗網(wǎng)格生成方法是將細(xì)網(wǎng)格沿網(wǎng)格方向稀疏一倍［2］。常用的多重網(wǎng)格循環(huán)有V 型、W 型等，其均由基礎(chǔ)循環(huán)組成［18］。基礎(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)格的個(gè)數(shù)；下標(biāo)h、2h分別表示細(xì)網(wǎng)格和其稀疏一倍所得粗網(wǎng)格。粗網(wǎng)格的體積|Ω|2h可表示為

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

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

粗網(wǎng)格所采用的時(shí)間推進(jìn)格式與細(xì)網(wǎng)格保持一致。粗網(wǎng)格上的迭代步數(shù)可任意給定，文獻(xiàn)［13］指出當(dāng)以LU-SGS 格式推進(jìn)求解時(shí)在粗網(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)格式求解對(duì)流通量時(shí)，文獻(xiàn)［18］推薦粗網(wǎng)格通過(guò)式（11）獲得新解后，先迭代一步再將其修正量插值給更密網(wǎng)格以獲得最佳的多重網(wǎng)格性能。

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

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

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

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

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

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

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

1.4.1 建立動(dòng)態(tài)計(jì)算域

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

1.4.2 增大對(duì)流動(dòng)態(tài)域

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

1.4.3 縮小對(duì)流動(dòng)態(tài)域

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

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

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

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

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

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

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

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

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

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

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

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

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

對(duì)于一個(gè)結(jié)構(gòu)網(wǎng)格，文獻(xiàn)［2］采用沿一網(wǎng)格方向?qū)⒁粚?duì)緊鄰單元合二為一的方式生成粗網(wǎng)格，即對(duì)于一個(gè)單元數(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；若不滿足，則會(huì)出現(xiàn)單個(gè)無(wú)法合并單元，影響多重網(wǎng)格方法的使用和效果。

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

式中：mod（?）代表取模。第i層網(wǎng)格的單元I（i）對(duì)應(yīng)于第i+1 層網(wǎng)格單元的標(biāo)號(hà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ù)無(wú)法整除2L時(shí)，本文方法僅在第1 層粗網(wǎng)格出現(xiàn)合并3 個(gè)緊鄰單元作為1 個(gè)粗網(wǎng)格的情況，如圖3 所示。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 算例驗(yàn)證

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

3.1 跨聲速無(wú)黏NACA0012 翼型

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

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

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

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

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

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

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

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

表2 NACA0012 翼型的氣動(dòng)力對(duì)比Table 2 Comparison of forces on NACA0012 airfoil

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

圖6 NACA0012 翼型問(wèn)題對(duì)流動(dòng)態(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)全局更新方法和擾動(dòng)域推進(jìn)方法采用0～4 層粗網(wǎng)格求解時(shí)的收斂曲線。由圖7（a）中不同粗網(wǎng)格層數(shù)的收斂曲線對(duì)比可知，多重網(wǎng)格方法有助于消除求解后期收斂曲線的振蕩，從而大幅降低總迭代步數(shù)。觀察圖7（b）中的收斂曲線可知，由于動(dòng)態(tài)計(jì)算域的縮小，擾動(dòng)域推進(jìn)方法在求解后期也可起到與多重網(wǎng)格相同的效果，因此也同樣有大幅降低總迭代步數(shù)的效果。

圖7 NACA0012 翼型問(wèn)題收斂曲線Fig.7 Convergence histories of NACA0012 airfoil

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

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

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

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

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

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

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

3.2 跨聲速湍流RAE2822 翼型

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

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

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

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

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

表5 RAE2822 翼型的氣動(dòng)力對(duì)比Table 5 Comparison of forces on RAE2822 airfoil

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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