周超红，姚清河，朱庆勇

(中山大学应用力学与工程系，广东 广州 510275)

浅水波方程是计算流体力学中的一类重要方程。浅水波方程的应用范围十分广泛，可用于河道流量、溃坝决堤、污染物输运扩散、河道输沙、洪水预报、河口潮汐、涌浪等大量问题。

由于浅水波方程的应用背景十分广泛，其数值方法的研究自然成为了计算流休力学和计算水动力学的重要课题。由于二维浅水波方程的一般形式比较复杂，本文中使用在无黏性、无摩擦、河底无倾斜的理想情况下的二维浅水波方程作为求解模型。群速度控制的思想是由傅德薰、马延文[1]提出的，朱庆勇[2-4]利用这种思想成功构造了求解空气动力学方程的高分辨率格式，之后他们[5]将这一思想引入浅水波方程，构造了求解一维浅水波方程的群速度控制格式。

利用群速度控制的思想，本文构造了多种求解二维浅水波方程的高分辨率格式。通过引入色散度量以及各向异性度量的概念，为评估差分算子对微分算子的逼近程度提供了客观标准。数值实验结果表明，这些格式均可以很好地控制非物理振荡并保证格式精度。

1 预备知识

考虑单波方程及其半离散方程

,f=f(u)

(1)

(2)

这里Fj/Δx是一阶导数∂f/∂x的差分逼近，定义差分算子如下

±(fj±1-fj)

Fj可以表示为

βm(fj + m + 1-fj + m)

下面我们利用傅里叶分析法对差分方法进行理论分析。

考虑f(u)=au，其中a为常数，初值条件u(x,0)=u0(x)，考虑以2π为周期的单波解，对于半离散方程(2)，我们设x方向网格步长为

设其初始条件为uj=exp(ikxj)，Fj=akeexp(ikxj)，ke=kr+iki。其中

βm{cos[k(m+1)Δx]-cos[kmΔx]}，

则半离散方程(2)的解是

(3)

而该问题的精确解为

uj(t)=exp[ik(xj-at)]

(4)

对比(3)和(4)，可知对于任意的k,kΔx<π，差分格式的相容性要求有kr/Δx→0,ki/kΔx→1，这里的akr和ki分别决定了差分格式的数值耗散与数值色散性质。akr>0说明格式是正耗散的，akr<0说明格式是负耗散的，而akr=0则格式是非耗散的，负耗散的格式是不稳定的。数值耗散主要影响振幅的大小,而数值色散效应则是产生虚假振荡、引起波传播速度不均一的根本原因。高分辨率算法首先要求能够抑制色散效应引起的振荡，其次也要尽量提高数值解的精度。

波的群速度，是指波振幅外形上的变化在空间中传播的速度。根据以上的讨论，我们可定义群速度为sg=dki/dα，其中α=kΔx。对于(1)的精确解(3)，ki=α，sg=1。

2 群速度修正格式的构造

在对群速度控制法的研究中我们发现，文[6]中构造三阶迎风紧致群速度控格式

(5)

其实就是针对三阶迎风紧致格式，引入了如下的群速度修正项

对半离散方程(2)进行修正。参考上式，我们考虑群速度修正项

其中ssj=sign[(uj+1-uj)(uj+2-uj+1-uj+uj+2)]，称为激波结构函数。φj满足在振荡区域φj=1，平滑区域φj=0，称为开关函数。

可取φj=

0<ω<1,κ=0.05

通过分析Lg我们发现，在平滑区域，Lg=0；在激波前后，Lg能使得振荡波向激波靠拢；在极值点，Lg只有耗散而无色散，有利于振荡的消除；而在拐点，Lg呈现负耗散性质，则有利于激波形状的保持。根据Lg的这些性质，我们得知，Lg其实可以对任意差分格式进行修正，而不仅仅限于三阶迎风紧致格式。

现在我们从算子修正的角度来看待方程(2)，即采用如下半离散格式

(6)

其中Fj为基本算子， 0<σ<1。

我们改称式(6)中的Lg为群速度修正算子，称σ为群速度修正系数。该格式与格式(5)不仅仅是在形式上有所区别，两者的意义完全不一样。前者表示群速度修正算子度基本算子的修正，而后者表示群速度修正项对半离散格式(2)的修正。利用算子修正的观点来构造群速度控制格式，有利于对构造的格式进行分析。

