于要杰 劉 鋒 高 超 馮 毅,2)

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

?(加州大學(xué)爾灣分校機(jī)械與航天工程系,美國(guó)加利福尼亞州 92697-3975)

**(成都流體動(dòng)力創(chuàng)新中心,成都 610071)

引言

高階精度數(shù)值格式因其高分辨率、低色散和低耗散等特點(diǎn)正日益受到關(guān)注[1-2].相比目前工業(yè)界最常用的二階精度格式,高階精度格式被認(rèn)為在花費(fèi)相同計(jì)算代價(jià)的前提下具備獲得更高精度數(shù)值結(jié)果的巨大潛力.當(dāng)前最常見的高階精度數(shù)值格式主要包括基于結(jié)構(gòu)網(wǎng)格的高階有限差分方法[3-4]、非結(jié)構(gòu)網(wǎng)格上的有限體積法[5-6]以及基于局部單元近似的間斷有限元(discontinuous Galerkin,DG)[7-8]、譜體積(spectral volume,SV)[9]和譜差分(spectral difference,SD)[10]等方法.

最近,一種基于非結(jié)構(gòu)網(wǎng)格的高階精度通量重構(gòu)格式(flux reconstruction,FR)因其構(gòu)造簡(jiǎn)單且更加通用而愈加流行.FR 格式最先于2007 年由Huynh[11-12]提出并證明了,至少對(duì)于一維線性方程,FR 格式通過不同修正多項(xiàng)式的選擇,不僅囊括了DG 和SD 這兩種格式,并且能夠產(chǎn)生新的高階格式.2009 年,Wang和Gao 等[13]提出“l(fā)ifting collocation penalty,LCP”格式并將其推廣到歐拉方程和Navier-Stokes 方程.由于FR 和LCP 這兩種格式具有高度的相似性和內(nèi)在聯(lián)系,因此后來也被Huynh 和王志堅(jiān)[14]聯(lián)合命名為FR/CPR (correction procedure via reconstruction) 格式.2011 年,Vincent 等[15]提出了一種針對(duì)一維線性對(duì)流方程的能量穩(wěn)定的FR 格式,即VCJH 格式,后續(xù)又在此格式的基礎(chǔ)上開發(fā)了三維并行CFD 求解器Hi-FiLES[16]和PyFR[17].如今,FR 已不僅是一種高階格式,它還為幾種現(xiàn)有的高階格式,例如SD 和DG,提供了一個(gè)統(tǒng)一的框架.相比DG 和SD 的原有形式,FR 在概念上更加簡(jiǎn)單且計(jì)算效率更高[18].關(guān)于FR和DG 間的聯(lián)系,Allaneau 和Jameson[19]和De Grazia等[20]已經(jīng)進(jìn)行了深入研究.

