, , , ,

(中國(guó)空氣動(dòng)力研究與發(fā)展中心, 四川 綿陽(yáng) 621000)

0 引 言

間斷Galerkin有限元方法(Discontinuous Galerkin finite element method，DG)[1]是一種適合于非結(jié)構(gòu)網(wǎng)格的高精度格式。Reed & Hill[2]在1973年將其用于求解中子輸運(yùn)方程問(wèn)題。DG方法具有計(jì)算精度高，計(jì)算模板緊致，并行能力強(qiáng)，自適應(yīng)方便等優(yōu)點(diǎn)，適用于非結(jié)構(gòu)網(wǎng)格，適合處理復(fù)雜外形及復(fù)雜邊界。作為一種內(nèi)自由度方法，DG方法通過(guò)提高單元內(nèi)變量多項(xiàng)式的階數(shù)實(shí)現(xiàn)高階精度，其內(nèi)存需求高、計(jì)算量大，隨著空間精度的提高，其內(nèi)存需求量和計(jì)算量非線性增加，因此發(fā)展高階精度DG方法的高效隱式時(shí)間迭代方法有重要的學(xué)術(shù)研究和工程應(yīng)用價(jià)值。

隱式迭代方法需要計(jì)算大型的非線性系統(tǒng)，如何高效求解是其重要研究?jī)?nèi)容。目前多種隱式迭代方法已成功應(yīng)用于間斷Galerkin有限元方法等高階精度方法，如Lower Upper-Symmetric Gauss-Seidel (LU-SGS)[3-6]，GMRES(Generalized Minimal Residual)[7-9]，Implicit Integration Factor(IIF)[10]等。非線性系統(tǒng)的計(jì)算方法和近似處理方法直接影響高階DG的計(jì)算穩(wěn)定性以及計(jì)算效率。

Landmann[11]、楊小權(quán)[12-13]、程劍[14]等在二維非結(jié)構(gòu)網(wǎng)格上通過(guò)精確求解右端項(xiàng)殘值對(duì)變量自由度的雅可比矩陣來(lái)提高計(jì)算的CFL，從而提高DG方法在黏性計(jì)算的計(jì)算效率。

DG方法是一種內(nèi)自由度方法，具有計(jì)算模板緊致的特點(diǎn)，這兩大特性帶來(lái)了如下優(yōu)點(diǎn)：

(1) 單元高斯積分點(diǎn)的變量值來(lái)自單元內(nèi)變量分布多項(xiàng)式，通過(guò)多項(xiàng)式求導(dǎo)能夠獲取變量值對(duì)變量自由度的偏導(dǎo)數(shù)，這使得精確計(jì)算無(wú)黏通量對(duì)變量自由度的雅可比矩陣較為容易。

(2) 單元邊界面上變量梯度只與當(dāng)前單元及其面相鄰單元有關(guān)，無(wú)需重構(gòu)，故能較為簡(jiǎn)單地精確給出變量梯度對(duì)邊界面左右變量自由度的偏導(dǎo)數(shù)。

基于以上優(yōu)點(diǎn)，本文在三維非結(jié)構(gòu)網(wǎng)格基礎(chǔ)上，發(fā)展了求解右端項(xiàng)殘值雅可比矩陣(簡(jiǎn)稱殘值雅可比)的精確計(jì)算方法，同時(shí)建立了簡(jiǎn)化計(jì)算、部分精確計(jì)算方法，實(shí)現(xiàn)雅可比矩陣元素的簡(jiǎn)化和精確計(jì)算，并在此基礎(chǔ)上結(jié)合PETSC(Portable, Extensible Toolkit for Scientific Computation)庫(kù)[15]建立了GMRES隱式迭代方法。采用NACA0012翼型首先研究了每步重啟步數(shù)及殘差收斂參數(shù)對(duì)基于殘值雅可比矩陣精確計(jì)算的GMRES計(jì)算效率影響；然后采用NACA0012翼型、ONERA-M6機(jī)翼比較了基于殘值雅可比矩陣精確計(jì)算的GMRES相比基于矩陣近似計(jì)算的GMRES以及LU-SGS/TVD-RK的計(jì)算效率；最后采用方腔流動(dòng)以及平板層流邊界層驗(yàn)證了本文方法在黏性流數(shù)值模擬中的計(jì)算效率。

1 數(shù)值方法

1.1 控制方程

