一種適用于三維混合網(wǎng)格的GMRES加速收斂新方法

張健， 鄧有奇， 李彬， 張耀冰

中國空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所， 綿陽 621000

一種適用于三維混合網(wǎng)格的GMRES加速收斂新方法

張健*， 鄧有奇， 李彬， 張耀冰

中國空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所， 綿陽 621000

為提高流場計(jì)算收斂效率，發(fā)展了一套適用于三維混合網(wǎng)格Naiver-Stokes方程求解的并行廣義最小殘差(GMRES)隱式時(shí)間推進(jìn)方法。該方法由科學(xué)計(jì)算可移植擴(kuò)展工具包(PETSc)中的Krylov子空間求解器實(shí)現(xiàn)，線性方程系統(tǒng)中的系數(shù)矩陣直接以顯式給出以提高算法的穩(wěn)定性。為進(jìn)一步提高GMRES方法的收斂速度，對非結(jié)構(gòu)網(wǎng)格的序號(hào)進(jìn)行了重排序，使得系數(shù)矩陣的非零元素盡量向主對角線靠近。利用所發(fā)展的GMRES方法，完成了對ONERA-M6機(jī)翼、AIAA阻力預(yù)測會(huì)議通用研究模型(CRM)等算例的計(jì)算，計(jì)算結(jié)果與試驗(yàn)結(jié)果吻合良好。通過與其他隱式推進(jìn)方法進(jìn)行比較，對算法的收斂特性進(jìn)行了研究。結(jié)果表明，所發(fā)展的GMRES方法計(jì)算更加穩(wěn)定，殘差下降速度相對LU-SGS(Lower-Upper Symmetric Gauss-Seidel)方法更快，尤其是氣動(dòng)力系數(shù)向著收斂解逼近的速度更加明顯，提高了計(jì)算效率。

混合網(wǎng)格； 收斂效率； 隱式時(shí)間推進(jìn)方法； GMRES； PETSc； 網(wǎng)格重排序

計(jì)算流體力學(xué)(Computational Fluid Dynamics, CFD)技術(shù)在近些年來取得了飛速發(fā)展，這使得求解復(fù)雜的三維流動(dòng)問題更加易于實(shí)現(xiàn)。但隨著工業(yè)界所研究問題的不斷深入，所關(guān)注部件細(xì)節(jié)的更加詳盡，一些工程外形也變得越來越復(fù)雜，從而導(dǎo)致相應(yīng)網(wǎng)格的生成工作變得更加困難，時(shí)間開銷也更大。盡管非結(jié)構(gòu)混合網(wǎng)格技術(shù)在處理復(fù)雜外形方面相對于結(jié)構(gòu)網(wǎng)格有一定優(yōu)勢，但通常為了取得更為精確的計(jì)算結(jié)果，其所需的網(wǎng)格量依然十分巨大。由此生成出一個(gè)大型的非線性耦合方程系統(tǒng)求解問題，需要大量的迭代時(shí)間才能得到一個(gè)收斂的流場。因此，對于解決大型工程問題找到一個(gè)快速有效的時(shí)間推進(jìn)方法十分必要。

常用的時(shí)間推進(jìn)技術(shù)包括顯式的Runge-Kutta法[1]，隱式的LU-SGS(Lower-Upper Symmetric Gauss-Seidel)方法[2]等，其應(yīng)用已較為成熟。隱式策略相對顯式策略能夠選取更大的時(shí)間步長，在很大程度上提高收斂速度?，F(xiàn)今，多重網(wǎng)格技術(shù)[3]作為一種流行的加速收斂技術(shù)，通過在不同疏密程度的網(wǎng)格之間迭代循環(huán)計(jì)算，來消除不同頻次的誤差，以達(dá)到加速收斂的目的。但在一些復(fù)雜情況下，通過全局融合法融合細(xì)網(wǎng)格生成的粗網(wǎng)格可能具有不規(guī)則性，致使疏密網(wǎng)格之間沒有好的插值、延拓算子，導(dǎo)致多重網(wǎng)格技術(shù)存在魯棒性問題。另一種更加復(fù)雜的隱式策略就是Krylov子空間方法，Krylov子空間方法最初用來求解大規(guī)模矩陣問題，在有些情況下能夠達(dá)到超線性甚至二次收斂的效果。Krylov子空間類型的方法有很多種，在CFD應(yīng)用中最廣泛的是廣義最小殘差(Generalized Minimal Residual, GMRES)方法[4], GMRES方法通過Arnoldi過程構(gòu)造Krylov子空間的正交基和求解一個(gè)最小二乘法問題，在Krylov子空間上選擇最優(yōu)解，使每一步子迭代時(shí)的殘差模最小。由于這個(gè)過程本身是面向全頻率的殘差，所以GMRES方法不需要結(jié)合多重網(wǎng)格加速技術(shù)，避免了多重網(wǎng)格不魯棒的問題。GMRES方法在結(jié)構(gòu)網(wǎng)格以及二維非結(jié)構(gòu)CFD代碼中已經(jīng)有過成功的應(yīng)用[5-6]，Xia等[7]將其應(yīng)用在間斷有限元算法(Discontinuous Galerkin, DG)上，燕振國等[8]將其應(yīng)用在高階耗散緊致格式上也都取得了不錯(cuò)的效果。本文的目的旨在將GMRES方法應(yīng)用于復(fù)雜三維混合網(wǎng)格流場求解計(jì)算上，探索出適用于大型三維工程問題的加速收斂技術(shù)，以達(dá)到減少收斂時(shí)間，提高計(jì)算效率的目的。