盡管FR 格式已經(jīng)取得了較大的進(jìn)步,但將其應(yīng)用于復(fù)雜外形流動(dòng)的數(shù)值模擬時(shí),仍然有一些難點(diǎn)問題需要解決,例如計(jì)算效率的問題、曲面邊界處理方法、間斷偵測(cè)方法與限制器的構(gòu)造等.其中,計(jì)算效率的問題在很大程度上制約了FR 格式在實(shí)際工程中的應(yīng)用.在當(dāng)前的FR 格式應(yīng)用中,時(shí)間離散多采用存儲(chǔ)需求小且易實(shí)現(xiàn)的顯式Runge-Kutta 類方法[21].但隨著格式精度的提高,FR 對(duì)應(yīng)的CFL 穩(wěn)定性條件越來越嚴(yán)格,時(shí)間推進(jìn)步長(zhǎng)會(huì)受到明顯的限制.對(duì)一些復(fù)雜流動(dòng)問題(流場(chǎng)中存在邊界層、剪切層等小尺度流動(dòng)結(jié)構(gòu))的計(jì)算更是如此.因此,發(fā)展FR 格式的加速收斂方法,具有十分重要的工程應(yīng)用價(jià)值.目前常用的加速收斂方法在算法層面主要有隱式時(shí)間推進(jìn)方法、當(dāng)?shù)貢r(shí)間步長(zhǎng)、殘差光順、多重網(wǎng)格以及物理層面的并行計(jì)算技術(shù)等.在這些方法中,隱式時(shí)間推進(jìn)方法的加速效果是十分明顯的[22].相比顯式時(shí)間方法,隱式方法受CFL 穩(wěn)定性條件的約束小,可以選取更大的時(shí)間推進(jìn)步長(zhǎng),進(jìn)而能夠有效地提高CFD 求解的效率.近年來已有學(xué)者提出了多種基于FR 格式的隱式時(shí)間方法,例如LU-SGS 方法[23]、預(yù)處理GMERS 方法[24]、無矩陣(Jacobianfree) Newton Krylov 方法[25]以及多色高斯-賽德爾(multi-colored Gauss Seidel,MCGS) 方法[26]等.這些文章中以數(shù)值實(shí)驗(yàn)證明,使用隱式方法的FR 格式相比顯式具有非常明顯的計(jì)算效率提升.以上這些隱式時(shí)間方法大都是在CPU 架構(gòu)上實(shí)現(xiàn),由于基于局部單元近似的FR 格式的絕大多數(shù)操作都是在局部單元中進(jìn)行,計(jì)算量較為集中,因此更適合GPU 計(jì)算.當(dāng)前,基于CPU 并行策略的有限體積求解器一般只能使用并行集群性能峰值的5%～10%,而使用GPU 計(jì)算的FR 格式計(jì)算速度最高可以達(dá)到性能峰值的58%左右[27].因此,基于GPU 并行的隱式時(shí)間推進(jìn)方法有望進(jìn)一步提高FR 格式的計(jì)算效率.

近年來,GPU 計(jì)算在CFD 領(lǐng)域中的應(yīng)用已經(jīng)越來越普遍.已有的研究結(jié)果表明[27-29],對(duì)于相同規(guī)模的問題,GPU 并行能夠較大程度地提升高階CFD 程序的計(jì)算效率.例如,Romero 等[30]利用顯式時(shí)間方法的VCJH 格式計(jì)算非定常流動(dòng)問題,發(fā)現(xiàn)GPU 并行相比基于MPI 并行的CPU 節(jié)點(diǎn)具有更高的計(jì)算效率.Vermeire 等[31]使用PyFR 求解器在多個(gè)GPU 卡上以三維Taylor-Green 渦和SD7003 翼型為例進(jìn)行隱式大渦模擬計(jì)算,發(fā)現(xiàn)與使用CPU 節(jié)點(diǎn)并行的商業(yè)CFD 軟件Star-CCM+相比,在計(jì)算成本相當(dāng)?shù)那疤嵯?5 階FR 計(jì)算得到的模擬結(jié)果更精細(xì).最近,Jourdan 和Wang 等[32]在美國(guó)“Summit”超級(jí)計(jì)算機(jī)上使用FR/CPR 高階求解器,分別利用了800 個(gè)節(jié)點(diǎn)的CPU 和4800 個(gè)GPU 卡進(jìn)行了大渦模擬計(jì)算,發(fā)現(xiàn)使用顯式時(shí)間推進(jìn)方法時(shí),在六面體和四面體這兩種網(wǎng)格單元類型上,GPU 相比CPU 并行分別有一個(gè)和兩個(gè)量級(jí)左右的加速效果.

在上述的大規(guī)模CFD 數(shù)值模擬中,基于FR 格式的GPU 計(jì)算使用的都為顯式時(shí)間方法.為了探索如何進(jìn)一步提高FR 格式的求解效率,本文以二維歐拉方程為例,提出一種基于高階FR 格式的單GPU并行隱式時(shí)間求解方法.通過塊Jacobi 迭代的方式,改變了經(jīng)過隱式離散后全局線性方程組左端矩陣的特征,克服了影響求解并行性的相鄰單元依賴.最終,將解全局線性方程組轉(zhuǎn)化解一系列局部單元線性方程組,進(jìn)而又通過LU 分解法并行求解這些局部線性方程組.數(shù)值實(shí)驗(yàn)表明,相比顯式方法而言,該隱式方法在計(jì)算效率方面至少有一個(gè)量級(jí)的提升.并且,該隱式方法不局限于FR 格式,對(duì)DG 和SD 等基于局部單元近似的高階格式也具有參考意義.