在直角坐標(biāo)系下，三維非定常可壓縮Euler/Navier-Stokes方程組的守恒形式為：

(1)

當(dāng)ivis=0時(shí)為Euler方程，ivis=1時(shí)為N-S方程。式中U為守恒變量，F(xiàn)(U)=(Fx,Fy,Fz)為無(wú)黏通量張量，F(xiàn)v(U,為黏性通量張量，其具體形式如下：

(2)

(3)

式中，ρ、p分別為密度和壓力；x、y、z分別為笛卡爾坐標(biāo)系下的坐標(biāo)分量；u、v、w分別為笛卡爾坐標(biāo)系下的速度分量；E為總能；τxx、τyy、τzz、Θx、Θy、Θz具體形式見(jiàn)參考文獻(xiàn)[16]。

1.2 DG空間離散

空間離散前首先將整個(gè)求解區(qū)域劃分為互不重疊的非結(jié)構(gòu)單元，對(duì)于三維問(wèn)題，包括三棱柱、四面體、六面體和金字塔單元。

令ivis=1，將方程(1)兩端乘以測(cè)試函數(shù)φi(x,y,z)，然后在單元內(nèi)積分，并使用分部積分方法可得DG的弱形式：

(4)

其中，n為單元邊界面的外法向，?Ω為單元Ω的邊界。假設(shè)變量在單元內(nèi)具有分段的多項(xiàng)式分布：

(5)

RHS=RHSv+RHSf

RHSv=RHSi_v+RHSv_v

RHSf=RHSi_f+RHSv_f

(6)

最后化簡(jiǎn)得：

(7)

其中：

(8)

M為質(zhì)量矩陣，其元素為mij，RHSv、RHSf分別為體積分項(xiàng)殘值和面積分項(xiàng)殘值，RHSi_v、RHSv_v分別為體積分項(xiàng)殘值中的無(wú)黏殘值和黏性殘值，RHSi_f、RHSv_f分別為面積分項(xiàng)殘值中的無(wú)黏殘值和黏性殘值，這四項(xiàng)采用高斯數(shù)值積分計(jì)算，如公式(8)。ω(l)、ωf(k)分別為高斯面積分和體積分加權(quán)系數(shù)，|Jl|、i|Jfk|為坐標(biāo)變化雅可比矩陣行列式。無(wú)黏通量和黏性通量面積分項(xiàng)中的F·n、Fv·n采用數(shù)值通量代替，無(wú)黏數(shù)值通量Hc采用通量差分分裂FDS方法(Flux Difference Splitting)和矢通量分裂FVS方法(Flux Vector Splitting)計(jì)算，黏性數(shù)值通量Hv計(jì)算方法與選取的黏性計(jì)算方法有關(guān)。

1.3 DG時(shí)間離散

隱式迭代方法具有較高的效率，尤其是在高精度計(jì)算時(shí)，目前隱式方法主要有直接求解法、近似因子分解法和迭代法三種。直接求解法在求解網(wǎng)格規(guī)模較大問(wèn)題時(shí)相當(dāng)困難，后兩種方法在近四十年間得到了大量的研究，并在CFD中廣泛應(yīng)用。

對(duì)方程(7)采用一階歐拉后差得到：

(9)

方程(9)的第三項(xiàng)為一大型線性系統(tǒng)Ax=B，以三維DG(P3)為例，矩陣A的維數(shù)為n_cell×(5×20)，其中n_cell為單元總數(shù)，5為方程組個(gè)數(shù)，20為P3次多項(xiàng)式的自由度，直接存儲(chǔ)該矩陣并且求逆，其存儲(chǔ)量和計(jì)算量太大，一般采用迭代法求解該線性系統(tǒng)。

將殘值的面積分項(xiàng)和體積分項(xiàng)帶入方程(9)中第二式的左端RHS得：

(10)

將殘值項(xiàng)帶入方程(10)，對(duì)其采用鏈?zhǔn)椒▌t得：

=RHSn

(11)

在迭代求解Ax=B過(guò)程中如果精確給出矩陣A中的矩陣元素，方程兩端的匹配性和相容性更好，這將提高求解過(guò)程的CFL數(shù)，增強(qiáng)求解過(guò)程穩(wěn)定性，提高計(jì)算效率。

1.3.1 右端項(xiàng)殘差雅可比矩陣精確求解方法