PETSc(Portable, Extensible Toolkit for Scientific Computation)科學(xué)計(jì)算可移植擴(kuò)展工具包是美國能源部ODE2000支持開發(fā)的20多個(gè)ACTS(Advanced Computational Testing and Simulation)工具箱之一，主要用于在分布式存儲(chǔ)環(huán)境高效求解偏微分方程組及相關(guān)問題[9]。對于大規(guī)模線性方程組，PETSc提供了豐富的Krylov子空間迭代方法和各種預(yù)條件子。PETSc具有高性能、可移植等優(yōu)點(diǎn)，并在CFD 領(lǐng)域的并行應(yīng)用獲得了相當(dāng)?shù)某晒?，其中較為典型的PETSc-FUN3D[10]獲得了1999年的戈登·貝爾特別獎(jiǎng)[11]，國內(nèi)屈崑等[12]將其應(yīng)用在CFD數(shù)學(xué)模型的線性化方法上。本文將在自主研發(fā)的三維混合網(wǎng)格流場求解器MFlow[13]的基礎(chǔ)上，結(jié)合PETSc數(shù)學(xué)工具包實(shí)現(xiàn)GMRES時(shí)間推進(jìn)算法。通過顯式計(jì)算線性系統(tǒng)的系數(shù)矩陣，來提高GMRES算法的穩(wěn)定性。同時(shí)，針對非結(jié)構(gòu)網(wǎng)格編號(hào)順序隨機(jī)的問題，對網(wǎng)格序號(hào)進(jìn)行重排序，得到帶寬減小的系數(shù)矩陣，以提升GMRES的計(jì)算效率。并通過對M6機(jī)翼及AIAA阻力預(yù)測會(huì)議(DPW)標(biāo)模等算例的計(jì)算驗(yàn)證，同時(shí)結(jié)合原有LU-SGS方法計(jì)算結(jié)果的對比分析，完成了對適用于混合網(wǎng)格流場求解加速收斂技術(shù)的研究。

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

三維可壓縮流Navier-Stokes方程表達(dá)式為

(1)

式中：U為守恒變量；FI為無黏通量；FV為黏性通量；Ω為控制單元體積；S為單元面面積。

1.1 隱式時(shí)間離散

對方程式(1)進(jìn)行空間離散后，得到如下半離散形式：

(2)

式中：M為質(zhì)量矩陣(Mass Matrix)，對于基于格心的求解策略，M用單位矩陣；RI為殘差，下標(biāo)I代表對應(yīng)的控制單元。

將式(2)改寫成差分形式得到隱式離散格式：

(3)

式中：Δt為時(shí)間步長；上標(biāo)n為當(dāng)前時(shí)間迭代步；β為0～1之間的系數(shù)，一般取1較為魯棒，取0則退化為顯式格式。

對于隱式求解策略，需要估算出新的時(shí)間步上的殘差：

(4)

式中：?R/?U一般稱為通量雅克比矩陣(Flux Jacobian)。將式(4)代入式(3)，整理簡化得到線性方程組的形式：

(5)

式中：

(6)

1.2 GMRES方法

本文的主要工作就是針對式(5)形式的線性系統(tǒng)進(jìn)行求解，不同的求解方法構(gòu)成了不同的隱式時(shí)間推進(jìn)方法。

(7)