1 控制方程

無黏流動(dòng)的控制方程為歐拉方程,其二維形式如下

其中,U=[ρ,ρu,ρv,ρE]為守恒變量,F和G分別為x和y方向的無黏通量

式中,ρ 是密度,u和v分別表示x和y方向的速度,p為壓強(qiáng),E代表單位質(zhì)量流體的總能量。本文中假設(shè)流體為理想氣體,則有

其中,γ=Cp/Cv為比熱容比,本文中其值為1.4.

2 高階通量重構(gòu)格式

本節(jié)將以二維標(biāo)量守恒律為例來簡(jiǎn)要地介紹FR格式的基本思想.在區(qū)域Ω ∈R2上的一般形式的二維標(biāo)量守恒律形式為

其中,u=u(x,y,t)為標(biāo)量守恒量,f=(f,g)為通量項(xiàng),f=f(u)和g=g(u)分別為f在x和y方向的分量.

下面將區(qū)域Ω 分割為N個(gè)互不重疊的四邊形子區(qū)域且子區(qū)域滿足如下條件

為了數(shù)學(xué)書寫的方便以及計(jì)算效率方面考慮,首先通過如圖1 所示的變換將物理域(x,y)上任意Ωn轉(zhuǎn)換到計(jì)算域(ξ,η)上的一個(gè)標(biāo)準(zhǔn)單元ΩS(為了敘述簡(jiǎn)潔,下面將省略下標(biāo)n).

圖1 將任意四邊形單元變換到一個(gè)標(biāo)準(zhǔn)單元Fig.1 Mapping of an arbitrary physical space quadrilateral element to the standard reference quadrilateral element

經(jīng)過以上變換后,子區(qū)域Ωn上的方程(4) 變換為

其中,頂標(biāo)“?”代表其為計(jì)算域(ξ,η)上的變量,計(jì)算域(ξ,η)和物理域(x,y)上的變量存在如下關(guān)系

其中,右上標(biāo)“δ”代表近似解,“D”表示是單元上的局部解,是求解點(diǎn)(ξi,ηj)處的近似解,li和lj是相對(duì)應(yīng)的拉格朗日插值多項(xiàng)式.在單元內(nèi)部是連續(xù)的,但在相鄰單元之間是間斷的.圖2 中給出了p=2 時(shí)在四邊形單元上求解點(diǎn)和后續(xù)用于通量修正的邊界面通量點(diǎn)的分布.

圖2 兩個(gè)四邊形單元上求解點(diǎn)(圓圈)和通量點(diǎn)(方塊)的分布(p=2)Fig.2 Arrangement of solution points(circles)and flux points(squares)on a pair of quadrilateral elements(p=2)

在計(jì)算散度項(xiàng)之前,首先進(jìn)行通量項(xiàng)的構(gòu)造和計(jì)算.參照方程(9) 中近似解的構(gòu)造,同樣可以獲取局部通量.但由于僅在單元內(nèi)部連續(xù),而在單元之間并不連續(xù).為了滿足格式守恒性的要求,FR 通過增加修正項(xiàng)的方式構(gòu)造在單元界面處連續(xù)的通量.因此,連續(xù)的全局通量由單元內(nèi)的局部通量和修正通量構(gòu)成上式中,局部通量通過p階的插值多項(xiàng)式進(jìn)行近似

綜上,經(jīng)過FR 格式空間離散之后的半離散方程為

上述常微分方程組可進(jìn)一步通過顯式或隱式的時(shí)間推進(jìn)方法進(jìn)行求解.

3 隱式時(shí)間推進(jìn)方法

二維歐拉方程經(jīng)過空間離散之后得到的半離散方程可進(jìn)一步寫作如下形式

