曹文博，劉溢浪，張偉偉

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

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

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

早期的外插方法幾乎都是基于線性迭代假設xj+1=Axj+b（其中A、b 未知），通過迭代過程中的解快照{ xn}近似該系統(tǒng)解s：s=As+b。這些方法主要分為2 大類：①epsilon 算法；②多項式方法。向量epsilon 算法[5］（Vector Epsilon Algorithm，VEA）和拓撲epsilon 算法[6］（Topological Epsilon Algorithm，TEA）是最早的外插方法，其借助求和方法和一些復雜的反演公式將緩慢收斂的序列轉化為快速收斂的序列，但這類方法的復雜性使得其應用較困難。多項式是一種更簡單的外插方法，最小多項式外插[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ā)展并應用于Euler 方程和其他迭代求解器的加速求解[11，14］。Djeddi 等[15］假設迭代線性收斂，將方程殘值表示為解的線性函數R(U )≈AU-b，并結合降階模型近似該系統(tǒng)的解，實現了雷諾平均Navier-Stokes 方程（Reynlods-Averaged Navier-Stokes，RANS）的加速收斂。

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

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

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

1 數值方法

1.1 控制方程

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

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

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

對于理想氣體，有

由于SA 方程的復雜性，本文不對SA 方程做降階處理，即求解降階模型時渦黏固定不變。

1.2 數值格式

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

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

1.3 基于POD 的降階模型

降階模型用于將控制方程轉化為給定POD模態(tài)下的低階常微分方程組，從而極大地減少了計算時間。首先收集CFD 計算得到的流場快照X=[U1，U2，…，UN]，對其進行奇異值分解得到POD 模態(tài)（通常對奇異值較小的模態(tài)截斷）。

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

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

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

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

梯度下降是一種用來最小化目標函數R(ξ)的方法，根據目標函數對參數ξ 的負梯度-?R(ξ)以給定學習率ε 更新參數的值。

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

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

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

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

2 典型算例驗證

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

2.1 NACA0012 翼型繞流無黏算例

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

圖1 NACA0012 翼型無黏網格Fig.1 Inviscid grid of NACA0012 airfoil

圖2 不同網格點數的加速方法收斂歷程Fig.2 Residual history with different numbers of grids

表1 不同網格點數的加速方法計算時間Table 1 Calculation time with different numbers of grids

由圖2 和表1 可知，所提的加速收斂方法對于無黏流動能達到約1.5 倍的加速收斂效果，從收斂歷程來看，該方法對于網格數的選取不敏感，但由于降階模型中計算的網格數越多，降階模型計算時間越長，因此網格數M不宜太大。圖3 給出了加速方法和初始CFD 方法的壓力系數Cp分布對比，可看出2 種方法計算得到的翼型表面壓力分布一致。事實上，因為降階模型中的快照來源于CFD 計算并且被不斷更新，所以加速收斂方法的精度仍然取決于CFD 計算過程，不會損失任何精度。

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

2.2 NACA0012 翼型繞流湍流算例

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

圖4 NACA0012 翼型黏性網格Fig.4 Viscous grid of NACA0012 airfoil

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

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

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

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

表2 不同采樣方式的加速方法計算時間Table 2 Calculation time with different sampling methods

2.3 S809 翼型繞流湍流算例

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

圖7 S809 翼型黏性網格Fig.7 Viscous grid of S809 airfoil

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

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

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

表3 不同快照參數的加速方法計算時間Table 3 Calculation time with different snapshot parameters

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

圖10 壓力云圖對比Fig.10 Comparison of pressure contours

2.4 ONERA M6 機翼繞流湍流算例

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

圖11 ONERA M6 機翼網格Fig.11 Grid of ONERA M6 wing

圖12 和表4 分別為該算例的收斂歷程和總用時，它們表明，對于復雜三維問題，加速收斂方法仍然有效，能夠實現殘值的快速下降。從氣動力系數收斂曲線也可以看到加速收斂方法波動更小。圖13 給出了20%展長處機翼截面的壓力分布，可以看到2 種方法具有相同的收斂精度。

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

表4 不同方法的計算時間Table 4 Calculation time of different methods

圖13 20%展長處機翼截面壓力分布對比Fig.13 Comparison of pressure coefficient of 20% wing cross-section

3 結 論

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

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