对于半离散格式(6)，使用不同的基本算子和群速度修正系数σ，我们就可以构造出各种各样的群速度控制格式。我们将使用格式(6)构造群速度控制格式的方法，称为群速度修正法。

图1 波数比例因子与色散显现因子

图2 三阶迎风群速度控制格式的各向异性效应

下面以三阶迎风群速度控制格式为例，利用色散度量函数，寻找最佳的群速度修正参数。由于三阶迎风格式的振荡主要在激波后，我们考虑，激波后三阶迎风群速度控制格式的差分算子

σ(-fj+1+3fj-3fj-1+fj-2)，

σ(-4cos2α+2cosα+2)

给定不同的σ，将sg代入Dmf右端，对Dmf的右端求数值积分(可取αn=nπ/10，n=0,1,…,10)，再画出σ-Dmf的曲线，Dmf最低点对应的σ，即为最优的修正参数。由于群速度控制格式的耗散随着σ的增大而增大，因而若最低点偏左侧点的纵坐标与最低点相差不大，最优的修正参数则应取左侧点对应的σ。

图3 群速度控制算子的色散度量

对于二维问题，我们还需要考虑算子的各项异性效应。考虑二维单波方程

，a,b为常量

(7)

初值u(x,y,0)=exp(iK·X)，其中K=[k1,k2]T，X=[x,y]T，k1,k2分别为x,y方向上的波数。定义

并将它重写为l=[cosθ,sinθ]T，这里的θ就是方位角。该问题的精确解可以表示为

(8)

方程(7)的半离散格式为

(9)

其中δxuij/Δx、δyuij/Δy分别为∂u/∂x、∂u/∂x的差分逼近。方程(9)的精确解可以表示为

(10)

其中

，，

α=k1Δx，β=k2Δy

对‖l‖的逼近程度。

为度量差分算子的各项异性，仿照前面色散度量的定义，定义

τθdω

为各项异性度量。其中τ= ecos(ω)-e-1称为波成分比例因子，表示数值解中ω与对应的波的比例。ν=1/σ4/3称为各向异性显现因子，用于用于度量耗散对各向异性的压制力度。

图4 群速度控制算子的各向异性度量

根据本节提供的方法，我们可以得到各种高阶群速度控制法的最佳群速度修正参数，具体为表1。

表1 基本算子最佳群速度修正参数对应表

由于开关函数的存在，平直区域将不进行群速度控制，因而只需考虑主要振荡区域对应的最佳群速度修正参数的值。对于一维问题只需考虑色散度量的值。例如，三阶迎风算子的振荡主要出现在激波后，因而推荐取值为激波后色散度量下的最佳值0.20。对于二维问题，要同时考虑色散度量和各向异性度量。如四阶迎风算子的振荡主要出现在激波前，色散度量最佳值为0.20，各向异性度量下的最佳值为0.10，推荐取两者的平均值0.15。对于三阶迎风紧致算子，实验发现σ取0.05时色散比较明显，故推荐值取为0.10。

3 二维浅水波方程的群速度修正格式

二维浅水波方程在无黏性、无摩擦、河底无倾斜的理想情况下，其数学模型的齐次守恒形式如下

(11)

其中

，，

h为水深，u为x方向的速度大小，v为y方向的速度大小，g为重力加速度。

为适应一般物型计算。引入随时间变化的贴体曲线坐标变换

τ=t,ξ=ξ(t,x,y),η=η(t,x,y)

则由方程(11)可导出二维曲线坐标系下的守恒型浅水波方程

(12)

其中

V=JU=(v1,v2,v3)T,

对于双曲型方程，为了保证计算我们常常需要使用迎风方法，但迎风方法要求对流通矢量F*、G*进行分裂[7]。文[8]提出了针对二维浅水波方程的有效分裂方法，参照该文中的方法，我们给出曲线坐标下二维浅水波方程流通矢量的分裂格式。

下边构造曲面坐标下二维浅水波方程的群速度修正格式。记ωi,j=ω(iΔξ,jΔη)，定义算子

ωi,j=(ωi+1,j-ωi-1,j)/2,

并令

其中φξ、φη为ξ方向、η方向的开关函数，λξ(i,j)、λη(i,j)分别为A*、B*特征值中绝对值最大的特征值。假设δξ、λξ分别为ξ方向、η方向的基本算子，则曲坐标系下二维浅水波方程的群速度修正格式为