在上式中為M×N×4 維度的全局解向量,其中,M和N分別代表網(wǎng)格單元總數(shù)以及每個(gè)單元內(nèi)求解點(diǎn)的個(gè)數(shù),UM(t)為第M個(gè)單元上的局部解,是第M個(gè)單元第N個(gè)求解點(diǎn)上包含守恒變量的解向量.

代表右端殘差項(xiàng),RM(U)是第M個(gè)單元上的局部殘差,為第M個(gè)單元第N點(diǎn)上殘差.

由于本文計(jì)算的是定常流動(dòng),無需考慮時(shí)間精度,因此時(shí)間項(xiàng)的離散使用一階向后隱式差分格式.方程(15)經(jīng)過該方法離散后可得

式中,Un=U(tn)和Un+1=U(tn+1)分別為tn和tn+1時(shí)刻的解.ΔTn為一個(gè)包括時(shí)間步長(zhǎng)的對(duì)角矩陣,其與Un具有相同的維數(shù),其具體形式為

由于方程(16)中的R(Un+1)是非線性的,通過泰勒展開對(duì)其線性化,則有

上式中,?R/?U|Un為全局雅可比矩陣.將方程(18)代入方程(16)并整理得

上式中,tn+1時(shí)的解Un+1需要先求解線性方程組才能獲取.

3.1 塊雅可比迭代

線性方程組的求解可分為直接法和迭代法.當(dāng)網(wǎng)格規(guī)模較大時(shí),使用直接法求解效率低下且內(nèi)存占用較大.因此本文通過迭代法求解方程(19),所用的迭代方法為雅可比迭代.

首先對(duì)方程(19) 左端的矩陣進(jìn)行分解.令A(yù)=[(ΔTn)-1-?R/?U|Un].則A=D-N可看作是由分塊對(duì)角矩陣D和剩余非對(duì)角部分N組成.其中,分塊對(duì)角矩陣D為

與此同時(shí),N的矩陣形式是

式中的N為一稀疏矩陣.經(jīng)過上述分解之后,對(duì)全局線性方程組(19)進(jìn)行迭代求解,可得

在上式中,k為雅可比迭代指標(biāo).當(dāng)k→∞時(shí),雅可比迭代解Uk+1趨近于下一個(gè)時(shí)刻的解Un+1.對(duì)于定常問題的計(jì)算,k只需要設(shè)置為有限的值,當(dāng)n→∞,Uk+1和Un最終都將趨向于所求的定常解.

為了避免非對(duì)角部分N的計(jì)算.對(duì)方程(23)兩側(cè)同時(shí)減去D(Uk-Un),則有

上式經(jīng)進(jìn)一步整理后可得

和方程(18)類似,將R(Uk)在R(Un)處進(jìn)行泰勒展開

上式可進(jìn)一步寫成

將方程(27)代入到方程(25)中消去全局矩陣相關(guān)的一項(xiàng)?R/?U|Un(Uk-Un),并最終整理可得

式中,D=diag{D1,D2,···,DM}.

方程(28)和方程(19)相比,它的左端矩陣D是一個(gè)分塊對(duì)角矩陣,這使得求解全局線性方程組(19)的問題解耦為求解一系列局部單元上的小型線性方程組.

值得注意的是,全局矩陣A的分解方式并不唯一,如將其分解為分塊上三角矩陣U和分塊下三角矩陣L,即為Gauss-Seidel 迭代法.此時(shí)需要求解的方程的左端矩陣為分塊下三角矩陣L.由于此時(shí)Uk+1存在數(shù)據(jù)依賴問題,無法同時(shí)更新,為了并行地求解該方程,需要使用圖著色技術(shù),所得方法即為更復(fù)雜的MCGS 法[26].

3.2 局部線性方程組的求解

