基于推廣的Fu-Shu壞單元指示子的h自適應(yīng)RKDG方法

韓文秀, 王海云, 朱洪強(qiáng)

(南京郵電大學(xué)理學(xué)院, 南京 210023)

求解雙曲守恒律方程的Runge-Kutta間斷Galerkin方法(the Runge-Kutta discontinuous Galerkin method, RKDG)[1]具有精度高、易于處理復(fù)雜計(jì)算區(qū)域和邊界條件、求解高效、并行效率高、易于進(jìn)行hp自適應(yīng)等優(yōu)點(diǎn), 它已成為當(dāng)今科學(xué)計(jì)算領(lǐng)域高分辨方法的研究前沿和熱點(diǎn),也是計(jì)算流體力學(xué)的主流方法之一．非線性雙曲守恒律方程的解經(jīng)常含有間斷,如果使用固定網(wǎng)格來求解,為了提高間斷的求解效果,只能在全部計(jì)算區(qū)域上使用密網(wǎng)格,從而導(dǎo)致計(jì)算量大幅增加．自適應(yīng)方法是解決這個(gè)問題的一種有效策略,它能根據(jù)解的情況按一定要求以自適應(yīng)的方式重新生成網(wǎng)格,從而節(jié)約計(jì)算成本．h自適應(yīng)RKDG方法目前已有不少研究工作, 如Remacle[2]、Hartmann[3]和Dedner[4]等人分別基于誤差估計(jì)生成自適應(yīng)網(wǎng)格．另一種h自適應(yīng)策略基于所謂的壞單元指示子, 包括Devine[5]、Zhu[6-7]等人的研究工作, 這種方法因?yàn)椴恍枰褂幂^難得到的誤差估計(jì)而相對(duì)簡(jiǎn)單, 同時(shí)也能實(shí)現(xiàn)良好的自適應(yīng)效果,故近些年仍不斷出現(xiàn)新的研究成果[8-11]．

RKDG方法需要限制器[1]來控制間斷附近的數(shù)值偽振蕩, 以保持格式穩(wěn)定．壞單元指可能含有間斷,需要限制器發(fā)揮作用的單元,而壞單元指示子通常來自各種限制器或激波探測(cè)技術(shù)．Zhu等[6]利用壞單元指示子將每個(gè)單元標(biāo)記為壞單元或好單元,采用加密壞單元、合并好單元的策略,實(shí)現(xiàn)了在間斷區(qū)域使用密網(wǎng)格、連續(xù)區(qū)域使用粗網(wǎng)格的自適應(yīng)效果,有效提高了數(shù)值解的質(zhì)量．區(qū)別于普通無懸點(diǎn)的網(wǎng)格,基于單元加密和合并的二維h自適應(yīng)網(wǎng)格將不可避免地出現(xiàn)懸點(diǎn),而且不同單元可能會(huì)有不同數(shù)量和大小的直接鄰居單元, 這種幾何復(fù)雜性導(dǎo)致絕大多數(shù)壞單元指示子無法直接用于h自適應(yīng)網(wǎng)格．常用的KXRCF壞單元指示子[12]是少數(shù)能直接應(yīng)用的指示子, 但文獻(xiàn)[6]中的數(shù)值結(jié)果表明,隨著逼近空間次數(shù)的增加, KXRCF指示子傾向于標(biāo)記越來越多的壞單元,導(dǎo)致連續(xù)區(qū)域也會(huì)使用密網(wǎng)格,顯著影響自適應(yīng)的效果．Fu等[13]提出1個(gè)新的壞單元指示子,該指示子做法簡(jiǎn)單,模板緊湊,不含有依賴于問題的敏感參數(shù)且數(shù)值表現(xiàn)良好; Zhu等[14]將它推廣到二維h自適應(yīng)網(wǎng)格,同時(shí)保持了原來的優(yōu)點(diǎn)．本文擬將這個(gè)推廣的Fu-Shu壞單元指示子用于h自適應(yīng)RKDG方法, 從而進(jìn)一步提高數(shù)值解的質(zhì)量,節(jié)約計(jì)算成本．

1 二維RKDG方法及h自適應(yīng)方法

RKDG方法對(duì)時(shí)間使用非線性穩(wěn)定的顯式RK方法離散, 對(duì)空間使用DG離散, 采用分片多項(xiàng)式空間作為近似解和測(cè)試函數(shù)空間．考慮二維標(biāo)量雙曲守恒律方程

(1)

該方法中限制器的壞單元指示子采用推廣的Fu-Shu壞單元指示子, 限制器的數(shù)值解重構(gòu)采用WENO重構(gòu)方法[6]．

2 推廣的Fu-Shu壞單元指示子

圖1 目標(biāo)單元K0有7個(gè)直接鄰居單元的示例圖Fig.1 An example of the target cell K0 with 7 immediate neighbors

3 數(shù)值結(jié)果

為了研究推廣的Fu-Shu壞單元指示子用于生成h自適應(yīng)網(wǎng)格, 本文對(duì)二維可壓Euler方程組

(2)