對于大型問題，為減小內(nèi)存開銷，一般采用重啟型GMRES方法。通過在有限的Krylov子空間下反復(fù)重啟迭代GMRES方法，直到殘差收斂到預(yù)定值，算法的具體過程如下所示：

Forl=1,mDom次重啟迭代

r0=P-1v0預(yù)處理

v1=r0/β

Forj=1,kDo 內(nèi)迭代

wj+1=P-1yj+1預(yù)處理

Fori=1,jDo Gram-Schmidt

hi,j=wj+1·vj

wj+1=wj+1-hi,jvj

End Do

hj+1,j=wj+1Hessenberg矩陣

vj=wj+1/hj+1,jKrylov向量

End Do

ΔU0=ΔU重啟

End Do

1.3 隱式系統(tǒng)矩陣

在一般的GMRES方法實(shí)現(xiàn)過程中，雅克比矩陣都是以某一向量內(nèi)積的形式出現(xiàn)，所以矩陣向量乘都以差分近似的無矩陣方法實(shí)現(xiàn)，即

(8)

然而，這個(gè)方法卻存在兩個(gè)問題。首先，小量ε的大小取值依賴于不同的問題，如果流動(dòng)變化劇烈很有可能導(dǎo)致求解的不穩(wěn)定。其次，在每一個(gè)迭代中需要額外計(jì)算一個(gè)殘差向量，增加了計(jì)算量[14]。因此，本文嘗試?yán)蔑@式設(shè)置雅可比矩陣來解決上述問題。

通量雅克比矩陣可以分解為對流通量雅可比矩陣和黏性通量雅可比矩陣，即

(9)

式中：NF為單元面?zhèn)€數(shù)；ΔSIJ為相鄰單元I、J之間單元面面積。黏性通量雅克比計(jì)算較為復(fù)雜，一般利用其譜半徑簡化得到

(10)

式中：ρIJ為密度；γIJ為比熱比；μL和μT分別為層流黏性系數(shù)和湍流黏性系數(shù)；PrL和PrT分別為層流和湍流Prandtl數(shù)。對于無黏通量，根據(jù)Roe通量差分裂格式：

(11)

(12)

(13)

1.4 網(wǎng)格重排序

步驟1從某一個(gè)單元開始(盡量是物面邊界單元)，定義其網(wǎng)格單元所屬超平面序號(hào)Np為1。

步驟2尋找超平面Np所包括網(wǎng)格的相鄰網(wǎng)格, 若相鄰網(wǎng)格還無所屬超平面, 則將其超平面序號(hào)定義為Np+1。

步驟3重復(fù)步驟2直到所有網(wǎng)格都有其對應(yīng)的超平面為止。

步驟4從序號(hào)為1的超平面起按序號(hào)從小到大依次對每個(gè)超平面上的網(wǎng)格進(jìn)行編號(hào), 即得到了排序后的網(wǎng)格編號(hào)。

圖1給出了一個(gè)簡單的二維非結(jié)構(gòu)網(wǎng)格示例圖，右側(cè)給出了其對應(yīng)的系數(shù)矩陣。很明顯，網(wǎng)格編號(hào)的無序性導(dǎo)致矩陣非零元素的分布無特定規(guī)律。圖2給出了經(jīng)過上述重排序算法之后的結(jié)果，可以看到，矩陣的非零元素更加向?qū)蔷€靠近，矩陣的帶寬明顯減小。

圖1 二維非結(jié)構(gòu)網(wǎng)格(左)及相應(yīng)隱式算子矩陣(右)Fig.1 2D unstructured grid(left) and associatedimplict operator matrix(right)

圖2 重排序后的網(wǎng)格及帶寬減小后的矩陣Fig.2 Reordered grids and reduced bandwidth matrix

1.5 結(jié)合PETSc庫并行實(shí)現(xiàn)GMRES方法

PETSc中提供了多個(gè)庫，類似于C++中的類，每個(gè)庫都有相應(yīng)的接口來進(jìn)行數(shù)據(jù)的設(shè)置或者方程的求解，其框架圖[9]如圖3所示。

圖3 PETSc庫框架圖[9]Fig.3 Organization of PETSc libraries[9]

圖4 網(wǎng)格并行劃分示意圖Fig.4 Example of grid’s parallelization

圖5 裝配好的并行塊稀疏矩陣示意圖Fig.5 Assembled sparse parallel matrix in block

2 數(shù)值模擬結(jié)果與分析