在精確求解右端項(xiàng)殘差雅可比矩陣時(shí)，需確保殘值雅可比的無(wú)黏通量和黏性通量完全與右端項(xiàng)殘值中的通量一致，因此本文針對(duì)Roe[17]格式以及BR2[18]黏性計(jì)算格式建立對(duì)應(yīng)的雅可比矩陣精確求解方法。

文獻(xiàn)[19]給出了Roe格式的具體表達(dá)形式，由方程(8)的第三式及方程(11)的第二式，通過(guò)鏈?zhǔn)椒▌t求導(dǎo)，

(12)

(13)

(14)

(15)

(16)

公式(16)中各符號(hào)及具體表達(dá)式可參看文獻(xiàn)[19]。

(1) 首先計(jì)算出Roe格式中邊界面左右的原始變量ρ,u,v,w,p以及V,c,H對(duì)邊界面左右守恒變量偏導(dǎo)數(shù)；

(3) 其次再計(jì)算出a1,a2,a3,a4,a5,a6,a7,a8對(duì)邊界面左右守恒變量的偏導(dǎo)數(shù)；

黏性通量的雅可比矩陣由于包含了u,v,w,T的一階導(dǎo)數(shù)，求解過(guò)程相對(duì)更復(fù)雜。BR2格式的黏性通量計(jì)算表達(dá)式如下：

Hv(UL,UR,UL,UR,n)=

Fv(UR,UR+ηerRf))·n

(17)

對(duì)于當(dāng)前單元L，進(jìn)一步寫(xiě)為：

(18)

(19)

(2) 計(jì)算守恒變量及其方向?qū)?shù)對(duì)其守恒變量自由度的偏導(dǎo)數(shù)。

(20)

(3) 計(jì)算u,v,w,T對(duì)守恒變量的偏導(dǎo)數(shù)。

(21)

(4) 結(jié)合(19)、(20)、(21)，計(jì)算u,v,w,T的一階導(dǎo)數(shù)對(duì)守恒變量自由度的偏導(dǎo)數(shù)。

(22)

(23)

(5) 結(jié)合(23)計(jì)算出層流黏性系數(shù)μL對(duì)守恒變量自由度的偏導(dǎo)數(shù)。

(24)

(8) 最后按照公式(18)求出黏性通量的雅可比矩陣。

同理求出當(dāng)前單元相鄰一側(cè)單元的黏性通量的面積分項(xiàng)的雅可比矩陣。對(duì)于邊界單元的邊界面，將其按內(nèi)部邊界來(lái)獲得無(wú)黏通量雅可比矩陣，并由邊界黏性通量獲得黏性雅可比矩陣。

(25)

1.3.2 右端項(xiàng)殘差雅可比矩陣近似求解方法

為提高計(jì)算效率，方程(11)左端項(xiàng)中的RHSi_f采用簡(jiǎn)單的LF數(shù)值通量格式。

(26)

在有限體積方法中，殘差中黏性通量的雅可比矩陣采用黏性譜半徑簡(jiǎn)化代替，本文同樣采用類似的方法來(lái)處理。

(27)

(28)

其中Ω為當(dāng)前單元的體積，Ss為當(dāng)前單元邊界面S的面積，nface為當(dāng)前單元體總的面數(shù),γ為比熱比，Pr為普朗特?cái)?shù)。

(29)

1.3.3 GMRES方法

GMRES算法是以Galerkin原理為基礎(chǔ)的Krylov子空間投影法，通過(guò)Arnoldi過(guò)程構(gòu)造Krylov子空間的正交基，同時(shí)求解一個(gè)最小二乘法問(wèn)題在Krykov子空間上選擇最優(yōu)解，使得每一步子迭代時(shí)的殘值模最小。理想狀態(tài)下GMRES具有平方收斂速度，計(jì)算效率高。該方法在基于結(jié)構(gòu)網(wǎng)格和非結(jié)構(gòu)網(wǎng)格的有限體積和有限差分方法中得到大量應(yīng)用。線性系統(tǒng)的收斂速度與左端項(xiàng)系數(shù)矩陣的條件數(shù)有關(guān)，其矩陣條件數(shù)越小，收斂速度越快，一般采用預(yù)處理方法來(lái)改善系數(shù)矩陣條件數(shù)。本文采用PETSc科學(xué)計(jì)算可移植擴(kuò)展工具包調(diào)用GMRES以及預(yù)處理技術(shù)。

Forl=1,mDom次重啟迭代

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

b=‖r0‖2

v1=r0/b

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

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