進(jìn)行廣泛的數(shù)值試驗(yàn), 其中ρ表示密度,u和v分別為x方向和y方向的速度,E是總能量,p是壓力, 滿足p=(γ-1)[E-0.5ρ(u2+v2)], 式中γ=1.4．文獻(xiàn)[6]中使用KXRCF壞單元指示子生成h自適應(yīng)網(wǎng)格, 在k=1時(shí)得到了質(zhì)量較高的網(wǎng)格, 但在k=2時(shí)因?yàn)閴膯卧袛噙^多導(dǎo)致加密太多,自適應(yīng)效果仍有改進(jìn)的空間．在本文的數(shù)值試驗(yàn)中,k=1時(shí)推廣的Fu-Shu壞單元指示子得到了與KXRCF指示子效果相當(dāng)?shù)淖赃m應(yīng)網(wǎng)格, 故下面僅討論k=2時(shí)的數(shù)值解和自適應(yīng)網(wǎng)格圖,以展示推廣的Fu-Shu指示子在高階時(shí)的優(yōu)勢(shì)．在所有的數(shù)值試驗(yàn)中均取Lmax=4．

例1黎曼問題．求解帶有以下初值條件的方程(2):

圖2給出了初始單元尺寸Δx0=Δy0=1/80,t=0.25時(shí)的密度等高線圖和網(wǎng)格圖．從圖2中可以看出, 推廣的Fu-Shu壞單元指示子精準(zhǔn)捕捉到了間斷, 并在間斷區(qū)域使用了最密網(wǎng)格,而其他區(qū)域使用了粗網(wǎng)格, 達(dá)到了預(yù)期的自適應(yīng)效果．對(duì)比文獻(xiàn)[6]圖6中k=2時(shí)的密度等高線圖和網(wǎng)格圖, 推廣的Fu-Shu指示子得到的結(jié)果在間斷處的加密區(qū)域更窄, 也沒有連續(xù)區(qū)域被誤加密的問題．

圖2 黎曼問題在t=0.25時(shí)的密度等高線圖(a)和網(wǎng)格圖(b)Fig.2 Density contours (a) and mesh (b) at t=0.25 for the Riemann problem

例2雙馬赫反射問題．取計(jì)算區(qū)域[0,4]×[0,1], 底邊1/6≤x≤4區(qū)域?yàn)楣瘫? 初始時(shí)與底邊交于(1/6,0)并與底邊成60°夾角的激波以10馬赫的速度撞擊固壁．波前未受干擾的空氣密度為1.4, 壓力為1．底邊0≤x≤1/6區(qū)域使用精確波后邊界條件, 其余區(qū)域使用反射邊界條件．頂部的邊界條件使用描述10馬赫激波的精確邊界條件, 左右邊界分別使用入流和出流邊界條件．圖3給出了Δx0=Δy0=1/60,t=0.2時(shí)的密度等高線圖和網(wǎng)格圖．圖3顯示,最終網(wǎng)格和間斷線相吻合, 自適應(yīng)方法在光滑區(qū)域生成了初始粗網(wǎng)格,而在間斷附近生成了充分加密的網(wǎng)格, 自適應(yīng)效果良好．對(duì)比文獻(xiàn)[6]圖9中k=2時(shí)的密度等高線圖和網(wǎng)格圖, 推廣的Fu-Shu壞單元指示子只在必要的間斷附近加密網(wǎng)格, 被加密區(qū)域小很多, 相應(yīng)的計(jì)算成本也更低．

圖3 雙馬赫反射問題在t=0.2時(shí)的密度等高線圖(a)和網(wǎng)格圖(b)Fig.3 Density contours (a) and mesh (b) at t=0.2 for the double Mach reflection problem

例3前臺(tái)階問題．問題始于長(zhǎng)度為3, 高度為1的風(fēng)洞中, 風(fēng)洞底部的臺(tái)階高0.2,位于距離風(fēng)洞左側(cè)0.6處并一直延伸到風(fēng)洞右端．初始時(shí), 速度為3馬赫的均勻空氣向右流動(dòng),空氣密度為1.4,壓力為1．風(fēng)洞壁和臺(tái)階處為反射邊界,左右邊界分別為入流和出流邊界．臺(tái)階角為奇點(diǎn),采用與文獻(xiàn)[13]同樣的方法處理．圖4給出了Δx0=Δy0=1/40,t=4.0時(shí)的密度等高線圖和網(wǎng)格圖．從圖4中的結(jié)果以及對(duì)比文獻(xiàn)[6]圖7中k=2時(shí)的結(jié)果可知, 推廣的Fu-Shu壞單元指示子捕捉間斷更加精準(zhǔn), 并生成了質(zhì)量更高的自適應(yīng)網(wǎng)格．

圖4 前臺(tái)階問題在t=4.0時(shí)的密度等高線圖(a)和網(wǎng)格圖(b)Fig.4 Density contours (a) and mesh (b) at t=4.0 for the forward-facing step problem

例4激波擴(kuò)散問題．計(jì)算域?yàn)閇0,1]×[6,11]∪[1,13]×[0,11]．激波初始時(shí)位于x=0.5, 6

圖5 激波擴(kuò)散問題在t=2.3時(shí)的密度等高線圖(a)和網(wǎng)格圖(b)Fig.5 Density contours (a) and mesh (b) at t=2.3 for the shock diffraction problem

表1 自適應(yīng)網(wǎng)格數(shù)據(jù)表