基于LU-SGS迭代的DGM隱式方法研究

郝海兵, 張強, 楊永, 梁益華

間斷Galerkin方法[1](DGM)由于具備易于實現高精度以及求解間斷問題等優點而受到廣泛關注。該方法在任意單元內部通過提高逼近多項式的階次實現高精度;空間離散時允許單元邊界上存在間斷,具有靈活處理復雜區域邊值問題的能力;在單元邊界處可以采用有限體積法中的思路,基于Riemann問題來構造數值通量函數,實現逆風格式,從而更有利于求解間斷問題。目前,國外對該方法已經展開了近十多年的研究,并取得了較大的進展[2],而國內相關研究起步較晚,在近幾年才逐漸受到重視。國內主要進展包括邱建賢等[3]將WENO格式作為限制器引入到DGM中;張來平等[4]提出了靜動態混合重構的DG/FV混合格式等。當前DGM主要應用于空間離散,而時間方向推進多沿用了Shu等[2]構造的TVD-RKDG方法。因此,當求解更高精度或大規模計算問題時,穩定性條件變得越來越苛刻,從而導致當地時間步長取值越來越小,收斂率會急劇下降,計算非常耗時。針對該問題,目前主要采用p型多重網格方法[2]和隱式時間離散2種加速收斂技術來提高計算效率。本文主要對隱式格式離散進行研究,p型多重網格方法見參考文獻[2,5]。目前,國外很多學者對間斷Galerkin有限元隱式算法進行了多種試探性的研究,不過目前還沒有形成成熟穩定的算法。主要進展包括如下:Bassi和Rebay首次將GMRES隱式算法應用于DG求解二維可壓縮N-S方程[6];Hartmann等人使用GMRES-Newton算法求解二維Euler和N-S方程[7]等。

本文基于非結構網格,針對定常Euler方程組,發展相應的DGM隱式算法?？紤]到LU-SGS迭代法不需要存儲和處理大型矩陣,并且非常適合于非結構網格的無序性,具有顯著的計算效率,并且在有限體積法中得到了廣泛應用。本文借鑒其思想,將該方法推廣到高精度DGM隱式格式中,構造適合于DGM隱式求解的LU-SGS迭代格式。通過對NACA0012翼型和ONERA M6機翼跨聲速無粘流動進行數值模擬驗證其加速收斂性能。

1 控制方程

考慮非定常Euler方程在直角坐標系的守恒形式:

(1)

式中:

ρ、e、p分別表示氣體的密度、單位體積總內能及壓力,ui表示xi方向的速度。對于氣體動力學方程。

γ=1.4是比熱比。

2 間斷Galerkin方法

間斷Galerkin方法求解方程組(1),首先需要將計算區域劃分成互不重疊的子域。子域可以選取為任意形狀,對于二維空間,本文采用三角形非結構網格,三維空間,采用四面體非結構網格。定義有限元空間Vh:

Vh={vh∈L2(Ω):vh|K∈V(K);?K∈τh}

式中：τh為子域空間,V(K)為局部函數空間,取作p(p=1,2,…)次多項式的集合。

假設間斷函數在有限元空間中的近似解為Uh,將Euler方程(1)兩邊同乘以試驗函數v,寫成變分形式,再由Green公式,得到守恒方程組的弱形式:

(2)

式中:Ω為計算域;?Ω為Ω的邊界。在每個單元內:

(3)

式中:φj(x)是基函數。將(3)式帶入(2)式,得到半離散形式的守恒方程。

(4)

(4)式稱為p階間斷Galerkin有限元離散,p為(3)式中基函數的最大階數。不難看出Uh、vh在整個計算區域上不再連續,在單元邊界上是間斷的,而且整個流場的自由度已經轉變成求解插值系數Ui,不再是流場守恒變量。

為了便于計算,通常還需要進行坐標變換。針對三角型單元,我們采用面積坐標,針對四面體單元,采用體積坐標。將總體笛卡爾坐標轉換成局部的自然坐標,將空間中任意單元轉變成局部坐標系中的標準單元。(4)式變為:

(5)

式中：|J|為體坐標變換雅可比矩陣,|Js|面坐標變換雅可比矩陣。

對于(5)式中的積分,我們選用高斯數值積分。

3 隱式時間推進格式

TVD-RKDG時間推進格式中,由于受穩定性條件限制,計算中CFL數過小,導致效率比較低,而隱式推進格式一般是無條件穩定的,計算步長可以取得較大,從而整體效率較高。本文在借鑒非結構網格有限體積LU-SGS隱式計算格式的基礎上,構造出了適用于間斷Galerkin有限元法的LU-SGS隱式計算格式。

LU-SGS方法基本思想是運用通量線化假設和最大特征值方法進行雅可比矩陣分裂,把塊對角矩陣分解為上、下2個三角矩陣,這種分解可以避免繁雜的矩陣求逆運算,極大地提高了計算效率。該方法本質上仍然是一種近似因式分解方法,對于結構網格比較適合,而在非結構網格下,由于網格存儲的不規則性,該方法的計算效率要差一些。

對公式(5)在時間方向上采用一階向后歐拉積分:

(6)

將其展開:

(7)

由于方程(6)為一非線性系統,直接求解非常麻煩,對(7)式中的二、三項分別進行一階泰勒展開,進行線化處理:

(8)

式中:

(9)

為了簡化處理,本文將方程(6)中時間項和空間項的數值通量分開單獨處理,在時間項上無粘通量Fnum采用LF數值通量格式,可以分裂為

