基于高階物面近似的自適應(yīng)間斷有限元法歐拉方程數(shù)值模擬

孫 強，呂宏強，伍貽兆

(南京航空航天大學，航空宇航學院，江蘇南京 210016)

基于高階物面近似的自適應(yīng)間斷有限元法歐拉方程數(shù)值模擬

孫 強，呂宏強*，伍貽兆

(南京航空航天大學，航空宇航學院，江蘇南京 210016)

將高階間斷有限元與網(wǎng)格自適應(yīng)相結(jié)合，于非結(jié)構(gòu)網(wǎng)格上求解二維Euler方程。將數(shù)值解多項式的高階項貢獻用人工粘性項系數(shù)的形式進行量化，網(wǎng)格自適應(yīng)過程中以人工粘性項系數(shù)作為網(wǎng)格自適應(yīng)的指示器。在系數(shù)達到設(shè)定的上限閥值的區(qū)域進行網(wǎng)格加密，在系數(shù)達到設(shè)定的下限閥值的區(qū)域?qū)⒌^程中加密過的網(wǎng)格稀疏以減少網(wǎng)格量。所有自適應(yīng)均在高階曲線逼近真實物面的基礎(chǔ)上進行，以保證數(shù)值結(jié)果的精度。典型數(shù)值算例結(jié)果與實驗結(jié)果進行了對比，表明采用該自適應(yīng)間斷有限元法可以保證以盡可能小的計算量得到高精度結(jié)果。

高階間斷有限元;自適應(yīng)方法;Euler方程;人工粘性;物面高階近似

0 引言

間斷Galerkin有限元法(DGM)發(fā)展于20世紀70年代［1］。由于該方法能夠構(gòu)造高精度的格式，并且易于實現(xiàn)自適應(yīng)計算和并行計算，已成為計算流體力學領(lǐng)域研究熱點之一。Bassi等在1997年用DGM成功求解了二維Euler方程和N-S方程［2-3］;Cockburn等提出Runge-Kutta間斷Galerkin方法(RKDG)［4-5］; Nguyen等將S-A湍流模型成功應(yīng)用于DGM［6］。國內(nèi)呂宏強等在稀疏的非結(jié)構(gòu)網(wǎng)格上求解了二維Euler方程和線化Euler方程［7-9］;張來平等發(fā)展了靜動態(tài)混合重構(gòu)的DG/FV混合格式［10］;于劍等對人工粘性和NS方程進行了研究［11-12］。

自適應(yīng)方法根據(jù)研究對象的不同和解的性質(zhì)自動調(diào)整網(wǎng)格分布或者數(shù)值解的逼近精度。在流體力學計算中，解通常在局部變化較大，其他大部分區(qū)域是光滑的［12］，此時使用自適應(yīng)方法可以有效降低計算量。自適應(yīng)方法通常被分為三類:對局部網(wǎng)格加密或稀疏的h型自適應(yīng);調(diào)整局部單元內(nèi)逼近多項式次數(shù)的p型自適應(yīng);調(diào)整網(wǎng)格單元位置保持單元總數(shù)不變的r型自適應(yīng)。

由于間斷有限元方法通常將未知數(shù)表達為高階形式，在相同的網(wǎng)格上計算量與有限體積法和有限差分法相比成倍增加，故而自適應(yīng)方法與間斷有限元方法相結(jié)合勢在必行。Yang等提出了一種于非結(jié)構(gòu)網(wǎng)格上模擬界面的自適應(yīng)方法［13］;Flaherty和Remacle研究小組則首次將網(wǎng)格度量場引入自適應(yīng)間斷有限元方法中［14］。徐云、吳迪、蔚喜軍［15-16］等則根據(jù)后誤差估計提出了不同的h型網(wǎng)格自適應(yīng)方法。