方程(28)可進(jìn)一步寫成一系列局部單元線性方程組.對(duì)于第i(1 ≤i≤M)個(gè)單元,則有上式中,i為單元指標(biāo),n是時(shí)間迭代指標(biāo),k為雅可比迭代指標(biāo).為第i個(gè)單元對(duì)應(yīng)的單元雅可比矩陣,在下一小節(jié)中將討論其具體計(jì)算方法.以上獲取的M個(gè)局部線性方程組(29),其左端矩陣規(guī)模較小,在程序中通過調(diào)用cuBLAS 函數(shù)庫(kù)使用LU 分解法實(shí)現(xiàn)同步求解.

3.3 單元雅可比矩陣的計(jì)算

本文中所用是解析法與自動(dòng)微分相結(jié)合的方法.由式(13)可知,第i個(gè)單元上,FR 格式離散后的殘差Ri是不連續(xù)通量的散度和修正通量的散度之和.則單元雅可比矩陣為

4 數(shù)值實(shí)驗(yàn)

本節(jié)中以二維Bump 流動(dòng)和NACA0012 翼型無黏繞流兩個(gè)算例來展示隱式方法的精度和效率.為了進(jìn)行計(jì)算效率的對(duì)比研究,所有算例均同時(shí)使用顯式和隱式兩類方法進(jìn)行.顯式方法為具有TVD 性質(zhì)的三級(jí)Runge-Kutta(TVD RK3)方法[21],同時(shí)使用局部時(shí)間步長(zhǎng)和多重網(wǎng)格(pmulti-grid)方法加速其收斂.隱式方法為本文中所提出的塊雅可比迭代的方法.隱式方法的時(shí)間步長(zhǎng)ΔTn由CFL 數(shù)決定,CFL數(shù)初始值為2,并隨迭代步數(shù)指數(shù)增長(zhǎng),最大值為104.

所有算例均使用高階并行通量重構(gòu)軟件HAFR(heterogeneous architecture flux reconstruction) 在型號(hào)為英偉達(dá)GeForce RTX 2080ti 的單個(gè)GPU 卡上計(jì)算.為了保證FR 方法在邊界的精度不至于下降,算例中用到的網(wǎng)格是由開源網(wǎng)格生成軟件Gmsh 生成的二階曲邊網(wǎng)格.計(jì)算的收斂性通過密度殘差的L2范數(shù)判斷,當(dāng)其下降到10-14以下或平穩(wěn)而不再下降時(shí)便認(rèn)為計(jì)算已穩(wěn)態(tài)收斂并終止.

4.1 無黏亞聲速Bump 算例

本算例是 international workshop on high-order CFD methods[39]給定的一個(gè)二維無黏驗(yàn)證算例.該問題的物理域?yàn)閇-1.5,1.5]×[0,0.8],底邊的幾何描述為y=0.062 5e-25x2.左側(cè)為亞聲速入口邊界條件(給定入口總壓和總溫),右側(cè)出口為出口邊界條件(給定靜壓),上下邊界均為滑移壁面邊界條件.來流的馬赫數(shù)Ma=0.5 且攻角為零度.盡管該流動(dòng)的精確解未知,但由于流動(dòng)為亞聲速無黏流動(dòng)且流場(chǎng)光滑,因此可通過熵來衡量解的精度,流場(chǎng)中熵的計(jì)算公式如下

式中的積分在HAFR 程序中通過10 階高斯積分法進(jìn)行計(jì)算.

計(jì)算中所使用的網(wǎng)格由粗到細(xì)(即h-refinement)為：R1(12×4),R2(24×8),R3(48×16)和R4(96×32),圖3 是R2 網(wǎng)格的示例.對(duì)于同一網(wǎng)格,同時(shí)還將使用p=2,3,4 階的插值多項(xiàng)式(即p-refinement)進(jìn)行計(jì)算.

圖3 本算例中所用的24×8 的四邊形網(wǎng)格Fig.3 Quadrangular mesh(24×8)used in the case