Fori=1,jDo Gram-Schmidt

hi,j=wj+1×vi

wj+1=wj+1-hi,jvi

End Do

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

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

End Do

DU0=DU重啟

End Do

(30)

對(duì)于大型問(wèn)題，為減小內(nèi)存開(kāi)銷，一般采用帶預(yù)處理的重啟型GMRES算法。本文采用的預(yù)處理方法為ILU0(incomplete LU factorizations with zero fill-in)，預(yù)處理矩陣為左端項(xiàng)系數(shù)矩陣，算法的具體過(guò)程如上所示。

2 數(shù)值算例分析

為驗(yàn)證本文發(fā)展的基于殘值雅可比矩陣精確求解方法的GMRES的計(jì)算效率，本節(jié)采用二維、三維算例，通過(guò)比較殘差雅可比精確求解方法、殘值雅可比近似求解方法以及LU-SGS的收斂速度來(lái)驗(yàn)證殘值無(wú)黏項(xiàng)雅可比矩陣精確求解方法和黏性項(xiàng)雅可比矩陣精確求解方法的計(jì)算效率。驗(yàn)證算例包括：NACA0012亞聲速無(wú)黏繞流、ONERA-M6機(jī)翼亞聲速無(wú)黏繞流、方腔流動(dòng)以及平板層流流動(dòng)。整個(gè)測(cè)試在中國(guó)空氣動(dòng)力研究與發(fā)展中心的計(jì)算集群上開(kāi)展，計(jì)算節(jié)點(diǎn)為Intel(R) Xeon(R) E5-2660 V3，主頻2.60 GHz。為方便描述，以GMRES-Ex、GMRES-Ap、GMRES-Ap1表示迭代方法為GMRES且殘值雅可比分別為完全精確求解、完全近似求解以及雅可比中無(wú)黏項(xiàng)精確求解和黏性項(xiàng)近似求解；LU-SGS表示迭代方法為L(zhǎng)U-SGS且殘值雅可比都近似處理；RK表示顯式時(shí)間迭代。

2.1 NACA0012亞聲速無(wú)黏繞流

本小節(jié)以NACA0012翼型為例，首先考察本文的GMRES方法的重啟次數(shù)及每步收斂判定系數(shù)對(duì)計(jì)算效率的影響，在此基礎(chǔ)上針對(duì)本算例選擇最佳的重啟次數(shù)和每步收斂判定系數(shù)；然后對(duì)比分析不同隱式迭代方法、殘值雅可比不同計(jì)算方法及顯式時(shí)間迭代的計(jì)算效率。計(jì)算網(wǎng)格如圖1，物面共有398個(gè)單元，對(duì)稱面共7020個(gè)三角形單元，通過(guò)展向拉伸對(duì)稱面單元得到三棱柱單元。計(jì)算的來(lái)流馬赫數(shù)Ma=0.4，來(lái)流迎角0°，計(jì)算采用12核。分析GMRES參數(shù)影響時(shí)DG方法的計(jì)算精度為二階。

圖1 NACA0012物面附近計(jì)算網(wǎng)格Fig.1 Mesh near NACA0012 surface

圖2給出了GMRES每步不同重啟次數(shù)(Restarted iteration)對(duì)密度殘值收斂速度影響，此時(shí)每步收斂判定系數(shù)emax=0.1。從殘值隨CPU計(jì)算時(shí)間收斂曲線看到，Restarted iteration=5的收斂速度最慢，當(dāng)Restarted iteration≥10后不同Restarted iteration對(duì)收斂速度基本沒(méi)影響。

圖2 GMRES不同重啟次數(shù)的殘值收斂曲線Fig.2 Residual convergence history for different restarted iterations of GMRES

圖3給出了Restarted iteration=5、15時(shí)GMRES每步收斂判定系數(shù)對(duì)密度殘差收斂速度影響?？梢钥吹?，不同判定系數(shù)對(duì)收斂速度有較大影響，隨著系數(shù)增大，收斂速度明顯增大，emax=0.1收斂的總時(shí)間約為emax=0.001的1/2。同時(shí)看到在相同emax時(shí)，Restarted iteration=15的收斂速度都快于Restarted iteration=5。故在后面的GMRES隱式迭代中，取Restarted iteration=15，emax=0.1。

圖3 GMRES不同收斂判定系數(shù)的殘值收斂曲線Fig.3 Residual convergence history for different convergence coefficients of GMRES