本節(jié)通過對ONERA-M6機(jī)翼和AIAA阻力預(yù)測會(huì)議標(biāo)模等算例進(jìn)行計(jì)算，對比了改進(jìn)后的GMRES和其他隱式策略，并研究了多重網(wǎng)格、網(wǎng)格重排序和CFL數(shù)等因素對算法收斂速度的影響。所有算例的計(jì)算都是通過自主開發(fā)的基于混合網(wǎng)格的大型CFD計(jì)算軟件MFlow完成。MFlow軟件基于MPI并行化，使用標(biāo)準(zhǔn)C++語言編寫，采用Roe迎風(fēng)格式[17]計(jì)算通過控制體單元面的無黏通量，通過中心格式計(jì)算黏性項(xiàng)，湍流模型選取Spalart-Allmaras一方程模型[18]。

2.1 ONERA-M6機(jī)翼繞流流場模擬

ONERA-M6機(jī)翼是一個(gè)根弦約為0.8 m，半展長約為1.2 m的后掠機(jī)翼，剖面形狀為對稱翼型。ONERA-M6機(jī)翼繞流問題經(jīng)常用作考核外流流場解算程序的標(biāo)準(zhǔn)算例。

首先選取來流馬赫數(shù)為0.839 5，迎角為3.06° 的計(jì)算狀態(tài)對程序進(jìn)行驗(yàn)證。計(jì)算網(wǎng)格如圖6所示，網(wǎng)格總量大約為80萬，分24核并行計(jì)算完成。

圖6 ONERA-M6機(jī)翼計(jì)算網(wǎng)格Fig.6 Computational grids of ONERA-M6 wing

GMRES方法的重啟次數(shù)采用KSP求解器默認(rèn)的30次，預(yù)處理選取支持并行的Jacobi預(yù)處理器。圖7是本文發(fā)展的方法計(jì)算得到的機(jī)翼上表面等壓力系數(shù)Cp分布圖，由圖中可以很清晰地看到計(jì)算結(jié)果捕捉到了“λ”型激波結(jié)構(gòu)。圖8是機(jī)翼沿展向90%的剖面壓力系數(shù)分布，可以看到GMRES和LU-SGS兩種隱式方法計(jì)算的結(jié)果一致，前緣吸力峰值和激波位置的捕捉十分準(zhǔn)確，與試驗(yàn)數(shù)據(jù)[19]吻合較好，充分驗(yàn)證了新方法的正確性。

為測試不同方法對收斂特性的影響，針對此算例，分別利用LU-SGS和GMRES兩種隱式算法進(jìn)行計(jì)算。理論上，GMRES方法的CFL數(shù)可以取的很大，但由于殘差呈高度非線性，初始幾步迭代可能導(dǎo)致初場遠(yuǎn)遠(yuǎn)偏離穩(wěn)定解[20]。因此本文先在一系列粗網(wǎng)格上利用LU-SGS計(jì)算一定步數(shù)，再將結(jié)果插值到細(xì)網(wǎng)格上得到初始場。這個(gè)過程耗時(shí)較少，對于整個(gè)收斂過程可以忽略不計(jì)，并保證各個(gè)不同算法從相同的初場進(jìn)行啟動(dòng)，便于對比分析。在之后的計(jì)算中，分別應(yīng)用單重網(wǎng)格(SG)LU-SGS，多重網(wǎng)格V循環(huán)LU-SGS，GMRES(Origin)和網(wǎng)格重排序(Reorder)后的GMRES進(jìn)行計(jì)算。

圖7 ONERA-M6機(jī)翼壓力系數(shù)分布云圖Fig.7 Pressure coefficient contours of ONERA-M6

圖8 ONERA-M6機(jī)翼90%展向剖面壓力系數(shù)分布Fig.8 Pressure coefficient distributions at 90% span of ONERA-M6 wing

圖9給出了不同算法隨迭代步數(shù)以及迭代時(shí)間變化的收斂特性曲線。由圖可以明顯看到，殘差下降5個(gè)量級(jí)，GMRES和多重網(wǎng)格V循環(huán)LU-SGS所用步數(shù)差不多，GMRES稍微少一點(diǎn)。由于各個(gè)算法每一個(gè)時(shí)間步耗時(shí)不同，考查收斂效率重點(diǎn)需要看收斂所需時(shí)間?？梢钥吹剑肎MRES和多重網(wǎng)格的LU-SGS耗時(shí)相當(dāng)，比單重LU-SGS殘差下降5個(gè)量級(jí)大概要快1倍左右。而GMRES方法在經(jīng)過了網(wǎng)格重排序后有明顯加速收斂效果，特別是收斂過程的前期殘差的下降要更快一些。