在高階情況下(p≥2)，間斷有限元對于物面形狀的表達十分敏感［3］，若物面邊界以直線表達，則無法得到收斂解。本文設(shè)計了一種h型自適應(yīng)高階間斷有限元法，物面邊均采用高階曲線逼近真實物面以保證數(shù)值解的精度。將數(shù)值解多項式的高階項貢獻以人工粘性［17］系數(shù)ε的形式量化，以系數(shù)ε作為網(wǎng)格自適應(yīng)的指示器。數(shù)值模擬過程中產(chǎn)生的非線性方程組采用牛頓法進行迭代，每個迭代步產(chǎn)生的大型線性系統(tǒng)采用塊高斯-賽德爾方法求解［18］。

1 間斷有限元數(shù)值離散

二維守恒形式的Euler方程為

其中，U為守恒變量，對流通量F(U)=(f(U)，g (U))，表達式分別為:

式中，ρ為密度;v1和v2為x、y方向的速度分量;E為單位總能;P為壓強;h為單位總焓。

式(1)兩邊乘檢驗函數(shù)V，在計算域運用分部積分后，弱解形式為:

將計算域劃分為網(wǎng)格{e}，其中e表示網(wǎng)格單元。在每個單元e內(nèi):

式(4)中，φi(x，y)為基函數(shù)。將式(4)代入式(3)中，得

式(6)稱為p階間斷有限元離散［7］，p為式(4)中基函數(shù)的階數(shù)。

由于間斷有限元法允許解在單元的邊界處不連續(xù)，因而數(shù)值通量的定義不是唯一的。此處的數(shù)值通量直接借用傳統(tǒng)有限體積法的處理方法，即將式(6)中的通量函數(shù)F(Uh)·n用數(shù)值通量代替，其中和分別表示該單元和其相鄰單元在該邊界處的值。本文所有計算都采用LLF格式數(shù)值通量:

在物面和遠場分別采用無穿透和無反射邊界條件。式(6)最終的離散形式可寫為:

其中，M是質(zhì)量矩陣，如果間斷有限元采用正交基函數(shù)，則是對角矩陣;R(u)是殘值;u包含Uh中所有的未知系數(shù)。對式(8)的求解采用牛頓-塊高斯賽德爾方法求解［18］。

2 人工粘性

人工粘性法［18］的思想就是在原Euler方程中加入一個耗散修正項，從而抑制激波附近的非物理振蕩。方程(1)變?yōu)槿缦滦问?

式(9)中，參數(shù)ε控制粘性項的大小。

式中，(·，·)e表示L2(Ωe)空間的標準內(nèi)積。這里預估Se的值能達到1/p4數(shù)量級。

單元的光順指示器確定后，粘性項系數(shù)通過以下函數(shù)定義:

式(11)中，se= log10Se，參數(shù) ε0= De/p，s0= log10(1/p4)，κ是人為定義常數(shù)。其中，De是計算元單元e的特征尺度。

由式(10)和式(11)可以看出，se是由U中高階項決定，表征了U的高階項的貢獻大小。解的光順程度越差se越大。εe只是作為se的分段函數(shù)，在數(shù)值表征上意義相同，它抹平光順區(qū)域的較小的se，突出非光順區(qū)域，使計算和表達更簡便。

3 網(wǎng)格自適應(yīng)

3.1 網(wǎng)格加密

在計算過程中單元的ε達到設(shè)定好的上限閥值則認為該單元需要進行加密，采用三角形單元四分法對網(wǎng)格單元進行加密，如圖1所示。

圖1 單元加密Fig.1 Element refinement

圖中實線單元為非物面單元，可以直接取三條邊的中點連線剖分加密。

Bassi指出高階間斷有限元對于物面形的表達精度十分敏感［3］，因而本文中的物面單元的物面邊均采用高階曲線逼近真實物面:在計算過程中設(shè)計物面邊

的高階曲線(本文中為6階)，即使網(wǎng)格十分稀疏，構(gòu)造好的物面曲邊依然可以十分精確的表達真實物面，如圖1中左圖虛線所示。

在物面單元需要進行剖分時，利用已構(gòu)造好的高階曲線找到逼近真實物面的物面邊的中點，然后連接另外兩條邊的中點生成四個新的網(wǎng)格單元，并且新生成的兩個物面單元的物面邊也用高階曲線擬合，如圖1中右圖虛線所示。