圖4給出了DG(P1)在CFL=1000時(shí)不同迭代方法的收斂曲線，CFL在1000步以內(nèi)從1逐漸增大到1000。GMRES_Ex、GMRES_Ap和LU-SGS三種方法的殘值收斂所需的迭代步數(shù)遠(yuǎn)小于顯式迭代，GMRES_Ex隱式迭代所需的迭代步數(shù)最小。相同CFL下GMRES_Ex殘值收斂的迭代步數(shù)約只有GMRES_Ap的1/5，GMRES_Ap殘值收斂的迭代步數(shù)只有LU-SGS的1/3。三種隱式時(shí)間迭代的單步計(jì)算時(shí)間GMRES_Ex > GMRES_Ap > LU-SGS，從殘值隨計(jì)算時(shí)間曲線圖5看，GMRES_Ex隱式迭代所需的計(jì)算時(shí)間最小。GMRES_Ex殘值收斂的計(jì)算時(shí)間約為GMRES_Ap的1/4，GMRES_Ap殘值收斂的計(jì)算時(shí)間只有LU-SGS的1/2。相比LU-SGS，基于殘值雅可比精確求解方法的GMRES計(jì)算效率提高了約8倍。

圖4 CFL=1000，不同迭代方法殘值收斂曲線(殘差與迭代步)Fig.4 Residual convergence history for different iteration methods (Res. vs. step)

圖5 CFL=1000，不同時(shí)間迭代方法殘值收斂曲線(殘差與時(shí)間)Fig.5 Residual convergence history for different iteration methods (Res. vs. time)

圖6給出了不同隱式時(shí)間迭代方法得到的升力系數(shù)隨計(jì)算時(shí)間變化曲線。由圖看到，三種隱式方法得到的升力系數(shù)相同，GMRES計(jì)算得到的升力系數(shù)經(jīng)過(guò)前期短時(shí)間的震蕩后幾乎單調(diào)收斂，其氣動(dòng)力收斂速度遠(yuǎn)高于LU-SGS。

圖6 二階精度下不同方法的升力系數(shù)收斂曲線Fig.6 Lift coefficient convergence history for different iteration methods for DG(P1)

圖7 不同計(jì)算精度下，不同迭代方法殘值收斂曲線Fig.7 Residual convergence history for different iteration method for different orders of accuracy (Res. vs. time)

圖7給出了不同精度下，不同隱式離散方法的密度殘值的收斂曲線。計(jì)算的最大CFL數(shù)取100，CFL在1000步內(nèi)從1逐漸增大到100。從圖看到，隨著計(jì)算精度的提高，收斂所需的時(shí)間急劇增大，不同計(jì)算精度下，GMRES_Ex的計(jì)算效率都最高。

表1給出了不同計(jì)算精度下，不同隱式時(shí)間離散方法的計(jì)算時(shí)間，其中耗時(shí)比以GMRES_Ex計(jì)算時(shí)間為基準(zhǔn)。首先分析相同雅可比矩陣計(jì)算方法下不同計(jì)算精度的收斂時(shí)間，從表看到隨著DG計(jì)算精度的提高，收斂所需計(jì)算時(shí)間非線性增加，以GMRES_Ex為例，DG(P2)和DG(P3)的計(jì)算時(shí)間分別約為DG(P1)的8.9倍和74.3倍。其次分析相同計(jì)算精度下，不同隱式時(shí)間迭代方法的計(jì)算時(shí)間，對(duì)于DG(P1)、DG(P2)、DG(P3)，GMRES_Ap的計(jì)算時(shí)間分別約為GMRES_Ex的3.0倍、5.4倍、6.3倍，LU-SGS的計(jì)算時(shí)間分別約為GMRES_Ex的5.8倍、5.6倍、7.0倍。從以上分析看到，相比殘值雅可比近似求解方法，本文發(fā)展的殘值雅可比精確求解技術(shù)能夠顯著提高GMRES的計(jì)算效率，這對(duì)計(jì)算量巨大的DG方法來(lái)說(shuō)具有重要意義。

表1 不同精度以及不同隱式迭代方法下NACA0012的計(jì)算時(shí)間Table 1 Computation time of different iteration methods with different orders of accuracy for NACA0012 airfoil

2.2 ONERA-M6機(jī)翼亞聲速繞流