在工程實(shí)踐中，判斷計(jì)算收斂與否主要看氣動(dòng)力是否收斂到穩(wěn)定值。圖10是升力系數(shù)隨計(jì)算時(shí)間的變化曲線，重排序后的GMRES方法計(jì)算升力系數(shù)向穩(wěn)定值逼近速度最快，原始的GMRES方法次之，隨后是多重網(wǎng)格LU-SGS，最慢的是單重的LU-SGS。

圖9 不同推進(jìn)方法殘差收斂特性比較 Fig.9 Comparison of residual convergence property between different methods

圖10 不同方法升力系數(shù)收斂過程比較 Fig.10 Comparison of lift coefficient convergence property between different methods

2.2 CRM標(biāo)模計(jì)算

AIAA阻力預(yù)測會(huì)議的計(jì)算外形是NASA的通用研究模型(Common Research Model，CRM)[21]。該模型由機(jī)身和機(jī)翼組成，其中機(jī)身是典型的商業(yè)運(yùn)輸機(jī)機(jī)身，機(jī)翼是跨聲速超臨界機(jī)翼。該外形在NASA Langley National Transonic Facility和NASA Ames 11-ft風(fēng)洞都進(jìn)行了詳盡的試驗(yàn)[22]。本文所使用的MFlow軟件的計(jì)算結(jié)果也曾參與了第五次阻力預(yù)測會(huì)議(DPW5)[23]的數(shù)據(jù)對比[24]。

DPW提供了多套標(biāo)準(zhǔn)網(wǎng)格，本文選取了DPW4[25]提供的一套三棱柱/四面體非結(jié)構(gòu)網(wǎng)格，以及DPW5的一套中等大小的純六面體網(wǎng)格進(jìn)行計(jì)算，DPW5的外形相比DPW4少了水平尾翼，其表面網(wǎng)格及對稱面網(wǎng)格分別如圖11和圖12 所示。

為了考察隱式算法的收斂特性，對上述兩個(gè)算例進(jìn)行了計(jì)算，計(jì)算狀態(tài)均選取來流馬赫數(shù)為0.85、迎角為2°和雷諾數(shù)為5.0×106，其中DPW4三棱柱/四面體網(wǎng)格總數(shù)大約為2 128萬，分120個(gè)核并行計(jì)算完成，DPW5六面體網(wǎng)格總數(shù)大約為511萬，分48個(gè)核并行計(jì)算完成，兩套網(wǎng)格y+均為1。

圖11 DPW4三棱柱/四面體網(wǎng)格Fig.11 Prism/tetra grid of DPW4

圖12 DPW5六面體網(wǎng)格Fig.12 Hex grid of DPW5

首先對DPW4外形進(jìn)行計(jì)算，之前已經(jīng)考察了網(wǎng)格重排序和多重網(wǎng)格技術(shù)的加速收斂效果，所以這里直接采用重排序后的GMRES和多重網(wǎng)格的LU-SGS方法進(jìn)行對比，同時(shí)還對比了利用式(12)簡化后對角占優(yōu)的系統(tǒng)矩陣(以下簡稱GMRES_D方法)的效果。為考察CFL數(shù)大小對算法收斂性的影響，分別選取CFL=50，100，200進(jìn)行了計(jì)算。圖13是計(jì)算殘差隨計(jì)算時(shí)間下降圖，可以看到GMRES方法整體收斂要比LU-SGS方法快，而且通過增大CFL數(shù)能夠明顯加速收斂。對于GMRES_D方法，殘差下降到10-3后下降速度減緩，可能是因?yàn)橄到y(tǒng)矩陣簡化導(dǎo)致。對于本文所發(fā)展的GMRES方法，CFL數(shù)最多只能取到100，但殘差下降速度依然是最快的。

圖13 不同方法及CFL數(shù)殘差下降曲線(DPW4) Fig.13 Residual histories of different methods and different CFL numbers (DPW4)

圖14給出的是升力系數(shù)收斂曲線對比，通過局部放大圖可以明顯看到采用GMRES方法，尤其是CFL=100時(shí)，升力系數(shù)很快收斂為穩(wěn)定值，而LU-SGS方法計(jì)算的升力系數(shù)依然在劇烈振蕩。