網(wǎng)格加密需要遵守以下原則:

1)相鄰單元的被剖分次數(shù)相差最多為1。

2)當單元的幾何尺寸細化到設(shè)定好的閥值之后就不可繼續(xù)剖分。

3.2 網(wǎng)格稀疏

在計算過程中，部分單元被加密后可能在進一步的迭代過程中解變得十分光滑，需要對其進行粗化以減少網(wǎng)格量。當加密過的單元的四個子單元ε均達到設(shè)定好的下限閥值時，則這四個子單元進行合并。如圖2所示。

圖2 單元的粗化Fig.2 Element coarsening

圖2中，實線表示非物面單元的合并過程。若粗化后的單元為物面單元，物面邊依舊采用高階曲線逼近真實物面，如圖中虛線所示。

網(wǎng)格粗化是網(wǎng)格加密的逆向操作，需要遵守以下原則:

1)初始單元不可被粗化。

2)只有被加密過的單元才可以被粗化。

3)當鄰單元已經(jīng)比判斷需要粗化的四個單元剖分次數(shù)大1時，不可進行粗化。

3.3 網(wǎng)格數(shù)據(jù)結(jié)構(gòu)

為使程序的可移植性更高，將網(wǎng)格自適應(yīng)作為模塊進行操作。雖然二叉樹結(jié)構(gòu)存儲加密和粗化網(wǎng)格直觀高效，但自適應(yīng)過程中增加與減少的點和線編號復雜且對于程序的結(jié)構(gòu)不利。因而本文所有單元、點、線的存儲結(jié)構(gòu)均為鏈表結(jié)構(gòu)，加密過程中單元的增加操作如圖3所示。

E表示需要加密的單元，生成的四個子單元中間單元依舊為編號E，其余三個則排在總單元后三個，點與邊的操作類似。在這種操作下，網(wǎng)格模塊與程序流場計算模塊的接口不受自適應(yīng)操作的影響，從而提高程序的可移植性。

同樣的，網(wǎng)格稀疏的過程單元減少過程如圖4所示。

圖3 單元增加Fig.3 Elements increase

圖4 單元減少Fig.4 Elements decrease

其中E和三個fi表示由加密操作生成的四個子單元，合并后生成的單元依舊編號為E，而fi后的單元編號依次減3。

在網(wǎng)格加密過程中，每個子單元均設(shè)置兩個標記，用于記錄所在的初始單元號和本身的加密次數(shù)。在后續(xù)的網(wǎng)格稀疏過程中，根據(jù)這兩個信息和單元相鄰關(guān)系進行單元合并。

3.4 數(shù)據(jù)傳遞

網(wǎng)格加密和粗化的過程中，使各物理量在子單元內(nèi)的積分值等于合并單元內(nèi)的積分值，從而保證物理量的守恒。假設(shè)合并單元為E，各子單元為ei，uE(x，y)為合并單元的守恒變量的函數(shù)，uei(x，y)為子單元守恒變量的函數(shù)，則需滿足以下方程:

其中m為合并單元剖分后的子單元個數(shù)。

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

4.1 圓柱繞流測試算例

為驗證物面構(gòu)造與網(wǎng)格自適應(yīng)結(jié)合的效果，本文首先對圓柱繞流(Ma=0.38、α=0)這一經(jīng)典算例進行了數(shù)值模擬。

前人［3，8］做圓柱繞流計算為獲得精度較高的結(jié)果，網(wǎng)格均人為設(shè)計為規(guī)則分布。本文為驗證物面構(gòu)造和自適應(yīng)間斷有限元的優(yōu)點，初始網(wǎng)格使用了分布較隨意、單元尺寸較極端的計算網(wǎng)格。

圖5給出了圓柱繞流網(wǎng)格，該網(wǎng)格共有80個單元，物面僅有8個點。

圖5 圓柱繞流計算網(wǎng)格Fig.5 Mesh around a circle

