周超紅，姚清河，朱慶勇

(中山大學(xué)應(yīng)用力學(xué)與工程系，廣東 廣州 510275)

淺水波方程是計(jì)算流體力學(xué)中的一類重要方程。淺水波方程的應(yīng)用范圍十分廣泛，可用于河道流量、潰壩決堤、污染物輸運(yùn)擴(kuò)散、河道輸沙、洪水預(yù)報(bào)、河口潮汐、涌浪等大量問題。

由于淺水波方程的應(yīng)用背景十分廣泛，其數(shù)值方法的研究自然成為了計(jì)算流休力學(xué)和計(jì)算水動(dòng)力學(xué)的重要課題。由于二維淺水波方程的一般形式比較復(fù)雜，本文中使用在無黏性、無摩擦、河底無傾斜的理想情況下的二維淺水波方程作為求解模型。群速度控制的思想是由傅德薰、馬延文[1]提出的，朱慶勇[2-4]利用這種思想成功構(gòu)造了求解空氣動(dòng)力學(xué)方程的高分辨率格式，之后他們[5]將這一思想引入淺水波方程，構(gòu)造了求解一維淺水波方程的群速度控制格式。

利用群速度控制的思想，本文構(gòu)造了多種求解二維淺水波方程的高分辨率格式。通過引入色散度量以及各向異性度量的概念，為評(píng)估差分算子對微分算子的逼近程度提供了客觀標(biāo)準(zhǔn)。數(shù)值實(shí)驗(yàn)結(jié)果表明，這些格式均可以很好地控制非物理振蕩并保證格式精度。

1 預(yù)備知識(shí)

考慮單波方程及其半離散方程

,f=f(u)

(1)

(2)

這里Fj/Δx是一階導(dǎo)數(shù)?f/?x的差分逼近，定義差分算子如下

±(fj±1-fj)

Fj可以表示為

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

下面我們利用傅里葉分析法對差分方法進(jìn)行理論分析。

考慮f(u)=au，其中a為常數(shù)，初值條件u(x,0)=u0(x)，考慮以2π為周期的單波解，對于半離散方程(2)，我們設(shè)x方向網(wǎng)格步長為

設(shè)其初始條件為uj=exp(ikxj)，F(xiàn)j=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分別決定了差分格式的數(shù)值耗散與數(shù)值色散性質(zhì)。akr>0說明格式是正耗散的，akr<0說明格式是負(fù)耗散的，而akr=0則格式是非耗散的，負(fù)耗散的格式是不穩(wěn)定的。數(shù)值耗散主要影響振幅的大小,而數(shù)值色散效應(yīng)則是產(chǎn)生虛假振蕩、引起波傳播速度不均一的根本原因。高分辨率算法首先要求能夠抑制色散效應(yīng)引起的振蕩，其次也要盡量提高數(shù)值解的精度。

波的群速度，是指波振幅外形上的變化在空間中傳播的速度。根據(jù)以上的討論，我們可定義群速度為sg=dki/dα，其中α=kΔx。對于(1)的精確解(3)，ki=α，sg=1。

2 群速度修正格式的構(gòu)造

在對群速度控制法的研究中我們發(fā)現(xiàn)，文[6]中構(gòu)造三階迎風(fēng)緊致群速度控格式

(5)

其實(shí)就是針對三階迎風(fēng)緊致格式，引入了如下的群速度修正項(xiàng)

對半離散方程(2)進(jìn)行修正。參考上式，我們考慮群速度修正項(xiàng)

其中ssj=sign[(uj+1-uj)(uj+2-uj+1-uj+uj+2)]，稱為激波結(jié)構(gòu)函數(shù)。φj滿足在振蕩區(qū)域φj=1，平滑區(qū)域φj=0，稱為開關(guān)函數(shù)。

可取φj=

0<ω<1,κ=0.05

通過分析Lg我們發(fā)現(xiàn)，在平滑區(qū)域，Lg=0；在激波前后，Lg能使得振蕩波向激波靠攏；在極值點(diǎn)，Lg只有耗散而無色散，有利于振蕩的消除；而在拐點(diǎn)，Lg呈現(xiàn)負(fù)耗散性質(zhì)，則有利于激波形狀的保持。根據(jù)Lg的這些性質(zhì)，我們得知，Lg其實(shí)可以對任意差分格式進(jìn)行修正，而不僅僅限于三階迎風(fēng)緊致格式。

現(xiàn)在我們從算子修正的角度來看待方程(2)，即采用如下半離散格式

(6)

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

我們改稱式(6)中的Lg為群速度修正算子，稱σ為群速度修正系數(shù)。該格式與格式(5)不僅僅是在形式上有所區(qū)別，兩者的意義完全不一樣。前者表示群速度修正算子度基本算子的修正，而后者表示群速度修正項(xiàng)對半離散格式(2)的修正。利用算子修正的觀點(diǎn)來構(gòu)造群速度控制格式，有利于對構(gòu)造的格式進(jìn)行分析。