對DPW5外形計(jì)算過程中，同樣采用不同方法和不同CFL數(shù)進(jìn)行了計(jì)算。但是利用多重網(wǎng)格對LU-SGS進(jìn)行加速時(shí)，計(jì)算不穩(wěn)定殘差發(fā)散，所以不得不改用單重網(wǎng)格進(jìn)行計(jì)算。圖15給出了不同算法的殘差收斂曲線。從圖中看到GMRES方法的殘差下降速度明顯比LU-SGS方法快。對于GMRES方法，CFL數(shù)從50提高到100收斂加速明顯。而對于LU-SGS方法，將CFL數(shù)提高到100對殘差下降沒有明顯影響。相同CFL數(shù)下，GMRES方法比GMRES_D方法收斂更快，CFL取100時(shí)殘差振蕩較為明顯，但整體呈快速下降趨勢。

圖14 升力系數(shù)收斂曲線Fig.14 Histories of lift coefficient convergence

圖15 不同算法及CFL數(shù)殘差下降曲線(DPW5) Fig.15 Residual histories of different methods and different CFL numbers (DPW5)

圖16給出了不同算法計(jì)算的阻力系數(shù)的收斂曲線圖，對于GMRES方法，尤其是CFL=100時(shí)，在經(jīng)過初期的振蕩之后很快穩(wěn)定到收斂解，而相同CFL數(shù)條件下用LU-SGS方法計(jì)算至少需要數(shù)倍以上的CPU時(shí)間。

圖16 阻力系數(shù)收斂曲線Fig.16 Histories of drag coefficient convergence

3 結(jié) 論

1) GMRES方法是一個(gè)高效穩(wěn)定的時(shí)間推進(jìn)方法，能夠適應(yīng)不同復(fù)雜外形、不同類型網(wǎng)格的計(jì)算。

2) 相比于多重網(wǎng)格加速后的LU-SGS方法，GMRES方法能夠取得相當(dāng)甚至更快的殘差下降速度，特別是在相同CFL數(shù)下氣動(dòng)力收斂加速效果更加明顯。GMRES方法能夠適應(yīng)大的CFL數(shù)計(jì)算，提高CFL數(shù)能夠明顯提高收斂效率。

3) 對非結(jié)構(gòu)網(wǎng)格編號(hào)進(jìn)行重排序有利于提高GMRES方法的收斂效率。

4) 利用對流通量雅克比矩陣譜半徑對通量雅克比矩陣簡化，保證系統(tǒng)矩陣主對角占優(yōu)能夠增加算法的魯棒性，但可能對于加速收斂效果產(chǎn)生不利影響。

本文所發(fā)展的GMRES加速收斂新方法，通過顯式計(jì)算出通量雅克比矩陣實(shí)現(xiàn)，加速收斂效果明顯，旨在為解決復(fù)雜三維流場計(jì)算中，尤其對于混合網(wǎng)格解算器，多一種時(shí)間推進(jìn)算法選擇，以提高工程效率。雖然GMRES方法在復(fù)雜三維流動(dòng)問題求解上展現(xiàn)出高效、魯棒的特點(diǎn)，但是隱式算法的收斂效率受到多方面因素影響，例如網(wǎng)格質(zhì)量、流動(dòng)特性等。此外，本文程序針對個(gè)別算例的穩(wěn)定性受限于CFL數(shù)，在下一步工作中將針對此作進(jìn)一步研究。

[1] JAMESON A, SCHMIDT W, TURKEK E. Numerical solutions of the euler equations by finite volume methods using Runge-Kutta time-stepping schemes： AIAA-1981-1259[R]. Reston: AIAA, 1981.

[2] YOON S, JAMESON A. Lower-upper symmetric gauss-seidel method for the euler and naiver-stokes equations[J]. AIAA Journal, 1988, 26(9): 1025-1026.

[3] JAMESON A. Multigrid algorithms for compressible flow calculations, multigrid method II[M]. Berlin: Springer Verlag, 2006： 166-201.

[4] SAAD Y, SCHULZ M H. GMRES: A generalized minimum residual algorithm for solving nonsymmetric linear systems[J]. SIAM Journal on Scientific and Statistcal Computing, 1986, 7: 856-869.

[5] WIGTON L B, YOUNG D P, YU N J. GMRES acceleration of computational fluid dynamics codes: AIAA-1985-1494[R]. Reston: AIAA, 1985.