為進(jìn)一步驗(yàn)證本文發(fā)展的基于殘值雅可比矩陣精確求解方法的GMRES方法在三維外形的計(jì)算效率，本節(jié)開(kāi)展了ONERA-M6機(jī)翼[20]的亞聲速無(wú)黏繞流流場(chǎng)的數(shù)值模擬。計(jì)算條件為Ma=0.4,α=0°。網(wǎng)格如圖8，網(wǎng)格共約14.8萬(wàn)個(gè)四面體單元，為較好模擬機(jī)翼前后緣，機(jī)翼前后緣及空間分別采用各向異性和各向同性四面體網(wǎng)格，采用64核并行計(jì)算。

圖9～圖12分別給出了DG(P1)、DG(P3)下密度殘值的收斂曲線，同時(shí)給出了RKDG顯式收斂曲線。不同計(jì)算精度下GMRES_Ex隱式計(jì)算的CFL=1000，GMRES_Ap和LU-SGS的CFL數(shù)采用能夠穩(wěn)定計(jì)算的最大值，P1階時(shí)CFL=100，P3階時(shí)CFL=20，RK的CFL=0.3。以殘差下降到10-12量級(jí)為標(biāo)準(zhǔn)，對(duì)于DG(P1)、DG(P3)，GMRES_Ap所需要的計(jì)算時(shí)間約為GMRES_Ex的5倍和10倍，LU-SGS所需要的計(jì)算時(shí)間約為GMRES_Ex的16倍和14倍。

圖8 ONERA-M6機(jī)翼計(jì)算網(wǎng)格Fig.8 Tetrahedral mesh for ONERA-M6

圖9 M6機(jī)翼DG(P1)的密度殘值收斂曲線(殘差與迭代步)Fig.9 Residual convergence history of DG(P1) for M6 wing(Res. vs. step)

圖10 M6機(jī)翼DG(P1)的密度殘值收斂曲線(殘差與時(shí)間)Fig.10 Residual convergence history of DG(P1) for M6 wing(Res. vs. time)

圖11 M6機(jī)翼DG(P3)的密度殘值收斂曲線(殘差與迭代步)Fig.11 Residual convergence history of DG(P3) for M6 wing (Res. vs. step)

圖13給出了不同精度DG方法采用GMRES_Ex獲得的升力系數(shù)隨迭代步數(shù)變化曲線。從圖看到，不同精度下的升力系數(shù)都能500～1000迭代步內(nèi)收斂，顯示了本文的基于殘值雅可比精確求解的GMRES方法在計(jì)算效率方面的優(yōu)勢(shì)。

表2給出了不同計(jì)算精度下，不同隱式方法模擬ONERA-M6機(jī)翼無(wú)黏繞流的時(shí)間，其中耗時(shí)比以GMRES_Ex計(jì)算時(shí)間為基準(zhǔn)。從表可知基于殘值雅可比精確求解的GMRES計(jì)算效率遠(yuǎn)高于雅可比矩陣近似求解的GMRES和LU-SGS，在DG(P3)時(shí)，GMRES_Ex的計(jì)算效率相比另外兩種方法提高了一個(gè)數(shù)量級(jí)以上。

圖12 M6機(jī)翼DG(P3)的密度殘值收斂曲線(殘差與時(shí)間)Fig.12 Residual convergence history of DG(P3) for M6 wing(Res. vs. time)

圖13 DG不同精度的M6機(jī)翼升力系數(shù)收斂曲線Fig.13 Lift coefficient convergence history of DG with different orders of accuracy for ONERA-M6

計(jì)算精度迭代方法收斂判據(jù)計(jì)算時(shí)間(s)耗時(shí)比DG(P1)GMRES_Ex1.0×10-116371GMRES_Ap1.0×10-1128524.5LU-SGS1.0×10-11947314.9DG(P3)GMRES_Ex1.1×10-11198861GMRES_Ap1.1×10-111895939.5LU-SGS1.1×10-1128866214.5

2.3 方腔流動(dòng)

方腔流動(dòng)是一個(gè)經(jīng)典層流標(biāo)準(zhǔn)算例。它描述的是一類頂蓋驅(qū)動(dòng)流動(dòng)，頂蓋以給定速度u=1 m/s驅(qū)動(dòng)方腔內(nèi)流動(dòng)，不同雷諾數(shù)下方腔內(nèi)有不同的流態(tài)。本節(jié)采用該算例驗(yàn)證建立的殘值雅可比精確求解方法，分析精確計(jì)算殘值中無(wú)黏項(xiàng)和黏性項(xiàng)的雅可比矩陣對(duì)計(jì)算效率影響。