對于半離散格式(6)，使用不同的基本算子和群速度修正系數(shù)σ，我們就可以構(gòu)造出各種各樣的群速度控制格式。我們將使用格式(6)構(gòu)造群速度控制格式的方法，稱為群速度修正法。

圖1 波數(shù)比例因子與色散顯現(xiàn)因子

圖2 三階迎風(fēng)群速度控制格式的各向異性效應(yīng)

下面以三階迎風(fēng)群速度控制格式為例，利用色散度量函數(shù)，尋找最佳的群速度修正參數(shù)。由于三階迎風(fēng)格式的振蕩主要在激波后，我們考慮，激波后三階迎風(fēng)群速度控制格式的差分算子

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

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

給定不同的σ，將sg代入Dmf右端，對Dmf的右端求數(shù)值積分(可取αn=nπ/10，n=0,1,…,10)，再畫出σ-Dmf的曲線，Dmf最低點(diǎn)對應(yīng)的σ，即為最優(yōu)的修正參數(shù)。由于群速度控制格式的耗散隨著σ的增大而增大，因而若最低點(diǎn)偏左側(cè)點(diǎn)的縱坐標(biāo)與最低點(diǎn)相差不大，最優(yōu)的修正參數(shù)則應(yīng)取左側(cè)點(diǎn)對應(yīng)的σ。

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

對于二維問題，我們還需要考慮算子的各項(xiàng)異性效應(yīng)。考慮二維單波方程

，a,b為常量

(7)

初值u(x,y,0)=exp(iK·X)，其中K=[k1,k2]T，X=[x,y]T，k1,k2分別為x,y方向上的波數(shù)。定義

并將它重寫為l=[cosθ,sinθ]T，這里的θ就是方位角。該問題的精確解可以表示為

(8)

方程(7)的半離散格式為

(9)

其中δxuij/Δx、δyuij/Δy分別為?u/?x、?u/?x的差分逼近。方程(9)的精確解可以表示為

(10)

其中

，，

α=k1Δx，β=k2Δy

對‖l‖的逼近程度。

為度量差分算子的各項(xiàng)異性，仿照前面色散度量的定義，定義

τθdω

為各項(xiàng)異性度量。其中τ= ecos(ω)-e-1稱為波成分比例因子，表示數(shù)值解中ω與對應(yīng)的波的比例。ν=1/σ4/3稱為各向異性顯現(xiàn)因子，用于用于度量耗散對各向異性的壓制力度。

圖4 群速度控制算子的各向異性度量

根據(jù)本節(jié)提供的方法，我們可以得到各種高階群速度控制法的最佳群速度修正參數(shù)，具體為表1。

表1 基本算子最佳群速度修正參數(shù)對應(yīng)表

由于開關(guān)函數(shù)的存在，平直區(qū)域?qū)⒉贿M(jìn)行群速度控制，因而只需考慮主要振蕩區(qū)域?qū)?yīng)的最佳群速度修正參數(shù)的值。對于一維問題只需考慮色散度量的值。例如，三階迎風(fēng)算子的振蕩主要出現(xiàn)在激波后，因而推薦取值為激波后色散度量下的最佳值0.20。對于二維問題，要同時(shí)考慮色散度量和各向異性度量。如四階迎風(fēng)算子的振蕩主要出現(xiàn)在激波前，色散度量最佳值為0.20，各向異性度量下的最佳值為0.10，推薦取兩者的平均值0.15。對于三階迎風(fēng)緊致算子，實(shí)驗(yàn)發(fā)現(xiàn)σ取0.05時(shí)色散比較明顯，故推薦值取為0.10。

3 二維淺水波方程的群速度修正格式

二維淺水波方程在無黏性、無摩擦、河底無傾斜的理想情況下，其數(shù)學(xué)模型的齊次守恒形式如下

(11)

其中

，，

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

為適應(yīng)一般物型計(jì)算。引入隨時(shí)間變化的貼體曲線坐標(biāo)變換

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

則由方程(11)可導(dǎo)出二維曲線坐標(biāo)系下的守恒型淺水波方程

(12)

其中

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

對于雙曲型方程，為了保證計(jì)算我們常常需要使用迎風(fēng)方法，但迎風(fēng)方法要求對流通矢量F*、G*進(jìn)行分裂[7]。文[8]提出了針對二維淺水波方程的有效分裂方法，參照該文中的方法，我們給出曲線坐標(biāo)下二維淺水波方程流通矢量的分裂格式。

下邊構(gòu)造曲面坐標(biāo)下二維淺水波方程的群速度修正格式。記ωi,j=ω(iΔξ,jΔη)，定義算子

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

并令