首先從數(shù)值精度的角度分析計(jì)算結(jié)果.表1 為HAFR 程序使用塊雅可比迭代隱式方法在四組網(wǎng)格上,插值多項(xiàng)式為p=2,3,4 階時(shí)的熵以及數(shù)值精度.從表中可知,對(duì)任意一階插值多項(xiàng)式階數(shù)而言,網(wǎng)格加密的過程中熵也隨之降低.同時(shí),給定某組網(wǎng)格,提高插值多項(xiàng)式的階數(shù)也同樣會(huì)導(dǎo)致熵的降低.正如圖4 所示,在R2 網(wǎng)格上使用p=3 和p=4 階插值多項(xiàng)式得到的下壁面馬赫數(shù)等值線相比p=2 的更為光滑.此外還觀察到,對(duì)于p=2,3,4 階插值多項(xiàng)式,在4 組網(wǎng)格上獲取的數(shù)值精度,都近似達(dá)到p+1階這一理論值,這也驗(yàn)證了HAFR 的計(jì)算精度.

圖4 在24×8 的網(wǎng)格上使用p=2,3,4 階插值多項(xiàng)式計(jì)算所得的馬赫數(shù)云圖Fig.4 Mach number contours calculated on the 24×8 mesh with p=2,3,4 polynomial orders

表1 在四組網(wǎng)格上,使用p=2,3,4 階多項(xiàng)式時(shí)的熵和數(shù)值精度Table 1 Entropies and numerical orders on four different grids with p=2,3,4 polynomial orders

下面從計(jì)算效率的角度進(jìn)行分析.圖5 是HAFR程序使用顯式和隱式兩種方法在4 組網(wǎng)格上、插值多項(xiàng)式為p=2,3,4 階時(shí)密度殘差隨迭代次數(shù)的收斂情況.

在圖5 中,左右兩列子圖分別為使用塊雅可比迭代隱式和TVD RK3 顯式時(shí)間方法的計(jì)算結(jié)果.從迭代次數(shù)來看,對(duì)同一階插值多項(xiàng)式,顯式方法和隱式方法收斂所需的迭代次數(shù)都隨著網(wǎng)格的加密而增大.對(duì)于同一網(wǎng)格,使用更高階的插值多項(xiàng)式都會(huì)導(dǎo)致迭代次數(shù)的增多.但同時(shí)發(fā)現(xiàn),隱式方法的迭代次數(shù)隨自由度(degree of freedom,DoF)增幅較小,并且對(duì)于同一網(wǎng)格或同一階插值多項(xiàng)式,隱式方法所需的迭代次數(shù)要遠(yuǎn)遠(yuǎn)小于顯式方法.例如,在R1 網(wǎng)格上,隱式方法收斂的的迭代次數(shù)均在80 步以內(nèi),而顯式方法最快也需要1000 步左右才能收斂.

圖5 在四組網(wǎng)格上,插值多項(xiàng)式為p=2,3,4 階時(shí),使用隱式和顯式方法所得到的密度殘差||Res(ρ)||2 隨迭代次數(shù)的變化Fig.5 Density residual||Res(ρ)||2 versus the iteration numbers for the explicit and implicit time-marching schemes on four different grids with p=2,3,4 polynomial orders

圖5 在四組網(wǎng)格上,插值多項(xiàng)式為p=2,3,4 階時(shí),使用隱式和顯式方法所得到的密度殘差||Res(ρ)||2 隨迭代次數(shù)的變化(續(xù))Fig.5 Density residual||Res(ρ)||2 versus the iteration numbers for the explicit and implicit time-marching schemes on four different grids with p=2,3,4 polynomial orders(continued)

在密度殘差收斂方面,使用隱式方法的計(jì)算的密度殘差均能收斂到10-14這一量級(jí).但顯式方法的收斂性隨著自由度的增大變差,其值下降到某個(gè)量級(jí)后將無法下降,例如使用p=2 階插值多項(xiàng)式時(shí),在R1 網(wǎng)格上的密度殘差能夠收斂10-14量級(jí),但在R4 網(wǎng)格上僅能下降到10-13這一量級(jí).