圖6 ε分布和自適應(yīng)后的網(wǎng)格Fig.6 ε distribution and mesh after three times adaption

圖6(a)為在原始網(wǎng)格上p=4時ε分布，可以看出主要分布在流場數(shù)值可能出現(xiàn)數(shù)值解變化比較劇烈的區(qū)域。對這一區(qū)域進行網(wǎng)格自適應(yīng)，2次自適應(yīng)后的網(wǎng)格如圖6(b)所示。可以看出網(wǎng)格主要在ε集中區(qū)域自適應(yīng)，從而僅需對這一區(qū)域的少量網(wǎng)格進行自適應(yīng)操作。自適應(yīng)前后的網(wǎng)格對比，發(fā)現(xiàn)自適應(yīng)操作確實基于真實物面曲線，自適應(yīng)后的網(wǎng)格更加逼近真實物面。

圖7給出了p=4情況下的圓柱繞流等馬赫線圖和壓強系數(shù)分布，可以看出壓強系數(shù)分布光滑并且對稱。

圖7 圓柱繞流等馬赫線圖和Cp分布Fig.7 Mach isolines and Cpdistribution

4.2 跨聲速NACA0012翼型繞流

采用上述自適應(yīng)間斷有限元法給出了經(jīng)典算例跨聲速(Ma=0.8、α=1.25)下NACA0012翼型繞流的數(shù)值結(jié)果。

圖8給出的NACA0012翼型的網(wǎng)格共有448個單元，而物面單元只有32個。

圖8 NACA0012翼型計算網(wǎng)格Fig.8 Mesh around NACA0012 airfoil

圖9給出了引入人工粘性未加入網(wǎng)格自適應(yīng)的粘性項系數(shù)ε的分布與馬赫云圖。

圖9 未加自適應(yīng)的ε分布及等馬赫線圖Fig.9 ε distribution and Mach isolines without adaption

從圖中可以看出雖然網(wǎng)格十分稀疏，但依然可以在一個網(wǎng)格尺度內(nèi)捕捉到激波，并且ε的分布集中在激波區(qū)域。

圖10給出了1次、2次自適應(yīng)后的網(wǎng)格變化。圖11給出了網(wǎng)格自適應(yīng)過程中局部網(wǎng)格的細化和粗化的動態(tài)變化情況。

從圖中可以看出，自適應(yīng)過程中激波區(qū)域的網(wǎng)格加密，而遠離激波區(qū)域在精度低的時候加密了的單元則合并降低網(wǎng)格量。

圖12給出了自適應(yīng)結(jié)束時的網(wǎng)格分布和ε分布。

圖13給出了添加自適應(yīng)后的五階精度(p=4)的等馬赫線圖和表面壓強系數(shù)分布。

圖10 1次和2次自適應(yīng)后的網(wǎng)格Fig.10 Mesh after once and twice adaption

圖11 1次和2次自適應(yīng)后的局部網(wǎng)格Fig.11 Local mesh after once and twice adaption

圖12 自適應(yīng)結(jié)束的網(wǎng)格和ε分布Fig.12 Mesh and ε distribution after adaption

圖13 等馬赫線圖和Cp分布曲線Fig.13 Mach isolines and Cpdistribution

從圖中可以看出，自適應(yīng)后的數(shù)值結(jié)果較原始網(wǎng)格上的數(shù)值精度有顯著提高，并且激波捕捉更精確接近實驗值。

4.3 超聲速NACA0012翼型繞流

為驗證方法的適應(yīng)性，采用上述自適應(yīng)間斷有限元法給出了超聲速(Ma=5.0、α=0)下NACA0012翼型繞流的數(shù)值結(jié)果(p=4)。

初始網(wǎng)格如圖14，初始網(wǎng)格數(shù)為660。自適應(yīng)計算結(jié)束后的馬赫云圖及網(wǎng)格、ε分布如圖15所示。

從圖中可以看出，自適應(yīng)操作主要集中在頭部激波區(qū)域。