其中φξ、φη為ξ方向、η方向的開關(guān)函數(shù)，λξ(i,j)、λη(i,j)分別為A*、B*特征值中絕對值最大的特征值。假設(shè)δξ、λξ分別為ξ方向、η方向的基本算子，則曲坐標(biāo)系下二維淺水波方程的群速度修正格式為

其中為群速度修正參數(shù)。若基本算子是迎風(fēng)算子，則上式中

,

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

若基本算子取為四階中心差分算子，δξ、δη滿足

確定基本算子及其對應(yīng)的群速度修正系數(shù)后，我們可以利用三階龍格庫塔法求解半離散方程(13)，進(jìn)而通過U=V/J,求得U。

4 數(shù)值算例

算例1 傾斜水躍問題

一條40 m長的平底河道，上游入口寬度為30 m，來流Froude數(shù)Fr=2.5，在河道一側(cè)10 m處開始以角度α=5°收縮,收縮的邊與超臨界流相互作用形成附體涌波。入口處采用超臨界流固定邊界條件,出口處采用傳輸邊界條件,河道兩邊采用滑移無穿透邊界條件。計(jì)算網(wǎng)格數(shù)為120×90。此算例有精確解，涌波線與X軸的夾角為β=28.32°，波后水深與波前水深比為h2/h1=1.250。

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

當(dāng)流體達(dá)到穩(wěn)定狀態(tài)時(shí),據(jù)圖5 可得到底端激波線與出口相交處點(diǎn)的坐標(biāo)為(40,13.46,1.001),據(jù)此可算得波角β=28.87°,數(shù)值結(jié)果很好地模擬了附體涌波結(jié)構(gòu)。波后水深與波前水深比為h2/h1=1.250～1.251。采用TVD方法[9]得到的結(jié)果為β=29.01°，h2/h1=1.246～1.254，顯然本文使用的四階群速度控制法結(jié)果比較準(zhǔn)確。

圖5 四階群速度控制法

算例2 二維圓形潰壩問題

一個(gè)50 m×50 m方形區(qū)域，區(qū)域被一圓心在區(qū)域中心(25,25)，半徑r=11 m的圓柱面形壩分為兩部分，壩內(nèi)水深10 m，壩外水深1 m。假定壩體突然潰破，計(jì)算網(wǎng)格數(shù)取50×50。

本算例算法使用三階迎風(fēng)群速度控制法σ=0.15。圖6和圖7給出了t=0.7 s時(shí)水波運(yùn)動(dòng)的狀態(tài)。二維圓形潰壩問題的真實(shí)解應(yīng)當(dāng)滿足在任意時(shí)刻水深等值線是一系列以(25，25)為中心的圓。圖7(b)中所示的水深等值線與標(biāo)準(zhǔn)圓相當(dāng)接近，說明三階迎風(fēng)群速度控制格式能夠有效地對各個(gè)方向激波進(jìn)行群速度控制，各向異性效應(yīng)極小，能很好地模擬二維淺水波的運(yùn)動(dòng)。

圖6 水深曲面

圖7 三階迎風(fēng)群速度控制法

5 結(jié)論

本文得到了幾種常用高階差分算子在色散度量和各向異性度量下的最佳群速度修正參數(shù)的值。對這些數(shù)據(jù)分析我們發(fā)現(xiàn)，三階迎風(fēng)緊致群速度控制格式具有最好的性質(zhì)。而四階中心緊致群速度控制格式，由于各向異性效應(yīng)大，不太適合模擬二維問題。四階中心群速度控制格式也是不錯(cuò)的格式，這個(gè)格式將是大量流通矢量不易分裂的方程構(gòu)造差分格式的首選。數(shù)值實(shí)驗(yàn)結(jié)果表明，利用本文方法構(gòu)造的差分格式均可以很好地控制非物理振蕩并能保證精度，數(shù)值結(jié)果是令人滿意的。

參考文獻(xiàn)：

[1]馬延文，傅德薰.群速度直接控制四階迎風(fēng)緊致格式[J].中國科學(xué)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]朱慶勇.一個(gè)求解雙曲守恒律方程的高分辨率GVC格式[J].高等學(xué)校計(jì)算數(shù)學(xué)學(xué)報(bào)，2000(2):169-174.

[4]朱慶勇，馬延文.求解雙曲型守恒律方程的高精度迎風(fēng)緊致群速度控制法[J].計(jì)算物理，1998,15(5):531-536.

[5]汪曉東，朱慶勇.求解帶源項(xiàng)雙曲守恒律方程的IGVC格式[J].中山大學(xué)學(xué)報(bào): 自然科學(xué)版，2008,11(7):131-135.

[6]朱慶勇，李岳生.數(shù)值求解Euler方程的UCGVC差分格式的注記[J].計(jì)算數(shù)學(xué)，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]王如云，張東生，張東寬.曲線坐標(biāo)網(wǎng)格下二維涌波數(shù)值模擬的TVD型格式[J].水利學(xué)報(bào)，2002: 72-77.