圖14、圖15給出了不同隱式迭代方法及殘值雅可比不同求解方法下DG(P1)的密度殘值收斂曲線，GMRES-Ex的CFL數(shù)在200步內(nèi)從1增加到1000。LU-SGS、GMRES-Ap、GMRES-Ap1的CFL數(shù)均采用能夠穩(wěn)定計(jì)算的最大值，CFL數(shù)在1000步內(nèi)從1增加到100。數(shù)值模擬發(fā)現(xiàn)，對(duì)于LU-SGS、GMRES-Ap、GMRES-Ap1這三種方法，當(dāng)CFL數(shù)大于100后CFL數(shù)對(duì)收斂速度影響很小，因此采用CFL數(shù)=100來(lái)比較計(jì)算效率是合適的。由于LU-SGS及GMRES-Ap收斂需要較長(zhǎng)時(shí)間，本文并未給出其完全收斂的收斂曲線。從圖看到GMRES-Ex及GMRES-Ap1的密度殘值很快下降到10-17次方量級(jí)，GMRES-Ex的收斂時(shí)間最短，GMRES-Ap1的收斂時(shí)間次之，LU-SGS的收斂時(shí)間最長(zhǎng)。

圖14 方腔DG(P1)的密度殘值收斂曲線(殘差與迭代步)Fig.14 Residual convergence history of DG(P1) for square cavity(Res. vs. step)

圖15 方腔DG(P1)的密度殘值收斂曲線(殘差與時(shí)間)Fig.15 Residual convergence history of DG(P1) for square cavity(Res. vs. time)

表3統(tǒng)計(jì)了二階計(jì)算精度下，不同隱式迭代方法及殘值雅可比不同求解方法的模擬時(shí)間。以GMRES-Ex模擬時(shí)間為基準(zhǔn)，當(dāng)密度殘值都下降到2.12×10-13時(shí)，GMRES-Ap1、GMRES-Ap、LU-SGS的計(jì)算時(shí)間分別是GMRES-Ex的3.4倍、12.3倍、37.4倍。從表可知GMRES-Ex的計(jì)算效率遠(yuǎn)高于LU-SGS，體現(xiàn)了本文發(fā)展的殘值雅可比精確求解方法的計(jì)算優(yōu)勢(shì)，同時(shí)看到即使只對(duì)殘值中無(wú)黏項(xiàng)的雅可比矩陣精確求解同樣能夠提高計(jì)算效率，精確求解黏性項(xiàng)雅可比矩陣能夠帶來(lái)計(jì)算效率提升。

表3 DG(P1)不同隱式時(shí)間離散方法的方腔計(jì)算時(shí)間Table 3 Computation time of DG(P1) with different iteration methods for square cavity

圖16、圖17給出了不同隱式迭代方法及殘值雅可比不同求解方法下DG(P2)的密度殘值收斂曲線，GMRES-Ex的CFL數(shù)在200步內(nèi)從1增加到1000，GMRES-Ap、GMRES-Ap1、LU-SGS的CFL數(shù)在5000步內(nèi)從1增加到100。由于GMRES-Ex殘值雅可比精確求解，保證迭代求解過(guò)程中能取較大的CFL數(shù)，且基本不受計(jì)算精度影響，而基于雅可比近似求解的GMRES-Ap、GMRES-Ap1的CFL數(shù)受計(jì)算精度影響較大。

圖16 方腔DG(P2)的密度殘值收斂曲線(殘差與迭代步)Fig.16 Residual convergence history of DG(P2) for square cavity(Res. vs. step)

2.4 二維平直平板層流邊界層

最后本文將殘值雅可比精確求解方法用于二維平板層流邊界層的數(shù)值模擬。平板為厚度為0的絕熱壁，來(lái)流馬赫數(shù)為0.2，來(lái)流雷諾數(shù)為105，來(lái)流溫度為288.15 K。計(jì)算區(qū)域?yàn)閤=(-0.5，1)，y=(0，1)，平板前緣點(diǎn)位于(0，0)，后緣點(diǎn)位于(1，0)。網(wǎng)格點(diǎn)為61×17，平板前方分布20個(gè)單元，平板表面分布40個(gè)單元，y方向共布置16個(gè)單元，展向一個(gè)網(wǎng)格單元。平板前后緣X方向的尺寸為0.001，0.11。第一層網(wǎng)格高度為0.0005。計(jì)算包括兩套網(wǎng)格，分別為六面體網(wǎng)格以及由六面體剖分而得的四面體網(wǎng)格。圖18、圖19展示了計(jì)算網(wǎng)格。

