基于WENO重構(gòu)的熵穩(wěn)定格式求解淺水方程

程曉晗，聶玉峰，蔡 力

（西北工業(yè)大學(xué)應(yīng)用數(shù)學(xué)系，陜西西安 710072）

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

基于WENO重構(gòu)的熵穩(wěn)定格式求解淺水方程

程曉晗，聶玉峰，蔡 力

（西北工業(yè)大學(xué)應(yīng)用數(shù)學(xué)系，陜西西安 710072）

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

淺水方程；熵穩(wěn)定通量；WENO重構(gòu)

0 引言

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

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

1 數(shù)值方法

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

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

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

式中的Fi＋1/2為數(shù)值通量函數(shù)，其對應(yīng)的熵穩(wěn)定條件變?yōu)槭剑?）的半離散形式.

1.1 熵穩(wěn)定格式

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

那么由相應(yīng)格式得到的解保持總熵不變，滿足離散熵等式，具有二階精度［1］.式（4）有兩個(gè)變量需要確定，卻只有一個(gè)方程.Fjordholm［6］等根據(jù)系統(tǒng)自身特點(diǎn)，構(gòu)造其熵守恒數(shù)值通量為

熵守恒通量式（5）在光滑區(qū)域表現(xiàn)良好，但在激波等間斷區(qū)域時(shí)，會產(chǎn)生偽振蕩，需要添加適當(dāng)?shù)臄?shù)值粘性來保持總熵耗散.將Roe格式的數(shù)值粘性項(xiàng)添加至熵守恒格式［6］，可得滿足離散熵不等式的熵穩(wěn)定數(shù)值通量為

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

1.2 基于WENO重構(gòu)的熵穩(wěn)定格式

式中

非線性權(quán)wk的計(jì)算方法為

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

按照上述步驟關(guān)于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是重力加速度.在空間上進(jìn)行矩形網(wǎng)格劃分，則求解式（9）的守恒型半離散格式為

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

式中，

為了提高格式的數(shù)值精度，分別在x，y方向進(jìn)行五階WENO重構(gòu)獲得，來代替

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

3 數(shù)值算例

在時(shí)間上的推進(jìn)采用三階強(qiáng)穩(wěn)定的Runge?Kutta方法

其中L（·）為空間離散算子，即式（3）或（10）的右端項(xiàng).取CFL條件數(shù)為0.4，下面給出不同算例的計(jì)算結(jié)果.

3.1 大型潰壩問題［6］

在區(qū)域［－1，1］上求解初值問題

取空間網(wǎng)格數(shù)為100，計(jì)算到t＝0.1.精確解由左行稀疏波和右行激波組成.圖1給出了水深的計(jì)算結(jié)果.由圖可知，ESW格式能準(zhǔn)確地捕捉解的結(jié)構(gòu)，比ES格式分辨率高.

3.2 微小水深問題［6］

在區(qū)域［－1，1］上求解初值問題

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

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

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

3.3 圓形潰壩問題［10］

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

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

3.4 部分潰壩問題［11］

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

圖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 結(jié)論

以熵穩(wěn)定數(shù)值通量為基礎(chǔ)，通過在單元交界面處進(jìn)行五階WENO重構(gòu)，建立了一種求解淺水方程的數(shù)值方法.通過數(shù)值模擬可知，該格式具有精度高、捕捉間斷能力強(qiáng)等特點(diǎn).

［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，F(xiàn)eng 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

國家自然科學(xué)基金（11171043，11471261）和西北工業(yè)大學(xué)博士論文創(chuàng)新基金（CX201426）資助項(xiàng)目

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

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