圖14 局部網(wǎng)格Fig.14 Local mesh

圖15 馬赫云圖及ε分布Fig.15 Mach contour and ε distribution

5 結(jié)論

(1)將數(shù)值解的高階項貢獻量化后作為自適應(yīng)指示器可行并直觀有效。

(2)即使在很稀疏的初始網(wǎng)格基礎(chǔ)上，使用基于真實物面逼近的網(wǎng)格自適應(yīng)方法也可以獲得高精度的數(shù)值解。

(3)本文設(shè)計的自適應(yīng)方法提高了數(shù)值解精度的同時有效地控制了網(wǎng)格量。

對于非定常問題和三維問題，自適應(yīng)方法在節(jié)省網(wǎng)格工作量和提高計算效率方面作用更為突出。同時具有加密和粗化的功能在非定常計算中會得到更有意義的應(yīng)用，該部分研究正在進行中。在粘性流動的情況下可以在不引入人工粘性項的前提下計算人工粘性項系數(shù)，將其作為網(wǎng)格自適應(yīng)判據(jù)，這部分工作也正在進行。

［1］ Reed W H，Hill T R.Triangular mesh methods for the neutron transport［R］.Los Alamos Scientific Laboratory Report LA-UR-73-479，1973.

［2］ Bassi F，Rebay S.A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations［J］.Journal of Computational Physics，1997，131 (2):267-279.

［3］ Bassi F，Rebay S.High-order accurate discontinuous finite element solution of the 2D Euler equations［J］.Journal of Computational Physics，1997，138(2):251-285.

［4］ Cockburn B，Shu C.The Runge-Kutta discontinuous Galerkin method for conservation laws V:multidimensional systems［J］.Journal of Computational Physics，1998，141(2):199-224.

［5］ Cockburn B，Shu C.The local discontinuous Galerkin method for time-dependent convection-diffusion systems［J］.SIAM Journal on Numerical Analysis，1998，35(6):2440-2463.

［6］ Nguyen N C，Persson P O，Peraire J.RANS solutions using high order discontinuous Galerkin methods［R］.AIAA 2007-0914.

［7］ Lu Hongqiang，Wu Yizhao，Zhou Chunhua，et al.High resolution of subsonic flows on coarse grids［J］.Acta Aeronautica et Astronautica Sinica，2009，30(2):200-204.(in Chinese).

呂宏強，伍貽兆，周春華，等.稀疏非結(jié)構(gòu)網(wǎng)格上的亞音速流高精度數(shù)值模擬［J］.航空學報，2009，30(2):200-204.

［8］ Lu Hongqiang，Zhu Guoxiang，Song Jiangyong，et al.High-order discontinuous Galerkin solution of linearized Euler equations［J］.Chinese Journal of Theoretical and Applied Mechanics，2011，43 (3):621-624.(in Chinese).

呂宏強，朱國祥，宋江勇，等.線化歐拉方程的高階間斷有限元數(shù)值解法研究［J］.力學學報，2011，43(3):621-624.

［9］ Lu H.High-order discontinuous Galerkin solution of low-Re viscous flows［J］.Modern Physics Letters B，2009，23(03):309-312.

［10］Zhang Laiping，Liu Wei，He Lixin，et al.A class of discontinuous Galerkin/finite volume hybrid schemes based on the“static re-construction”and“dynamic re-construction”［J］.Chinese Journal of Theoretical and Applied Mechanics，2010，42(6):1013-1022.(in Chinese).

張來平，劉偉，賀立新，等.基于靜動態(tài)混合重構(gòu)的DG/FV混合格式［J］.力學學報，2010，42(6):1013-1022.

［11］Yu J，Yan C.An artificial diffusivity discontinuous Galerkin scheme for discontinuous flows［J］.Computers＆Fluids，2013，75(20): 56-71.

［12］Yu Jian，Yan Chao.Study on discontinuous Galerkin method for navier-stokes equations［J］.Chinese Journal of Theoretical and Applied Mechanics，2010，42(5):962-969.(in Chinese).