(10)

式中：λij為雅克比矩陣的譜半徑。

則(9)式可以簡化成:

將公式(7)、(8)式帶入(6)式,省略高階項:

(11)

式中:

將方程(6)代入上式,得到等價的矩陣形式

(12)

(13)

(14)

式中：L為嚴格下三角矩陣,D為純對角線矩陣,U為嚴格上三角矩陣

(15)

將(15)式變形為

(D+L)D-1(D+U)Δu=Rn+(LD-1U)Δu

(16)

忽略上式右端第二項高階小量,可得:

(D+L)D-1(D+U)Δu=Rn

(17)

上式只需要通過向前掃描和向后掃描兩個過程即可求解自由度,其中

向前掃描:

(18)

向后掃描:

(19)

具體實現過程為:

(20)

(21)

4 算例與分析

為了驗證本文基于LU-SGS迭代法構造的DGM隱式算法的正確性和有效性,本文分別對繞NACA0012翼型和ONERA M6機翼的跨聲速無粘流動進行數值模擬,并和TVD-RKDG計算結果進行比較,來驗證隱式算法的加速收斂效果和精度。本文算例中使用的計算機基本配置CPU為酷睿I5 3.2GHz、內存為4G,所有網格均采用Delaunay方法生成。

NACA0012翼型的流動計算狀態為Ma∞=0.8,攻角α=1.25°。整個流場包括3 531個網格節點,6 789個網格單元。圖1給出了翼型表面壓力系數分布比較,從圖中可以看出,采用DGM隱式算法并沒有改變精度,并且翼型上下表面的激波均捕捉的很好。圖2分別給出了關于CPU時間和迭代步數的密度最大殘值收斂歷程比較,當采用TVD-RKDG時計算花費了116 min左右,采用隱式算法后,計算花費了14 min左右,計算效率提高8倍左右。

圖1 翼型表面壓力系數分布比較圖 圖2 殘值收斂歷程比較

ONERA M6機翼的流動計算狀態為Ma∞=0.84,攻角α=3.06°。整個流場包含65 115個網格節點和459 353個網格單元。圖3顯示了機翼上表面壓力等直線圖,圖中λ狀激波結構捕捉非常清晰,外弦和內弦激波大約在展長87%處相交,在94%處又分開,和實驗結果[15]基本吻合。圖4給出了M6機翼65%展長處的剖面壓力系數分布比較,結果和實驗值吻合較好,并且DG隱式算法和TVD-RKDG得到的壓力系數分布完全重合。圖5分別給出了關于CPU時間和密度最大殘值收斂歷程比較,從圖中可見,當采用TVD-RKDG時計算花費了2 300 min左右,采用隱式算法后,計算花費了500 min左右,計算效率提高5倍左右。

圖3 M6機翼上表面壓力等值線

圖4 65 %展向位置處剖面壓力系數分布的對比 圖5 殘值收斂歷程比較

5 結 論

本文主要構造了適合于DGM隱式求解的LU-SGS迭代格式,并發展了一套高效、實用的非結構網格Euler方程求解程序。通過數值模擬跨聲速Euler方程中來驗證其效率和精度。數值算例表明:采用DGM隱式算法能夠很好的維持DGM的精度,并能明顯的提高流場的收斂速度以及降低計算成本。在本文的二維算例中,計算效率達到8倍左右,三維算例中,計算效率達到5倍左右。雖然,本文的工作是基于DG(P1),但由于公式推導具有普遍性,可以比較容易推廣到更高精度。因此,今后工作還需要進一步考察本文構造的隱式算法在更高精度DGM中的性能。

參考文獻：

[1] Reed N H, Hill T R. Triangle Mesh Methods for the Neutron Transport Equation[R]. Los Almos Scientific Laboratory, Report LA-UR-73-479, 1973

[2] Fidkowski K J, Oliver T A, Lu J, Darmofal D L. P-Multigrid Solution of High-Order Discontinuous Galerkin Discretizations of the Compressible Navier-Stokes Equations[J]. J Comp Phys, 2005, 207: 92-113

[3] Qiu J, Shu C W. Hermite WENO Schemes and Their Application as Limiters for Runge-Kutta Discontinuous Galerkin Method: One Dimensional Case [J]. J Comp Phys, 2004, 193: 115-135

[4] Zhang L P, Liu W, He L X, Deng X G. A New Class of DG/FV Hybrid Schemes for One-Dimensional Conservation Law[C]∥The 8th Asian Conference on Computational Fluid Dynamics, Hong Kong, 2010: 10-14

[5] 郝海兵,楊永. 非結構網格上P型多重網格法流場數值模擬[J]. 計算力學學報, 2011, 28(3): 360-365

Hao Haibing, Yang Yong. The Research of P-Multigrid Solution for Discontinuous Galerkin Method [J]. Chinese Journal of Computational Mechanics, 2011, 28(3): 360-365 (in Chinese)

[6] Bassi F, Rebay S. GMRES for Discontinuous Galerkin Solution of the Compressible Navier-Stokes Equations∥Cockburn B, Karniadakis G E, Shu C W, Discontinuous Galerkin Method: Theory, Computations and Applications[M], Springer-Verlag, 2000

[7] Hartmann R, Houston P. Symmetric Interior Penalty DG Methods for the Compressible Navier-Stokes Equations I: Method Formulation[J]. Internatimal Journal of Numerical Analysis and Modeling, 2006(3): 1-20