基于WENO重構的熵穩定格式求解淺水方程

程曉晗，聶玉峰，蔡 力

（西北工業大學應用數學系，陜西西安 710072）

文章編號：1001?246X（2015）05?0523?06

基于WENO重構的熵穩定格式求解淺水方程

程曉晗，聶玉峰，蔡 力

（西北工業大學應用數學系，陜西西安 710072）

針對淺水方程，提出一種數值求解格式：空間方向采用滿足熵穩定條件的數值通量，并在單元交界面處進行高階WENO重構，時間上的推進采用強穩定的Runge?Kutta方法.模擬一維和二維經典問題，結果表明，該格式具有分辨率高、基本無振蕩性等特點.

淺水方程；熵穩定通量；WENO重構

0 引言

熵穩定格式是Tadmor等人提出的一種求解雙曲守恒律的數值格式并在近年來受到廣泛關注［1－4］.該格式的構造主要基于熵守恒通量和數值耗散項.熵守恒通量保持總熵不變，在光滑區域表現良好.但在激波等間斷區域，會產生偽振蕩，需要添加適當的數值粘性來保持總熵耗散.與現有的WENO，DG等方法相比，該格式從物理概念出發，滿足守恒律方程的額外條件——熵不等式，所獲得的解是方程具有物理意義的唯一解，能有效避免一些非物理現象的產生，具有廣闊的應用前景.

淺水方程屬于非線性雙曲守恒律，也可采用熵穩定格式求解.文［5］采用了系統本身的粘性作為耗散，計算了二維淺水方程，但需要精細的網格劃分才能得到滿意的數值解.之后，Fjordholm等設計了計算方法簡單的淺水方程熵守恒通量，并采用Roe格式的粘性項作數值耗散［6－7］.本文基于熵穩定數值通量，通過在單元交界面處進行高階WENO重構［8］，時間上的推進采用三階強穩定的Runge?Kutta格式［9］，給出了一種求解淺水方程的高分辨率數值方法.最后，將建立的格式應用于一維和二維的典型問題來驗證方法的有效性.

1 數值方法

考慮一維無粘淺水方程，其守恒形式為

式中U＝（h，hu）T是守恒型向量，F（U）＝（hu，hu2＋gh2/2）T是通量．h是水深，u是平均流速，g是重力加速度．給定熵對（E，H），其中E（U）是熵函數，其Hessian矩陣E″（U）對稱正定，H（U）滿足H′（U）T＝E′（U）TF′（U）是熵通量函數，則式（1）的具有物理意義的唯一弱解需要滿足熵不等式

在空間方向采用等距網格劃分，求解式（1）的守恒型半離散格式為

式中的Fi＋1/2為數值通量函數，其對應的熵穩定條件變為式（2）的半離散形式.

1.1 熵穩定格式

取熵函數為系統的總能，即E（U）＝（hu2＋gh2）/2，則可得相應熵通量函數為H（U）＝hu3/2＋guh2．定義熵變量為V＝?UE＝（gh－u2/2，u）T，熵勢為Ψ＝VTF－H＝guh2/2，則如數值通量滿足

那么由相應格式得到的解保持總熵不變，滿足離散熵等式，具有二階精度［1］.式（4）有兩個變量需要確定，卻只有一個方程.Fjordholm［6］等根據系統自身特點，構造其熵守恒數值通量為

熵守恒通量式（5）在光滑區域表現良好，但在激波等間斷區域時，會產生偽振蕩，需要添加適當的數值粘性來保持總熵耗散.將Roe格式的數值粘性項添加至熵守恒格式［6］，可得滿足離散熵不等式的熵穩定數值通量為

式中的R，Λ，Vi＋1，Vi分別為

1.2 基于WENO重構的熵穩定格式

式中

非線性權wk的計算方法為

其中d1＝3/10，d2＝3/5，d3＝1/10，ε＞0的引入是為了防止分母為零，一般取ε＝10－6，光滑因子βk定義為

按照上述步驟關于xi＋1/2作對稱即可獲得

2 二維推廣

考慮二維淺水方程

其中守恒型向量U＝（h，hu，hv）T，通量分別為F（U）＝（hu，hu2＋gh2/2，huv）T和G（U）＝（hv，huv，hv2＋gh2/2）T．h是水深，（u，v）分別是x，y方向的平均流速，g是重力加速度.在空間上進行矩形網格劃分，則求解式（9）的守恒型半離散格式為

取熵函數E（U）＝（hu2＋hv2＋gh2）/2和熵通量函數H（U）＝（hu3＋huv2）/2＋guh2，K（U）＝（hv3＋hu2v）/2＋gvh2，則熵變量和對應的熵勢分別為V＝（gh－（u2＋v2）/2，u，v）T，Ψ＝guh2/2，Φ＝gvh2/2．類似于一維情況的處理方式，可得求解二維淺水方程的熵穩定數值通量為

式中，

為了提高格式的數值精度，分別在x，y方向進行五階WENO重構獲得，來代替

Ui＋1/2±1/2，j，Ui，j＋1/2±1/2，于是可得高精度數值通量為

3 數值算例

在時間上的推進采用三階強穩定的Runge?Kutta方法

其中L（·）為空間離散算子，即式（3）或（10）的右端項.取CFL條件數為0.4，下面給出不同算例的計算結果.

3.1 大型潰壩問題［6］

在區域［－1，1］上求解初值問題