其中为群速度修正参数。若基本算子是迎风算子，则上式中

,

ωi+1,j+ 3ωi,j-6ωi-1,j+ωi-2,j),

若基本算子取为四阶中心差分算子，δξ、δη满足

确定基本算子及其对应的群速度修正系数后，我们可以利用三阶龙格库塔法求解半离散方程(13)，进而通过U=V/J,求得U。

4 数值算例

算例1 倾斜水跃问题

一条40 m长的平底河道，上游入口宽度为30 m，来流Froude数Fr=2.5，在河道一侧10 m处开始以角度α=5°收缩,收缩的边与超临界流相互作用形成附体涌波。入口处采用超临界流固定边界条件,出口处采用传输边界条件,河道两边采用滑移无穿透边界条件。计算网格数为120×90。此算例有精确解，涌波线与X轴的夹角为β=28.32°，波后水深与波前水深比为h2/h1=1.250。

本算例算法使用四阶群速度控制法，σ=0.2

当流体达到稳定状态时,据图5 可得到底端激波线与出口相交处点的坐标为(40,13.46,1.001),据此可算得波角β=28.87°,数值结果很好地模拟了附体涌波结构。波后水深与波前水深比为h2/h1=1.250～1.251。采用TVD方法[9]得到的结果为β=29.01°，h2/h1=1.246～1.254，显然本文使用的四阶群速度控制法结果比较准确。

图5 四阶群速度控制法

算例2 二维圆形溃坝问题

一个50 m×50 m方形区域，区域被一圆心在区域中心(25,25)，半径r=11 m的圆柱面形坝分为两部分，坝内水深10 m，坝外水深1 m。假定坝体突然溃破，计算网格数取50×50。

本算例算法使用三阶迎风群速度控制法σ=0.15。图6和图7给出了t=0.7 s时水波运动的状态。二维圆形溃坝问题的真实解应当满足在任意时刻水深等值线是一系列以(25，25)为中心的圆。图7(b)中所示的水深等值线与标准圆相当接近，说明三阶迎风群速度控制格式能够有效地对各个方向激波进行群速度控制，各向异性效应极小，能很好地模拟二维浅水波的运动。

图6 水深曲面

图7 三阶迎风群速度控制法

5 结论

本文得到了几种常用高阶差分算子在色散度量和各向异性度量下的最佳群速度修正参数的值。对这些数据分析我们发现，三阶迎风紧致群速度控制格式具有最好的性质。而四阶中心紧致群速度控制格式，由于各向异性效应大，不太适合模拟二维问题。四阶中心群速度控制格式也是不错的格式，这个格式将是大量流通矢量不易分裂的方程构造差分格式的首选。数值实验结果表明，利用本文方法构造的差分格式均可以很好地控制非物理振荡并能保证精度，数值结果是令人满意的。

参考文献：

[1]马延文，傅德薰.群速度直接控制四阶迎风紧致格式[J].中国科学A辑，2001,31(6): 554-561.

[2]ZHU Q Y,LI Y.An upwind compact approach with group velocity control for compressible flow fields[J].International Journal for Numerical Methods in Fluids,2004,44(5): 463-482.

[3]朱庆勇.一个求解双曲守恒律方程的高分辨率GVC格式[J].高等学校计算数学学报，2000(2):169-174.

[4]朱庆勇，马延文.求解双曲型守恒律方程的高精度迎风紧致群速度控制法[J].计算物理，1998,15(5):531-536.

[5]汪晓东，朱庆勇.求解带源项双曲守恒律方程的IGVC格式[J].中山大学学报: 自然科学版，2008,11(7):131-135.

[6]朱庆勇，李岳生.数值求解Euler方程的UCGVC差分格式的注记[J].计算数学，2000，22(2)：209-218.

[7]VAN LEER B.Towards the ultimate conservative difference scheme II[J].Journal of Computational Physics,1974(14):361-376.

[8]BERMUDEZ A,DERVIEUX A,DESIDERI J A,et al.Upwind schemes for the two-dimensional shallow water equations with variable depth using unstructured meshes [J].Computer Methods in Applied Mechanics and Engineering,1998,155(1-2): 49-72.

[9]王如云，张东生，张东宽.曲线坐标网格下二维涌波数值模拟的TVD型格式[J].水利学报，2002: 72-77.