基于降階模型和梯度優(yōu)化的流場(chǎng)加速收斂方法

曹文博，劉溢浪，張偉偉

西北工業(yè)大學(xué) 航空學(xué)院，西安 710072

計(jì)算流體力學(xué)（Computational Fluid Dynamics，CFD）在近40 年的時(shí)間內(nèi)發(fā)展迅速，成為流體力學(xué)理論研究及工程氣動(dòng)設(shè)計(jì)和分析的重要手段。如今，科學(xué)研究以及工程問(wèn)題的日益復(fù)雜化和精細(xì)化給CFD 算法的效率提出了更高要求。因此，發(fā)展高效的加速收斂方法仍是計(jì)算流體力學(xué)中至關(guān)重要的研究方向。

迄今為止，研究者已提出了許多經(jīng)典的加速收斂算法，例如局部時(shí)間步長(zhǎng)、隱式殘差光順[1］、多重網(wǎng)格[2］、焓阻尼方法[3］、低馬赫數(shù)預(yù)處理[4］等。最小化迭代求解器收斂殘差外插方法的一般過(guò)程是根據(jù)迭代求解過(guò)程中解快照外插得到殘差更小的解，而后基于該預(yù)測(cè)值繼續(xù)迭代求解生成新的解快照，以此循環(huán)。其優(yōu)勢(shì)在于應(yīng)用簡(jiǎn)單，不必修改求解器算法，并能和經(jīng)典加速收斂算法結(jié)合使用。

早期的外插方法幾乎都是基于線性迭代假設(shè)xj+1=Axj+b（其中A、b 未知），通過(guò)迭代過(guò)程中的解快照{(diào) xn}近似該系統(tǒng)解s：s=As+b。這些方法主要分為2 大類：①epsilon 算法；②多項(xiàng)式方法。向量epsilon 算法[5］（Vector Epsilon Algorithm，VEA）和拓?fù)鋏psilon 算法[6］（Topological Epsilon Algorithm，TEA）是最早的外插方法，其借助求和方法和一些復(fù)雜的反演公式將緩慢收斂的序列轉(zhuǎn)化為快速收斂的序列，但這類方法的復(fù)雜性使得其應(yīng)用較困難。多項(xiàng)式是一種更簡(jiǎn)單的外插方法，最小多項(xiàng)式外插[7］（Minimal Polynomial Extrapolation，MPE）和降階外插[8-9］（Reduced Rank Extrapolation，RRE）是其中主要的2 種方法。這些方法的比較可見Smith 等[10］的綜 述。Sidi[11］闡 述 了MPE、RRE 方 法 與Arnoldi方法[12］和GMRES 方法[13］的聯(lián)系。這些方法之后被發(fā)展并應(yīng)用于Euler 方程和其他迭代求解器的加速求解[11，14］。Djeddi 等[15］假設(shè)迭代線性收斂，將方程殘值表示為解的線性函數(shù)R(U )≈AU-b，并結(jié)合降階模型近似該系統(tǒng)的解，實(shí)現(xiàn)了雷諾平均Navier-Stokes 方程（Reynlods-Averaged Navier-Stokes，RANS）的加速收斂。

隨 著 動(dòng) 態(tài) 模 態(tài) 分 解[16-17］（Dynamic Mode Decom-position，DMD）成為流場(chǎng)分析的熱點(diǎn)，基于DMD 外插方法也引起了研究者的興趣。Liu等[18］使用DMD 預(yù)測(cè)雙時(shí)間步格式中內(nèi)迭代的初始條件從而加速RANS 方程求解。Andersson[19］使用DMD 預(yù)測(cè)了線性迭代假設(shè)下的收斂解。另外，基于DMD 的模態(tài)多重網(wǎng)格（Mode Multigrid，MMG）[20-21］也被驗(yàn)證是一種有效的加速收斂方法，這種方法認(rèn)為DMD 一階模態(tài)對(duì)應(yīng)于穩(wěn)態(tài)流動(dòng)，使用DMD 過(guò)濾振蕩模態(tài)以加速收斂。

降階模型（Reduced-Order Model，ROM）外插是另外一種以較小計(jì)算成本獲得合理外插結(jié)果的方法?；赑OD 的降階模型算法是近30 年發(fā)展起來(lái)的，它將高維控制方程轉(zhuǎn)化為給定POD模態(tài)下的一組低階方程，從而極大地減少了計(jì)算時(shí)間，該方法已廣泛應(yīng)用于各種領(lǐng)域[22-26］。然而，由于該方法需要用高保真的方法（通常是髙耗時(shí)的CFD 計(jì)算）預(yù)先獲取流場(chǎng)快照，早期這類方法只應(yīng)用于有較多變化參數(shù)且需要多次重復(fù)求解的CFD 問(wèn)題。后來(lái)，基于POD 的降階模型被發(fā)展并應(yīng)用于全隱式求解器中內(nèi)迭代過(guò)程的初值預(yù)測(cè)，從而加速了非定常問(wèn)題求解[27-29］。這些方法中，流場(chǎng)快照不是預(yù)先給定的，而是在非定常時(shí)間推進(jìn)求解過(guò)程中收集并不斷被更新，避免了已有工作[22-26］中預(yù)先計(jì)算流場(chǎng)快照帶來(lái)的昂貴代價(jià)，但這些加速方法均局限于隱式求解器。之后，Rapun 等[30］提出了當(dāng)?shù)豍OD 方法用于加速拋物型問(wèn)題的非定常求解過(guò)程，該方法使用數(shù)值格式和降階模型交替求解方程，并通過(guò)對(duì)降階模型近似誤差的先驗(yàn)估計(jì)確定降階模型的求解時(shí)長(zhǎng)，保證了魯棒性和有效性，并被應(yīng)用于各類穩(wěn)態(tài)[31］和瞬態(tài)[32-34］流動(dòng)的加速求解。

基于POD 降階模型盡管已經(jīng)被應(yīng)用于一些方程的加速求解，但CFD 算法的復(fù)雜性，特別是基于非結(jié)構(gòu)網(wǎng)格有限體積法的復(fù)雜性，且CFD 求解器構(gòu)造對(duì)應(yīng)的嵌入式降階模型較為困難，因此，該類方法并未被應(yīng)用于復(fù)雜問(wèn)題的RANS 方程加速求解。本文對(duì)已有方法進(jìn)行了多種改進(jìn)，使得降階模型中只包含關(guān)于流場(chǎng)快照的簡(jiǎn)單代數(shù)運(yùn)算，形成了一個(gè)加速RANS 方程穩(wěn)態(tài)求解過(guò)程的通用框架。具體的改進(jìn)主要包括：①流場(chǎng)快照中不僅包含網(wǎng)格格心處的流場(chǎng)值，還包含有限體積法離散格式中需要計(jì)算的面心流場(chǎng)值、面心左右兩側(cè)的流場(chǎng)值、面心流場(chǎng)梯度，避免了對(duì)POD 模態(tài)的空間離散；②不將控制方程轉(zhuǎn)化為常微分方程組而后時(shí)間推進(jìn)求解，取而代之的是通過(guò)梯度優(yōu)化尋找POD 子空間內(nèi)殘差最小解，這種方法具有更強(qiáng)的靈活性和魯棒性；③只計(jì)算少量網(wǎng)格殘值來(lái)優(yōu)化模態(tài)系數(shù)從而大大減少降階模型優(yōu)化求解的計(jì)算成本。

1 數(shù)值方法

1.1 控制方程

考慮了基于SA 湍流模型的二維雷諾平均Navier-Stokes方程。RANS 方程組的積分形式為

式 中：V=[ux，uy]T為 流 體 運(yùn) 動(dòng) 速 度；Ω 為 控 制體，其邊界為?Ω；n 為邊界的單位外法G(Q)線方向；Q 為守恒變量；F(Q)為無(wú)黏通量；G(Q)為黏性通量；ρ、e0、P、T 分別為流體的密度、單位體積的總能、壓強(qiáng)和溫度；nx、ny分別為邊界外法線方向n 在x、y 方 向 的 分 量；σij為 流 體 間 的 應(yīng) 力，i 為作用面的方向，j 為作用方向；qj為單位質(zhì)量的體積熱在3 個(gè)方向上的分布；γ 為氣體的比熱比，對(duì)于理想氣體γ=1.4；μ 為分子黏性系數(shù)，可利用Sutherland 公式得到，即

μt為湍流黏性系數(shù)，通過(guò)求解湍流模型或亞格子模型方程得到；傳熱項(xiàng)中Pr 為層流Prandtl 數(shù)，Prt為湍流Prandtl 數(shù)，文中分別取值為0.87 和0.9。

對(duì)于理想氣體，有

由于SA 方程的復(fù)雜性，本文不對(duì)SA 方程做降階處理，即求解降階模型時(shí)渦黏固定不變。

1.2 數(shù)值格式

有限體積法中，控制方程通常被離散為數(shù)值形式：

式中：m 為網(wǎng)格編號(hào)；Vm為網(wǎng)格m 的體積；Qm為網(wǎng)格m 格心處的守恒變量；F(m)為網(wǎng)格m 擁有的面的集合；p 為面的編號(hào)；Ump、?Ump、ULmp、URmp分別為面心處的流場(chǎng)值、流場(chǎng)梯度、面心左側(cè)的流場(chǎng)值、面心右側(cè)的流場(chǎng)值；ΔSmp為面積分別為面心無(wú)黏通量和黏性通量，可采用不同的通量格式求解。

1.3 基于POD 的降階模型

降階模型用于將控制方程轉(zhuǎn)化為給定POD模態(tài)下的低階常微分方程組，從而極大地減少了計(jì)算時(shí)間。首先收集CFD 計(jì)算得到的流場(chǎng)快照X=[U1，U2，…，UN]，對(duì)其進(jìn)行奇異值分解得到POD 模態(tài)（通常對(duì)奇異值較小的模態(tài)截?cái)啵?/p>

在本文中，將構(gòu)建N-S 方程的降階模型，實(shí)現(xiàn)加速收斂目的，降階模型的數(shù)值格式必須與CFD 一致，因此將基于有限體積法的數(shù)值格式式（8）投 影 到 POD 模 式 上 。 注 意 到Ump、?Ump、ULmp、URmp均為不同網(wǎng)格格心值Um的線性疊加，且POD 快照均為流場(chǎng)快照的線性疊加，若

式 中：Xm、Xmp、X?mp、XLmp、XRmp和 Φm、Φmp、Φ?mp、ΦLmp、ΦRmp分 別 為Um、Ump、?Ump、ULmp、URmp的 流 場(chǎng)快照和POD 模態(tài)。式（16）表明，僅需對(duì)格心流場(chǎng)值Um做POD 分析，其余變量的POD 模態(tài)均可以通過(guò)式（16）得到，其中T 通過(guò)式（17）計(jì)算：

通過(guò)式（16）和式（17）計(jì)算得到POD 模態(tài)后，則可得到Galerkin 展開式為

1.4 梯度優(yōu)化方法

梯度下降是一種用來(lái)最小化目標(biāo)函數(shù)R(ξ)的方法，根據(jù)目標(biāo)函數(shù)對(duì)參數(shù)ξ 的負(fù)梯度-?R(ξ)以給定學(xué)習(xí)率ε 更新參數(shù)的值。

文中，降階模型部分在Python 中實(shí)現(xiàn)，使用自動(dòng)微分計(jì)算梯度?R(ξ)，使用AdamW 算法進(jìn)行梯度優(yōu)化。特別地，最后一個(gè)流場(chǎng)快照對(duì)應(yīng)的模態(tài)系數(shù)顯然較逼近最優(yōu)解，因此可將其作為優(yōu)化參數(shù)的初值ξ0以加速優(yōu)化過(guò)程。

需要注意的是，盡管相比于CFD，降階模型采用式（18）計(jì)算Ump、?Ump、ULmp、URmp減少了計(jì)算量，但由于降階模型是使用梯度優(yōu)化求解，降階模型求解的計(jì)算量仍會(huì)較大，特別是當(dāng)網(wǎng)格量較大時(shí)。本文在求解式（19）時(shí)，只計(jì)算少量網(wǎng)格上的殘差以優(yōu)化參數(shù)ξ，這大大減少了降階模型求解的計(jì)算量，使得流場(chǎng)求解的整個(gè)過(guò)程中，降階模型的計(jì)算量不到CFD 求解的1/10，這是本文方法能夠有效的一個(gè)關(guān)鍵點(diǎn)。

加速收斂方法包含的步驟：①基于流場(chǎng)初值U0m進(jìn)行短時(shí)間內(nèi)的CFD 迭代，保存求解過(guò)程中Um、Ump、?Ump、ULmp、URmp的快照；②從快照中提取能量最高的POD 模態(tài)，并使用Galerkin 投影將數(shù)值格式投影到POD 模態(tài)上得到降階模型式（19）③使用梯度優(yōu)化最小化降階模型平均殘差R(ξ)，從而得到殘差更低的流場(chǎng)，更新U0m；④重復(fù)步驟①～步驟③，直至收斂。

本方法中主要超參數(shù)有：計(jì)算降階模型平均殘差R(ξ)時(shí)使用的網(wǎng)格數(shù)M 和網(wǎng)格的采樣方式、流場(chǎng)快照數(shù)目N、快照保存間隔K。

2 典型算例驗(yàn)證

本節(jié)將采用多個(gè)典型流場(chǎng)算例驗(yàn)證所提加速收斂方法在復(fù)雜流動(dòng)求解中的有效性以及參數(shù)敏感性。分別為NACA0012 翼型繞流的無(wú)黏算例和湍流算例、S809 翼型繞流湍流算例以及ONERA M6 機(jī)翼的亞聲速湍流算例。算例均采用二階有限體積法，無(wú)黏通量采用Roe 格式計(jì)算，黏性通量采用中心格式計(jì)算，時(shí)間推進(jìn)方法為隱式Gauss-Seidel 迭代。

2.1 NACA0012 翼型繞流無(wú)黏算例

NACA0012 翼型無(wú)黏算例中，來(lái)流馬赫數(shù)Ma=0.3，攻角α=0°，CFL=5，網(wǎng)格總數(shù)16 155，近壁面網(wǎng)格如圖1 所示。每隔200 步進(jìn)行一次降階模型求解，快照數(shù)目N=40，保存間隔K=5。分別隨機(jī)選取200、500、1 000、2 000 個(gè)網(wǎng)格點(diǎn)計(jì)算降階模型的平均殘差，以測(cè)試網(wǎng)格點(diǎn)數(shù)對(duì)加速方法的影響。計(jì)算得到的流場(chǎng)殘值收斂歷程如圖2 所示，不同方法所用的時(shí)間統(tǒng)計(jì)如表1所示。

圖1 NACA0012 翼型無(wú)黏網(wǎng)格Fig.1 Inviscid grid of NACA0012 airfoil

圖2 不同網(wǎng)格點(diǎn)數(shù)的加速方法收斂歷程Fig.2 Residual history with different numbers of grids

表1 不同網(wǎng)格點(diǎn)數(shù)的加速方法計(jì)算時(shí)間Table 1 Calculation time with different numbers of grids

由圖2 和表1 可知，所提的加速收斂方法對(duì)于無(wú)黏流動(dòng)能達(dá)到約1.5 倍的加速收斂效果，從收斂歷程來(lái)看，該方法對(duì)于網(wǎng)格數(shù)的選取不敏感，但由于降階模型中計(jì)算的網(wǎng)格數(shù)越多，降階模型計(jì)算時(shí)間越長(zhǎng)，因此網(wǎng)格數(shù)M不宜太大。圖3 給出了加速方法和初始CFD 方法的壓力系數(shù)Cp分布對(duì)比，可看出2 種方法計(jì)算得到的翼型表面壓力分布一致。事實(shí)上，因?yàn)榻惦A模型中的快照來(lái)源于CFD 計(jì)算并且被不斷更新，所以加速收斂方法的精度仍然取決于CFD 計(jì)算過(guò)程，不會(huì)損失任何精度。

圖3 翼型表面壓力分布對(duì)比Fig.3 Comparison of surface pressure coefficient

2.2 NACA0012 翼型繞流湍流算例

本節(jié)計(jì)算NACA0012 翼型的湍流算例用于測(cè)試加速收斂方法在復(fù)雜工程湍流問(wèn)題中的有效性。來(lái)流馬赫數(shù)Ma=0.3，攻角α=8°，雷諾數(shù)Re=3×10-6，CFL=2，網(wǎng)格總數(shù)為76 356，近壁面網(wǎng)格如圖4 所示。

圖4 NACA0012 翼型黏性網(wǎng)格Fig.4 Viscous grid of NACA0012 airfoil

該算例中，每隔200 步進(jìn)行一次降階模型求解，快照數(shù)目N=20，保存間隔K=10。求解降階模型時(shí)，只計(jì)算1 000 個(gè)網(wǎng)格點(diǎn)上的殘值。給出了3 種不同的采樣方式確定這1 000 個(gè)網(wǎng)格點(diǎn)，分別是：①隨機(jī)采樣；②選取殘差最大的網(wǎng)格點(diǎn)；③選取流場(chǎng)特征相差最大的一組網(wǎng)格點(diǎn)（通過(guò)計(jì)算不同網(wǎng)格流場(chǎng)值之間的相關(guān)性確定）。得到的流場(chǎng)殘值和阻力系數(shù)CD收斂歷程如圖5 所示。

圖5 不同采樣方式的加速方法收斂歷程Fig.5 Residual history with different sampling methods

圖5 和表2 表明，加速收斂方法對(duì)降階模型中網(wǎng)格點(diǎn)的采樣方式不敏感，不同的采樣方式均能達(dá)到約3 倍的加速效果，但采樣方式2 和3 均會(huì)帶來(lái)額外的計(jì)算量，因此在本文的其他算例中，均只采用隨機(jī)采樣網(wǎng)格的方法。另外，加速方法和初始的CFD 方法相比壓力分布和摩擦阻力系數(shù)Cf分布也基本完全一致，結(jié)果如圖6 所示。

圖6 NACA0012翼型表面壓力分布和摩擦阻力分布對(duì)比Fig.6 Comparison of surface pressure coefficient and skin friction coefficient （NACA0012）

表2 不同采樣方式的加速方法計(jì)算時(shí)間Table 2 Calculation time with different sampling methods

2.3 S809 翼型繞流湍流算例

本節(jié)采用S809 翼型湍流算例用于測(cè)試加速收斂方法在流場(chǎng)中存在大分離時(shí)的加速效果，來(lái)流 馬 赫 數(shù)Ma= 0.2，攻 角α= 14.2°，雷 諾 數(shù)Re= 2 × 10-6，CFL = 5。網(wǎng)格總數(shù)為34 205，近壁面網(wǎng)格如圖7 所示。

圖7 S809 翼型黏性網(wǎng)格Fig.7 Viscous grid of S809 airfoil

該算例中測(cè)試了4 種不同快照參數(shù)（快照數(shù)目N和保存間隔K）對(duì)加速效果的影響，得到的流場(chǎng)殘值收斂歷程如圖8 所示。

圖8 不同快照參數(shù)的收斂歷程Fig.8 Residual history with different snapshot parameters

由圖8 和表3 可知，對(duì)于流場(chǎng)中存在大分離的湍流問(wèn)題，加速收斂方法仍然有效，且加速收斂方法對(duì)降階模型中快照數(shù)目和保存間隔不敏感，不同的參數(shù)均能達(dá)到約3 倍的加速效果。另外，加速方法和初始CFD 方法相比壓力分布和摩擦阻力分布也完全一致（圖9）。同時(shí)，2 種方法計(jì)算得到的壓力云圖和分離區(qū)也一致（圖10）。

表3 不同快照參數(shù)的加速方法計(jì)算時(shí)間Table 3 Calculation time with different snapshot parameters

圖9 S809 翼型表面壓力分布和摩擦阻力分布對(duì)比Fig.9 Comparison of surface pressure coefficient and skin friction coefficient （S809）

圖10 壓力云圖對(duì)比Fig.10 Comparison of pressure contours

2.4 ONERA M6 機(jī)翼繞流湍流算例

本節(jié)將加速收斂應(yīng)用于三維ONERA M6 機(jī)翼亞聲速繞流問(wèn)題，驗(yàn)證該方法對(duì)于三維問(wèn)題的有效性。ONERA M6 機(jī)翼網(wǎng)格如圖11 所示，網(wǎng)格總數(shù)30 萬(wàn)，來(lái)流馬赫數(shù)Ma=0.3，攻角α=8°，雷諾數(shù)Re=1.172×10-7，CFL=1。該算例中網(wǎng)格數(shù)M取5 000，流場(chǎng)快照數(shù)目N=20，快照保存間隔K=20。

圖11 ONERA M6 機(jī)翼網(wǎng)格Fig.11 Grid of ONERA M6 wing

圖12 和表4 分別為該算例的收斂歷程和總用時(shí)，它們表明，對(duì)于復(fù)雜三維問(wèn)題，加速收斂方法仍然有效，能夠?qū)崿F(xiàn)殘值的快速下降。從氣動(dòng)力系數(shù)收斂曲線也可以看到加速收斂方法波動(dòng)更小。圖13 給出了20%展長(zhǎng)處機(jī)翼截面的壓力分布，可以看到2 種方法具有相同的收斂精度。

圖12 收斂歷程（ONERA M6 機(jī)翼）Fig.12 Residual history （ONERA M6 wing）

表4 不同方法的計(jì)算時(shí)間Table 4 Calculation time of different methods

圖13 20%展長(zhǎng)處機(jī)翼截面壓力分布對(duì)比Fig.13 Comparison of pressure coefficient of 20% wing cross-section

3 結(jié) 論

本文提出了一種基于降階模型和梯度優(yōu)化的加速收斂方法。該方法收集CFD 計(jì)算過(guò)程的流場(chǎng)快照，從快照中提取能量最高的POD 模態(tài)而后構(gòu)建基于POD 模態(tài)降階模型。通過(guò)CFD 迭代和降階模型的交替求解實(shí)現(xiàn)加速收斂。特別地，本文采用梯度優(yōu)化最小化降階模型的殘差，相比于其他迭代方法更靈活，且不存在穩(wěn)定性問(wèn)題。在該方法的框架下，降階模型中只是關(guān)于流場(chǎng)快照的簡(jiǎn)單代數(shù)運(yùn)算，并不涉及不同網(wǎng)格之間的相互運(yùn)算。因此，降階模型部分很容易在額外模塊中實(shí)現(xiàn)，不需要對(duì)CFD 求解器源代碼進(jìn)行改動(dòng)，只需讓CFD 求解器在合適位置保存流場(chǎng)快照即可，大大降低了嵌入式降階模型實(shí)現(xiàn)難度。

本文中的算例表明，所提出的加速收斂方法是魯棒的和高效的，對(duì)于無(wú)黏、湍流流動(dòng)均有明顯的加速效果，是一種頗具潛力的加速收斂方法。在進(jìn)一步的工作中，該方法可考慮應(yīng)用于非定常流場(chǎng)的加速求解，使用降階模型預(yù)測(cè)雙時(shí)間步中內(nèi)迭代的初值以加速內(nèi)迭代步收斂。同時(shí)，該方法在較大規(guī)模網(wǎng)格的流場(chǎng)求解中也有潛在價(jià)值，也將在未來(lái)的工作中進(jìn)行測(cè)試和研究。