圖17 方腔DG(P2)的密度殘值收斂曲線(殘差與時(shí)間)Fig.17 Residual convergence history of DG(P2) for square cavity(Res. vs. time)

圖18 平板層流邊界層計(jì)算網(wǎng)格(六面體)Fig.18 Hexahedron mesh for laminar flow over flat plate

圖19 平板層流邊界層計(jì)算網(wǎng)格(四面體)Fig.19 Tetrahedral mesh for flat laminar flow over flat plate

表4給出了兩種網(wǎng)格在不同計(jì)算精度下得到的總摩阻系數(shù)，同時(shí)給出了Blasius解。從表可知，隨著計(jì)算精度增大，總的摩阻系數(shù)增大。對(duì)于六面體和四面體網(wǎng)格，隨著計(jì)算精度增加，摩阻系數(shù)逐步向Blasius解靠近。

圖20、圖21給出了采用不同精度DG方法的密度殘值收斂歷程，計(jì)算網(wǎng)格為六面體網(wǎng)格。不同計(jì)算精度都采用GMRES迭代，精確求解方程殘值雅可比，CFL數(shù)從1增加到104。從圖看到本文的數(shù)值方法在200步以內(nèi)將殘值降到10-15量級(jí)，進(jìn)一步驗(yàn)證了本文的殘差雅可比精確求解方法的準(zhǔn)確性以及隱式迭代方法的計(jì)算效率。

表4 不同網(wǎng)格和精度下平板的摩擦阻力系數(shù)Table 4 Skin drag coefficient of flat with different grids and orders of accuracy of DG

圖20 平板邊界層六面體網(wǎng)格收斂歷程(殘差與迭代步)Fig.20 Residual convergence history of DG with different orders of accuracy for flat plate(Res. vs. time)

圖21 平板邊界層六面體網(wǎng)格收斂歷程(殘差與時(shí)間)Fig.21 Residual convergence history of DG with different orders of accuracy for flat plate(Res. vs. time)

圖22 六面體網(wǎng)格下不同計(jì)算精度的平板摩阻系數(shù)分布Fig.22 Comparison of the skin friction calculated from DG with different orders of accuracy for flat plate on hexahedral mesh

圖22給出了不同精度DG方法得到的平板的摩阻分布，同時(shí)也給出了Blasius摩阻分布，X、Y坐標(biāo)都為對(duì)數(shù)坐標(biāo)。總的來(lái)看，除去平板前緣，不同精度的摩阻分布與Blasius解分布基本重合，計(jì)算精度越高，重合度越高。高精度的優(yōu)勢(shì)主要體現(xiàn)在平板前緣附近的摩阻系數(shù)。從圖看到，DG(P1)在x<10-2范圍內(nèi)與Blasius解有較大差別，且x越小，差別越大，DG(P2)、DG(P3)在x<10-3范圍內(nèi)才與Blasius解出現(xiàn)明顯區(qū)別。

3 結(jié) 論

本文在三維非結(jié)構(gòu)網(wǎng)格上建立了高階精度DG方法，針對(duì)高階精度DG方法計(jì)算量大這一難題，發(fā)展了一種基于右端項(xiàng)殘值雅可比矩陣精確計(jì)算的高效GMRES隱式迭代方法。采用NACA0012、ONERA-M6、方腔流動(dòng)、平板層流邊界層算例，通過(guò)對(duì)比基于殘值雅可比矩陣不同計(jì)算方法的GMRES、LU-SGS隱式時(shí)間迭代以及顯式TVD-RK的計(jì)算效率，體現(xiàn)了不同殘值雅可比矩陣計(jì)算方法對(duì)計(jì)算效率的影響以及GMRES和LU-SGS的計(jì)算效率。從對(duì)比結(jié)果看到：右端項(xiàng)殘值雅可比矩陣精確求解方法更好保證方程組左右兩端的匹配性和相容性，基于矩陣精確求解方法的GMRES能夠顯著提高不同精度下DG方法數(shù)值模擬的CFL數(shù)，大幅提高計(jì)算效率。

致謝:感謝上海大學(xué)楊小權(quán)對(duì)本文雅可比矩陣精確求解部分提供的幫助！