[6] KNOLL D A, MCHUGH P R. Inexact Newton’s method solutions to the incompressible Navier-Stokes and energy equations using standard and matrix-free implementations: AIAA-1993-3332[R]. Reston: AIAA, 1993.

[7] XIA Y D, LUO H, NOURGALIEV R. An implicit hermite WENO resconstruction-based discontinuous Galerkin method on tetrahedral grids[J]. Computers & Fluids, 2014, 96(5): 406-421.

[8] 燕振國, 劉化勇, 毛枚良, 等. 基于高階耗散緊致格式的GMRES方法收斂特性研究[J]. 航空學(xué)報(bào), 2014, 35(5): 1181-1192.

YAN Z G, LIU H Y, MAO M L, et al. Convergence property investigation of GMRES method based on high-order dissipative compact scheme[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(5): 1181-1192 (in Chinese).

[9] BALAY S, ABHYANKAR S, ADAMS M, et al. PETSc users manual[EB/OL]. (2014-06-30)[2015-07-01]. http://www.mcs.anl.gov/petsc/documentation/index.html.

[10] ANDERSON W K, GROPP W D, KAUSHIK D K, et al. Achieving high sustained performance in an unstructured mesh CFD application[C]//Proceedings of the 1999 ACM/IEEE Conference of Supercomputing. New York: IEEE, 1999.

[11] 鐘英, 張理論, 劉巍. PETSc及典型開源流體力學(xué)應(yīng)用[J]. 計(jì)算機(jī)應(yīng)用, 2013, 33(S2): 63-66.

ZHONG Y, ZHANG L L, LIU W. PETSc and its typical open-source computational fluid dynamics application[J]. Journal of Computer Applications, 2013, 33(S2): 63-66 (in Chinese).

[12] 屈崑, 李記超, 蔡晉生. CFD數(shù)學(xué)模型的線性化方法及其應(yīng)用[J]. 航空學(xué)報(bào), 2015, 36(10): 3218-3227.

QU K, LI J C, CAI J S. Method of linearizing computational fluid dynamics model and its applications[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(10): 3218-3227 (in Chinese).

[13] 張耀冰, 鄧有奇, 吳曉軍, 等. DLR-F6翼身組合體數(shù)值計(jì)算[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2011, 29(2): 163-169.

ZHANG Y B, DENG Y Q, WU X J, et al. Drag prediction of DLR-F6 using MFlow unstructured mesh solver[J]. Acta Aerodynamica Sinica, 2011, 29(2):163-169 (in Chinese).

[14] XIA Y D, LUO H, FRISBEY M, et al. A set of parallel, implicit methods for a reconstructed discontinuous Galerkin method for compressible flows on 3D hybrid grids[J]. Computers & Fluids, 2014, 98: 134-151.

[15] BLAZEK J. Computational fluid dynamics principles and applications[M]. Oxford: Elsevier Science Ltd., 2001: 419-420.

[16] 李勁杰, 王剛, 史愛明, 等. 基于非結(jié)構(gòu)網(wǎng)格流場計(jì)算的網(wǎng)格重排序[J]. 航空計(jì)算技術(shù), 2005, 35(3): 25-28.

LI J J, WANG G, SHI A M,et al. Reordering of 3D unstructured grids for computing efficiency improvement[J]. Aeronautical Computer Technique, 2005, 35(3): 25-28 (in Chinese).

[17] LUO H, BAUM J D, LOHNER R. An improved finite volume scheme for compressible flows on unstructured grids: AIAA-1995-0348[R]. Reston: AIAA, 1995.

[18] SPALART P R, ALLMARAS S R. A one equation turbulence model for aerodynamic flows: AIAA-1992-0439[R]. Reston: AIAA, 1992.

[19] SCHMITT V, CHARPIN F. Pressure distributions on the ONERA-M6 wing at transonic mach numbers, experimental data base for computer program assessment: AGARD AR-138-B1[R]. [S.l.]: AGARD, 1979.

[20] 康忠良, 閻超. 適用于混合網(wǎng)格的并行GMRES+LUSGS方法[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2013, 31(2): 225-230.

KANG Z L, YAN C. Parallel GMRES+LU-SGS method for mixed grids[J]. Acta Aerodynamica Sinica, 2013, 31(2): 225-230 (in Chinese).

[21] VASSBERG J C, DEHAAN M A, RIVERS S M, et al. Development of a common research model for applied CFD validation studies: AIAA-2008-6919[R]. Reston: AIAA, 2008.

[22] MELISSA B R, ASHLEY D. Experimental Investigation of the NASA common research model in the NASA langley national transonic facility and NASA Ames 11-ft transonic wind tunnel: AIAA-2011-1126[R]. Reston: AIAA, 2011.

[23] LEVY D W, LAFLIN K R, TINOCO E N, et al. Summary of data from fifth AIAA CFD drag prediction workshop: AIAA-2013-0046[R]. Reston: AIAA, 2013.

[24] GONG X, CHEN J, ZHOU N, et al. The effects of turbulence model corrections on drag prediction of NASA common research model: AIAA-2014-3171[R]. Reston: AIAA, 2014.

[25] VASSBERG J, TINOCO E, MANI M, et al. Summary of date from fourth AIAA CFD drag prediction workshop: AIAA-2010-4547[R]. Reston: AIAA, 2010.

張健男， 碩士， 研究實(shí)習(xí)員。主要研究方向： 計(jì)算空氣動(dòng)力學(xué)。

Tel: 0816-2463271

E-mail: j-zhang@foxmail.com

鄧有奇男， 博士， 研究員， 博士生導(dǎo)師。主要研究方向： 計(jì)算空氣動(dòng)力學(xué)。

Tel: 0816-2463005

E-mail: cai@cardc.cn

李彬男， 博士， 助理研究員。主要研究方向： 計(jì)算空氣動(dòng)力學(xué)。

Tel: 0816-2463276

E-mail: leebin2008@hotmail.com

張耀冰男， 博士， 副研究員。主要研究方向： 計(jì)算空氣動(dòng)力學(xué)。

Tel: 0816-2463276

E-mail: zhyb_super@sina.cn

URL：www.cnki.net/KCMS/detail/11.1929.V.20160302.1449.008.html

*Correspondingauthor.Tel.：0816-2463271E-mail:j-zhang@foxmail.com

AnewmethodtoaccelerateGMRES’sconvergenceapplyingtothree-dimensionalhybridgrid

ZHANGJian*,DENGYouqi,LIBin,ZHANGYaobing

ComputationalAerodynamicInstituteofChinaAerodynamicResearch&DevelopmentCenter,Mianyang621000,China

Toimprovetheconvergenceefficiencyofsolvingtheflowfield,aparallelimplicittimeintegrationmethodgeneralizedminimalresidual(GMRES)isappliedtoathree-dimentionalhybridgridNavier-Stokessolver.ThemethodisimplementedbasedonaKrylovsubspacesolverinthescientificcomputationtoolkitportable,extensibletoolkitforscientificcomputation(PETSc).Thecoefficientmatrixoflinearsystemisprovidedexplicitlytostabilizethescheme.Inordertoacceleratetheconvergencemorespecifically,thecellindexesofunstructuredgridarereorderedsuchthatthesystemmatrix’snonzeroelementsareclusteredclosetothemaindiagonal.ThemethodisappliedtosimulationsofONERA-M6wingandAIAADragPredictionWorkshopmodelCRM.Theresultsshowgreatagreementwithexperimentaldata.Comparisonsaremadebetweendifferentimplicitschemes.TheGMRESmethoddevelopedinthispapershowsmorerobustnessandtheresidual’sconvergencehasasignificantspeedupcomparedwithLU-SGSmethod.Moreover,themethodhasafasterspeedapproximatingtothesteadystateofaerodynamiccoefficient.Itgreatlyimprovesthecomputationalefficiency.

hybridgrid;convergenceefficiency;implicittimeintegrationmethod;GMRES;PETSc;gridreorder

2015-11-25；Revised2015-12-27；Accepted2016-01-24；Publishedonline2016-03-021449

2015-11-25；退修日期2015-12-27；錄用日期2016-01-24； < class="emphasis_bold">網(wǎng)絡(luò)出版時(shí)間

時(shí)間：2016-03-021449

www.cnki.net/KCMS/detail/11.1929.V.20160302.1449.008.html

*

.Tel.：0816-2463271E-mailj-zhang@foxmail.com

張健， 鄧有奇， 李彬， 等． 一種適用于三維混合網(wǎng)格的GMRES加速收斂新方法J． 航空學(xué)報(bào),2016,37(11):3226-3235.ZHANGJ,DENGYQ,LIB,etal.AnewmethodtoaccelerateGMRES’sconvergenceapplyingtothree-dimensionalhybridgridJ.ActaAeronauticaetAstronauticaSinica,2016,37(11):3226-3235.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0038

V211.3; O355

A

1000-6893(2016)11-3226-10