取空間網格數為100，計算到t＝0.1.精確解由左行稀疏波和右行激波組成.圖1給出了水深的計算結果.由圖可知，ESW格式能準確地捕捉解的結構，比ES格式分辨率高.

3.2 微小水深問題［6］

在區域［－1，1］上求解初值問題

取空間網格數為100，計算到t＝0.1.初始速度使得水流向兩邊擴展，導致中間區域水深很小，接近干床（水深為0），許多格式會產生負水深從而使得計算崩潰.由圖2可知，ESW格式能成功地模擬微小水深問題，且對左右稀疏波的計算更為精確.

圖1 大型潰壩問題水深Fig.1 Water heights for the big dam break problem

圖2 微小水深問題水深Fig.2 Water heights for the expansion problem

3.3 圓形潰壩問題［10］

在區域［0，50］×［0，50］上求解初值問題

取空間網格數為100×100，計算到t＝0.69.精確解由向內稀疏波和向外激波組成.圖3依次給出了水流曲面圖、水深等值線圖和y＝25時的水深截面圖.結果顯示，ESW能準確地計算圓形潰壩問題，等值線圖中圓的對稱性保持得很好.

3.4 部分潰壩問題［11］

該問題由Fennema和Chaudhry首次提出.考慮一個正方形（200×200）平底水庫，上下游為靜水，水深分別為10和5.上下游邊界采用入流和出流條件，其余部分采用反射邊界條件.水壩從95到170之間自初始時刻起瞬時潰破，取空間網格數為40×40，計算到t＝7.2.圖4依次給出了水流曲面圖和水深等值線圖.可以看出，該結果與文［11］中的結果保持一致，ESW格式能準確模擬潰壩現象.

圖3 圓形潰壩問題的水流曲面、水深等值線及截面Fig.3 Surface，contours and slice plots ofwater height for the circular dam break problem

圖4 圓形潰壩問題的水流曲面和水深等值線Fig.4 Surface and contour plots of water height for the partial dam break problem

4 結論

以熵穩定數值通量為基礎，通過在單元交界面處進行五階WENO重構，建立了一種求解淺水方程的數值方法.通過數值模擬可知，該格式具有精度高、捕捉間斷能力強等特點.

［1］ Tadmor E.The numerical viscosity of entropy stable schemes for systems of conservation laws［J］.Math Comp，1987，49 （179）：91－103.

［2］ Tadmor E.Entropy stability theory for difference approximations of nonlinear conservation laws and related time?dependent problems［J］.Acta Numer，2003，12（1）：451－512.

［3］ Luo Li，Feng Jianhu，Tang Xiaojuan，Xiang Liang.High resolution entropy stable schemes for hyperbolic conservation laws ［J］.Chinese JComput Phys，2010，27（5）：671－678.

［4］ Hiltebrand A，Mishra S.Entropy stable shock capturing space?time discontinuous Galerkin schemes for systems of conservation laws［J］.Numer Math，2014，126（1）：103－151.

［5］ Tadmor E，ZhongW G.Energy?preserving and stable approximations for the two?dimensional shallow water equations［C］∥Mathematics and Computation，a Contemporary View，Proc of the Third Abel Symposium.Berlin：Springer，2008：67－94.

［6］ Fjordholm U S，Mishra S，Tadmor E.Energy preserving and energy stable schemes for the shallow water equations［C］∥Foundations of Computational Mathematics，Hong Kong 2008.London Mathematical Society Letures Note Series（363）. London：Cambridge University Press，2009：93－139.

［7］ Fjordholm U S，Mishra S，Tadmor E.Well?balanced and energy stable schemes for the shallow water equations withdiscontinuous topography［J］.JComput Phys，2011，230（14）：5587－5609.

［8］ Jiang G S，Shu CW.Efficient implementation of weighted ENO schemes［J］.JComput Phys，1996，126（1）：202－228.

［9］ Gottlieb S，Shu CW，Tadmor E.Strong stability?preserving high?order time discretizationmethods［J］.SIAM Rev，2001，43 （1）：89－112.

［10］ Chen Jianzhong，Shi Zhongke.Study on semidiscrete central?upwind scheme for the2D shallow water equations［J］.Advances in Water Science，2005，16（6）：853－857.

［11］ Fennema R J，Chaudhry M H.Explicitmethods for2?D transient free surface flows［J］.JHydraul Eng，1990，116（8）：1013 －1034.

WENO Based Entropy Stable Scheme for Shallow W ater Equations

CHENG Xiaohan，NIE Yufeng，CAILi （Department of Applied Mathematics，Northwestern Polytechnical University，Xi′an 710072，China）

A high resolution scheme is presented for shallow water equations.The scheme is based on entropy stable numerical flux with high order weighted essentially non?oscillatory（WENO）reconstruction at cell interfaces.A strong stability?preserving Runge?Kuttamethod is employed to advance in time.Severalbenchmark numerical examples demonstrate that the scheme is accurate and non?oscillatory.

shallow water equations；entropy stable flux；WENO reconstruction

TV131.4；O242

A

2014－08－26；

2014－12－01

國家自然科學基金（11171043，11471261）和西北工業大學博士論文創新基金（CX201426）資助項目

程曉晗（1987－），男，安徽樅陽，博士生，主要從事計算流體力學的研究，E?mail：chengxh168＠163.com

Received date： 2014－08－26；Revised date： 2014－12－01