于劍，閻超.Navier-Stokes方程間斷Galerkin有限元方法研究［J］.力學學報，2010，42(5):962-969.

［13］Yang X，James A J，Lowengrub J，et al.An adaptive coupled levelset/volume-of-fluid interface capturing method for unstructured triangular grids［J］.Journal of Computational Physics，2006，217 (2):364-394.

［14］Remacle J F，Li X，Shephard M S，et al.Anisotropic adaptive simulation of transient flows using discontinuous Galerkin methods［J］.International Journal for Numerical Methods in Engineering，2005，62(7):899-923.

［15］Xu Yun，Yu Xijun.Adaptive discontinuous Galerkin methods for hyperbolic conservation laws［J］.Chinese Journal of Computational Physics，2009，26(2):159-168.(in Chinese).

徐云，蔚喜軍.自適應(yīng)間斷有限元方法求解雙曲守恒律方程［J］.計算物理，2009，26(2):159-168.

［16］Wu Di，Yu Xijun.Adaptive discontinuous Galerkin method for euler equations［J］.Chinese Journal of Computational Physics，2010，27 (004):492-500(in Chinese).

吳迪，蔚喜軍.自適應(yīng)間斷有限元方法求解三維歐拉方程［J］.計算物理，2010，27(004):492-500.

［17］VonNeumann J，Richtmyer R D.A method for the numerical calculation of hydrodynamic shocks［J］.Journal of Applied Physics，1950，21(3):232-237.

［18］Xia Y D，Wu Y Z，Lv H Q，et al.Parallel computation of a highorder discontinuous Galerkin method on unstructured grids［J］.Acta Aerodynamica Sinica，2011，29(5):537-541.(in Chinese).

夏軼棟，伍貽兆，呂宏強，等.高階間斷有限元法的并行計算研究［J］.空氣動力學學報，2011，29(5):537-541.

［19］Lu H，Berzins M，Goodyer C E，et al.Adaptive high-order discontinuous Galerkin solution of elastohydrodynamic lubrication point contact problems［J］.Advances in Engineering Software，2012，45 (1):313-324.

Adaptive discontinuous Galerkin method to solve Euler equations based on high-order approximative boundary

Sun Qiang，Lu Hongqiang*，Wu Yizhao
(College of Aerospace Engineering，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，China)

A high-order discontinuous method(DGM)is integrated with adaptive method to solve Euler equations on unstructured mesh.Contribution of the polynomial’s highest-order terms is quantified in the form of artificial viscous coefficient.The coefficient is regarded as the indicator of h-adaptivity.Elements where the coefficients are greater than the upper limit are refined.Those where the coefficients are less than the lower limit are coarsened if they have been refined.A high-order geometric approximation of curved boundaries is adopted to ensure the convergence.Numerical results of test cases are consistent with corresponding experimental ones.High accurate numerical results can be obtained with the h-adaptive method at low expense.

high-order DGM;adaptive method;Euler equations;artificial viscosity;high-order boundary approximation

V211.3

Adoi:10.7638/kqdlxxb-2014.0027

0258-1825(2015)04-0446-08

2014-04-23;

2014-07-08

國家自然科學基金(11272152);航空科學基金(20101552018);江蘇高校優(yōu)勢學科建設(shè)工程資助項目

孫強(1987-)，男，山東萊陽市人，博士研究生，研究方向:計算流體力學.E-mail:submarineultra@126.com

呂宏強*(1977-)，男，博士，副教授，E-mail:hongqiang.lu@nuaa.edu.cn

孫強，呂宏強，伍貽兆.基于高階物面近似的自適應(yīng)間斷有限元法歐拉方程數(shù)值模擬［J］.空氣動力學學報，2015，33(4):446-453.

10.7638/kqdlxxb-2014.0027 Sun Q，Lu H Q，Wu Y Z.Adaptive discontinuous Galerkin method to solve Euler equations based on high-order approximative boundary［J］.Acta Aerodynamica Sinica，2015，33(4):446-453.