基于不同數值通量的間斷Galerkin方法性能的研究

李 廣，張小華

(三峽大學 理學院，湖北 宜昌 443002)

淺水波方程已經被廣泛應用于海洋和水利工程以及地球物理流，如河流、潮汐、流動的水庫和明渠等[1]，該系統將其描述為帶有附加源項的守恒律.然而對于雙曲守恒律，即使具有光滑的初始條件和邊界條件，也可能在計算域內產生不連續解[2].因此，淺水波方程的數值模型需要能夠捕獲這些不連續性，DG方法是一種解決這類問題的適當選擇.它具有精度高、并行效率高和易于處理邊界條件等優點，在計算流體力學的框架下，已經被廣泛地應用于空氣聲學、粒狀流、磁流體動力學、氣象學、海洋學、湍流、粘彈性流動、淺水建模以及天氣預報等領域.

最早的DG方法是由 Reed和Hill于1973年在一個報告中提出來的，并解出了與時間無關的中子輸運線性雙曲方程[3].Cockburn和舒其望等人對其發展作出了非常重要的貢獻，即建立了一個易于求解且與時間相關的非線性雙曲守恒律的框架[4-10].在這個框架中，時間離散使用的是高階龍格-庫塔(RK)格式，空間離散來自于有限元方法，數值通量一般使用的是Lax-Friedrichs通量，并且在時間步進格式中的每一步后都使用一次全局變分有界(TVB)限制器.

DG方法的一個重要組成部分是基于精確或近似黎曼求解器的數值通量，同時借鑒了有限差分法和有限體積法.一般情況下，通量的定義不是唯一的，取而代之的是一個數值通量.因此，必須選取一種正確的解，我們把這個解稱為數值通量[11].通量的作用是運用所求方程的信息流動來保證所構造的DG格式的穩定性，數值通量的選取也是構造各種DG格式的關鍵，因此必須保證數值通量是相容的.關于數值通量的另一要求是選取的通量必須是單調的，即既要對第一變量非遞減，又要對第二變量非遞增[12].在大多數研究DG的文獻中，Lax-Friedrichs 通量因其簡單性而被廣泛使用.本文主要研究使用基于各種不同數值通量的DG法在求解一維淺水波方程中的性能，比如Lax-Friedrichs通量[12-13]、Rusanov通量[14]、Harten-Lax-van Leer通量[2，15]以及Roe通量[2].

1 一維控制方程

在一維空間中考慮淺水波方程，其形式如下：

(1)

其中，h為水流深度，u為x方向上平均水深的速度，g是重力加速度，S0與Sf分別是底部床坡與斜坡摩擦，下標t和x分別表示關于時間t和坐標x的偏導數.該系統考慮了由于底部地形和斜坡摩擦引起的源項，還可以加入其他源項，比如粘性應力、科氏力、潮汐勢力和風表面應力等影響.先將方程(1)改寫成如下守恒形式為

其中，U=[h，hu]T為保守變量的向量，F=[f1，f2]T=[hu，hu2+0.5gh2]T是通量向量，源項向量是S=[0，gh(S0-Sf)]T，上標T表示矩陣的轉置.一般情況下，假定穩態流動的斜坡摩擦在非定常流態下是有效的，其中n為曼寧粗糙度系數，于是底部與斜坡摩擦項可以分別近似為

2 DG方法

2.1 空間離散

(2)

對K個單元都成立， 但是方程(2)還不能用于求解整體解. 假設試驗函數光滑但跨界面不連續，對方程(2)關于空間變量分部積分，所得格式就是經典弱形式下的節點DG方法

(3)

2.2 數值通量

2.2.1 Lax-Friedrichs (LF)通量LF通量是DG方法和高階有限體積法最簡單、應用最廣泛的基本部分之一.然而，LF通量的數值黏性在標量問題的單調通量中也是最大的.其定義和最大波速度分別為

2.2.2 Rusanov(RUS)通量定義為

2.2.3 Harten-Lax-van Leer (HLL)通量Harten，Lax和Van Leer于1983年提出了HLL黎曼求解器，該求解器具有三個恒定狀態，分別被兩個波隔開.淺水波方程的HLL通量為

2.2.4 ROE通量一維矩形河道淺水流動的Roe通量函數為

Roe通量函數中的變量為

2.3 時間離散

時間步長被記為Δt，且Δt=tn+1-tn，則Un+1=U(t=tn+1).此處采用文獻[13]中的三階三步總變差遞減-龍格庫塔格式，也被稱為保強穩定的龍格庫塔格式.這種格式既能夠保證高精度又能保證在時間積分過程中不會引進額外的振蕩.因此，適用于淺水波方程的具有最優三階精度的SSP-RK格式為

2.4 斜率限制器

即使在初始條件是光滑的情況下，非線性守恒律的解也會出現不連續甚至產生振蕩，這給數值求解帶來了很大的困難，所以在SSP-RK時間步進格式中每一步之后都使用一次斜率限制器是一種很好的選擇.限制器既能夠獲得高效、可靠的重建策略，又能夠獲得更好的數值解，還能節省計算成本.雖然可以在DG方法中使用多種限制器，但是對各種不同的問題，沒有一種確切的限制器明顯優于其他限制器[11].

2.4.1 Minmod梯度限制器(MD)梯度限制器重構的方法是基于守恒定律的單調上游格式和分段拋物線法.為尋求一個具體修改解的局部均值以保證TVDM性質，定義minmod函數m(·)為

如果k個變量取為k個單元上解的梯度，那么當所有梯度不同號時，minmod函數將梯度設置為零，這表明一個振蕩；否則，minmod函數取為最小的梯度.

2.4.2 改進的梯度限制器(MM)簡單利用梯度限制器可能會破壞光滑區域上的高階精度，而廣義梯度限制器雖然可以改進在解的光滑區域處的精度，但是不能克服在局部極值附近精度的損失.本文使用改進的一種minmod梯度限制器.其處理方法是放松全局變分的衰減條件，只要求平均的全局變分有界，即TVBM條件.對此，需要稍微修改minmod函數m(·)的定義為

則滿足TVBM條件.常數M是在局部極值處二階導數的一個上界，一般取M=20 .

3 數值實驗

表1 DG方法計算水深h所得的L1， L2和L∞誤差

表1分別提供了基于四種不同通量的DG格式關于水深h的L1，L2和L∞誤差.將通量為“X”的DG方案表示為DG—X，例如通量為LF的DG方案表示為DG—LF.當使用梯度限制器時，其他方案關于水深h的L1，L2和L∞誤差分別約是DG—LF方案的95%，96%，95%.當使用改進的梯度限制器時，其他方案關于水深h的L1，L2和L∞分別約是DG—LF方案的91.3%，88.7%，86.2%.相比之下，當限制器相同時，使用ROE通量計算所得的誤差在所有方案中最小.

可以看出，使用改進的梯度限制器計算出了水深h與水流速度u分別在DG—LF，DG—RUS， DG—ROE和DG—HLL方案下同一網格上的數值解(包含接觸不連續).可以看出，由DG—RUS， DG—HLL和DG—ROE方案計算的結果略優于DG—LF方案.關于不連續的解，以上幾種方案的計算效果類似.

4 結語

本文系統地研究和比較了DG方法的四種不同數值通量，并對一維淺水波方程進行了模擬，結果表明ROE通量和改進的梯度限制器可能是DG方法的較好選擇.本文僅對線性三角形單元和一維平底淺水波方程進行了測試，擴展到非平底河床和高維淺水波方程的研究，以及與現有模型相結合，將是下一步的工作.