在GPU 計(jì)算時(shí)間方面,表2 為HAFR 程序采用塊雅可比迭代隱式方法在4 組網(wǎng)格上,使用p=2,3,4 階插值多項(xiàng)式計(jì)算收斂所用的迭代次數(shù)和GPU 墻鐘時(shí)間.從表2 可知,隱式方法計(jì)算收斂所用時(shí)間隨著自由度的增大而增大.

表2 在四組網(wǎng)格上,插值多項(xiàng)式為p=2,3,4 階時(shí),使用隱式時(shí)間方法時(shí)的迭代次數(shù)和墻鐘時(shí)間Table 2 Iteration numbers and wall-clock times for the implicit time-marching scheme on four different grids with p=2,3,4 polynomial orders

圖6 是HAFR 程序分別使用塊雅可比隱式和TVD RK3 顯式時(shí)間方法,在不同網(wǎng)格上,使用p=2,3,4 階插值多項(xiàng)式時(shí)所用的墻鐘時(shí)間對(duì)比.從圖6中可看到,對(duì)每一組計(jì)算,隱式方法收斂所用的時(shí)間明顯小于顯式方法.在所有算例中,隱式方法相比顯式至少有一個(gè)量級(jí)的加速效果.

圖6 在四組網(wǎng)格上,插值多項(xiàng)式為p=2,3,4 階時(shí),使用隱式和顯式方法所用的計(jì)算時(shí)間Fig.6 Wall-clock times for the explicit and implicit time-marching schemes on four different grids with p=2,3,4 polynomial orders

4.2 NACA0012 翼型無黏繞流

在上一個(gè)Bump 算例中所使用的網(wǎng)格較為均勻,為了進(jìn)一步考察在非均勻網(wǎng)格條件下塊雅可比迭代隱式方法的計(jì)算效率,下面對(duì)NACA0012 翼型無黏繞流進(jìn)行數(shù)值模擬.在此算例中,自由來流的馬赫數(shù)為0.5 且攻角為零度.NACA0012 翼型表面為滑移壁面條件,遠(yuǎn)場(chǎng)為特征邊界條件,遠(yuǎn)場(chǎng)與翼型的距離大約為100 個(gè)弦長(zhǎng).

計(jì)算中用到的網(wǎng)格為Gmsh 生成的O 型網(wǎng)格,網(wǎng)格分為3 組：R1(32×32),R2(44×64)和R3(128×128).圖7 為R1(32×32) 網(wǎng)格的示意圖.對(duì)于同一網(wǎng)格,還將分別使用p=2,3,4 階的插值多項(xiàng)式(即prefinement)進(jìn)行計(jì)算.

圖7 本算例中所用的32×32 的O 型四邊形網(wǎng)格Fig.7 O-type quadrangular mesh(32×32)used in this case

下面對(duì) NACA0012 數(shù)值模擬結(jié)果進(jìn)行分析.表3 為HAFR 程序采用塊雅可比迭代隱式方法在三組網(wǎng)格上,使用p=2,3,4 階插值多項(xiàng)式計(jì)算收斂所用的迭代次數(shù)和GPU 墻鐘時(shí)間.和上個(gè)算例中類似,可觀察到無論是網(wǎng)格加密還是提高插值多項(xiàng)式階數(shù),這兩種增大自由度的方式都會(huì)導(dǎo)致迭代次數(shù)和GPU計(jì)算時(shí)間的增多.

表3 在三組網(wǎng)格上,插值多項(xiàng)式為p=2,3,4 階時(shí),使用隱式時(shí)間方法時(shí)的迭代次數(shù)和墻鐘時(shí)間Table 3 Iteration numbers and wall-clock times for the implicit time-marching schemes on three different grids with p=2,3,4 polynomial orders

圖8 是HAFR 程序分別使用隱式和顯式時(shí)間方法,在不同網(wǎng)格上,使用插值多項(xiàng)式為p=2,3,4 時(shí)所用的墻鐘時(shí)間對(duì)比.從圖8 中可知,對(duì)于每一組計(jì)算,使用塊雅可比迭代隱式方法收斂所用的時(shí)間相比TVD RK3 顯式方法至少有一個(gè)量級(jí)的加速效果.

圖8 在三組網(wǎng)格上,插值多項(xiàng)式為p=2,3,4 階時(shí),使用隱式和顯式方法所用的墻鐘時(shí)間Fig.8 Wall-clock times for the explicit and implicit time-marching schemes on three different grids with p=2,3,4 polynomial orders

圖9 為HAFR 程序使用顯式和隱式兩種方法在三組網(wǎng)格上、插值多項(xiàng)式為p=2,3,4 階時(shí)密度殘差隨迭代次數(shù)的收斂情況.在圖9 中,左右兩列子圖分別為使用塊雅可比迭代隱式和TVD RK3 顯式時(shí)間方法的結(jié)果.從迭代次數(shù)來看,對(duì)同一階插值多項(xiàng)式,顯式方法和隱式方法收斂所需的迭代次數(shù)都隨著網(wǎng)格的加密而增大.對(duì)于同一網(wǎng)格,使用更高階的插值多項(xiàng)式會(huì)導(dǎo)致迭代次數(shù)的增多.但同時(shí),相比顯式方法,隱式方法的迭代次數(shù)隨自由度增幅較小,并且對(duì)于同一網(wǎng)格或同一階插值多項(xiàng)式,隱式方法所需的迭代次數(shù)要遠(yuǎn)遠(yuǎn)小于顯式方法.

圖9 在三組網(wǎng)格上,插值多項(xiàng)式為p=2,3,4 階時(shí),使用隱式和顯式方法所得到的密度殘差||Res(ρ)||2 隨迭代次數(shù)的變化Fig.9 Density residual||Res(ρ)||2 versus the iteration numbers for the explicit and implicit time-marching schemes on three different grids with p=2,3,4 polynomial orders

在密度殘差收斂方面,在3 組網(wǎng)格計(jì)算中,使用塊雅可比迭代隱式方法計(jì)算的密度殘差均能收斂到10-14這一量級(jí).但TVD RK3 顯式時(shí)間方法的收斂性隨著自由度的增大變差,其值下降到某個(gè)量級(jí)后將無法下降.因此從殘差收斂角度,基于塊雅可比迭代的隱式方法要優(yōu)于TVD RK3 顯式方法.

綜上,隨著自由度的增大,顯式和隱式方法收斂所需的迭代次數(shù)都會(huì)增多.但對(duì)于同一網(wǎng)格或同一階插值多項(xiàng)式,隱式方法所需的迭代次數(shù)和計(jì)算時(shí)間要明顯小于顯式方法,且隱式方法的殘差下降幅度要比顯式方法更大.

5 結(jié)論

為了提高FR 方法的求解效率,本文提出一種基于雅可比塊迭代的高階FR 格式求解定常二維歐拉方程的單GPU 隱式時(shí)間推進(jìn)方法.通過雅可比塊迭代的方式,克服了影響求解并行性的相鄰單元依賴,使得只需存儲(chǔ)和計(jì)算對(duì)角塊矩陣.最終將求解全局線性方程組轉(zhuǎn)化為求解一系列局部單元線性方程組.這些局部的線性方程組在GPU 上又可利用LU 分解法并行求解.

文中通過二維無黏Bump 流動(dòng)和NACA0012 無黏繞流兩個(gè)算例表明:隨著自由度的增大,塊雅可比迭代隱式方法收斂所需的迭代次數(shù)和計(jì)算時(shí)間都會(huì)增多.但對(duì)于同一網(wǎng)格或同一階插值多項(xiàng)式,該隱式方法所需的迭代次數(shù)要明顯小于TVD RK3 顯式方法,且此隱式方法的殘差下降幅度要比TVD RK3 顯式方法更大.在計(jì)算時(shí)間方面,對(duì)于文中的每一組計(jì)算,雅可比塊迭代隱式方法收斂所用的GPU 計(jì)算時(shí)間明顯小于TVD RK3 顯式方法,且計(jì)算效率至少有一個(gè)量級(jí)的提高.

附錄

二維歐拉方程的雅克比矩陣為