求解二維淺水波方程的旋轉(zhuǎn)混合格式*

鄭素佩， 李 霄， 趙青宇， 封建湖

（長(zhǎng)安大學(xué) 理學(xué)院，西安 710064）

引 言

隨著計(jì)算機(jī)技術(shù)的不斷發(fā)展，對(duì)于計(jì)算流體動(dòng)力學(xué)中更精確的數(shù)值計(jì)算方法的需求也越來(lái)越迫切.淺水波方程是淺水環(huán)境下流體運(yùn)動(dòng)的一種數(shù)學(xué)描述，是研究河流、灌溉和海洋等波動(dòng)問(wèn)題的一種重要模型，對(duì)其的精確求解還存在很大困難.淺水波方程本質(zhì)上是一類非線性雙曲守恒方程，對(duì)于雙曲守恒律方程數(shù)值求解算法的研究一直是一個(gè)熱點(diǎn)，間斷解的存在對(duì)該類方程的數(shù)值求解格式提出了很高的要求.Tadmor[1]提出了熵穩(wěn)定格式，為得到穩(wěn)定的具有物理意義的解提供了一種簡(jiǎn)單、有效的方法.Roe[2]在自己提出的熵守恒格式的基礎(chǔ)上增加了數(shù)值黏性項(xiàng)，得到的熵穩(wěn)定格式能很好地捕捉激波，抑制間斷處出現(xiàn)的偽振蕩.程曉晗、蔡力等[3]提出了一種基于WENO 重構(gòu)的熵穩(wěn)定格式，該格式具有基本無(wú)振動(dòng)性和高分辨率的特點(diǎn).Liu 等[4]采用SPWENO/WENO-YC3 重構(gòu)構(gòu)造了高分辨率熵穩(wěn)定格式，該格式具有穩(wěn)定、低耗散的特點(diǎn).王令、鄭素佩[5]針對(duì)二維淺水波方程，提出了基于移動(dòng)網(wǎng)格的熵穩(wěn)定格式，有效提高了淺水波方程數(shù)值求解格式的分辨率.

對(duì)雙曲守恒律方程通量函數(shù)的離散方向，Levy 等[6]研究了多維控制方程的旋轉(zhuǎn)Riemann 求解器，討論了通量函數(shù)沿迎風(fēng)方向和界面法線方向的離散數(shù)值求解算法.Ren 提出了一種旋轉(zhuǎn)Roe 格式[7]，該格式不受激波不穩(wěn)定性的影響，從理論上分析了旋轉(zhuǎn)通量函數(shù)的耗散特征.Zhang 等[8]提出了一種旋轉(zhuǎn)迎風(fēng)格式求解二維問(wèn)題，其結(jié)果能夠明顯減少額外的耗散.鄭素佩等[9]基于移動(dòng)網(wǎng)格采用旋轉(zhuǎn)通量求解二維淺水波方程，該算法具有高分辨率的特點(diǎn).為了更好地捕捉激波，同時(shí)保留耗散通量的魯棒性，研究者們采用了混合格式.Nishikawa 和Kitamura[10]為了消除激波不穩(wěn)定性將Roe 格式和HLL 格式結(jié)合起來(lái)，其中，HLL 格式可以抑制激波不穩(wěn)定性，Roe 格式可以避免過(guò)度耗散.劉友瓊等[11]基于旋轉(zhuǎn)Riemann 求解器將HLLC 格式與HLL 格式進(jìn)行特定結(jié)合得到了一類混合型數(shù)值格式，此格式具有結(jié)構(gòu)簡(jiǎn)單、數(shù)值穩(wěn)定、分辨率高等優(yōu)點(diǎn).賈豆、鄭素佩[12]將熵穩(wěn)定格式與HLL 格式結(jié)合提出了一種旋轉(zhuǎn)混合格式，用于求解二維Euler 方程，數(shù)值結(jié)果分辨率高.

本文對(duì)由熵穩(wěn)定格式和HLL 格式結(jié)合得到的旋轉(zhuǎn)混合格式進(jìn)行了研究，介紹了二維淺水波方程、有限體積法與時(shí)間方向離散的Runge-Kutta 法的基礎(chǔ)知識(shí)，給出了旋轉(zhuǎn)通量混合格式的具體表達(dá)式.通過(guò)若干數(shù)值算例驗(yàn)證了該旋轉(zhuǎn)混合格式的性能，表明該混合格式對(duì)于二維淺水波方程數(shù)值求解魯棒性好、分辨率高、具有很好的激波捕捉能力.

1 控制方程

二維淺水波方程[13]為

其中

h表示水深，g表示重力加速度，u,v分 別表示x方 向和y方 向上的水流速度，S表示源項(xiàng).其中S可分為摩擦項(xiàng)和傾斜項(xiàng)：

zb表示底部地勢(shì)函數(shù)，gh(?zb/?x)和gh(?zb/?y) 表示水下底部作用力沿x方 向和y方向上的分力，ghSfx和ghSfy為水下底面摩擦力沿x方 向和y方 向的分量，Sfx和Sfy為x方 向和y方向上的摩阻比率，有

K表示摩擦因數(shù)，通常情況取K=50.

2 方程離散

本文采用任意四邊形單元對(duì)計(jì)算區(qū)域進(jìn)行離散[9]，即控制體為任意四邊形單元.

2.1 旋轉(zhuǎn)不變性[5,14]

對(duì)淺水波方程（1）左右兩邊同時(shí)進(jìn)行積分，得到

其中V表示控制體，A表示V的邊界，H=(F,G)表示通量的張量，n為 邊界面A的單位外法向量， dA表示一個(gè)面積元，H·ndA表示邊界A的通量法向分量.

方程（3）的半離散格式如下：

2.2 數(shù)值通量

2.2.1 熵穩(wěn)定格式

二維淺水波方程熵穩(wěn)定格式可表示為下面的形式[15-19]：

2.2.2 HLL 通量格式

在精確Riemann 求解器的基礎(chǔ)上，Harten 等提出了HLL 通量格式[20-23]，該格式是將兩種波分離成三個(gè)常數(shù)狀態(tài)，具有良好的魯棒性，能有效地消除紅斑現(xiàn)象.HLL 格式的具體表達(dá)式為

2.2.3 混合通量格式

2.3 時(shí)間離散

采用三階強(qiáng)穩(wěn)定Runge-Kutta 法[24]對(duì)時(shí)間進(jìn)行離散，有

3 數(shù)值算例

本節(jié)采用所構(gòu)造的混合格式對(duì)幾個(gè)淺水波方程數(shù)值算例進(jìn)行求解，取網(wǎng)格數(shù)為1 50×150，并對(duì)所得結(jié)果進(jìn)行分析、比較.

例1二維圓柱潰壩問(wèn)題

考慮方程（1）當(dāng)S=0 時(shí)在計(jì)算域[ 0,50]×[0,50]上滿足如下條件的問(wèn)題：

其中 C FL數(shù)為0.45 ， 時(shí)間t=1,g=0.98，采用周期性邊界條件.圖1～3 分別表示二維圓柱潰壩問(wèn)題的模擬圖、密度等值線圖和速度等值線圖.圖2(a)、3(a)和3(c)是僅由熵穩(wěn)定旋轉(zhuǎn)通量格式得到的數(shù)值結(jié)果，圖2(b)、3(b)和3(d)是用本文構(gòu)造的旋轉(zhuǎn)混合通量得到的結(jié)果，可以看出兩種格式均可以捕捉到激波，而旋轉(zhuǎn)混合格式得到的數(shù)值結(jié)果具有更高的分辨率.

圖1 圓柱潰壩模擬圖Fig.1 The simulation diagram for the cylindrical dam

圖2 二維圓柱潰壩問(wèn)題密度等值線圖：(a) 非混合格式密度結(jié)果；(b) 混合格式密度結(jié)果Fig.2 The density contour for the 2D cylindrical dam-break problem：(a) the non-mixed scheme density results; (b) the mixed scheme density results

例2二維區(qū)域上的圓形潰壩問(wèn)題

考慮方程（1）當(dāng)S=0 時(shí)在區(qū)域[ ?1,1]×[?1,1]上滿足如下條件的問(wèn)題：

其中C FL數(shù)為0.25 ， 時(shí)間t=0.2,g=0.98，采用周期性邊界條件.

圖3 二維圓柱潰壩問(wèn)題速度圖：(a) 非混合格式x 方向的速度u；(b) 混合格式x 方向的速度u；(c) 非混合格式y(tǒng) 方向的速度v；(d) 混合格式y(tǒng) 方向的速度vFig.3 The velocity contour for the 2D cylindrical dam-break problem: (a) the non-mixed scheme x direction velocity u results; (b) the mixed scheme x direction velocity u results; (c) the non-mixed scheme y direction velocity v results; (d) the mixed scheme y direction velocity v results

圖4～6 分別表示圓形潰壩問(wèn)題的模擬圖、密度等值線圖和速度等值線圖.圖5(a)、6(a)和6(c)是僅由熵穩(wěn)定旋轉(zhuǎn)通量格式得到的數(shù)值結(jié)果，圖5(b)、6(b)和6(d)是由旋轉(zhuǎn)混合格式得到的結(jié)果，通過(guò)對(duì)比可以看到旋轉(zhuǎn)混合格式得到的數(shù)值結(jié)果具有更高的分辨率，對(duì)稱性好，沒(méi)有出現(xiàn)非物理振蕩.

圖4 圓形潰壩模擬圖Fig.4 The simulation diagram for the circular dam break

圖5 二維圓形潰壩問(wèn)題密度等值線圖：(a) 非混合格式密度結(jié)果； (b) 混合格式密度結(jié)果Fig.5 The density contour for the 2D circular dam-break problem: (a) the non-mixed scheme density results; (b) the mixed scheme density results

圖6 二維圓形潰壩問(wèn)題速度圖：(a) 非混合格式x 方向的速度u；(b) 混合格式x 方向的速度u；(c) 非混合格式y(tǒng) 方向的速度v；(d) 混合格式y(tǒng) 方向的速度vFig.6 The velocity contour for the 2D circular dam-break problem: (a) the non-mixed scheme x direction velocity u results; (b) the mixed scheme x direction velocity u results; (c) the non-mixed scheme y direction velocity v results; (d) the mixed scheme y direction velocity v results

例3二維區(qū)域上的激波聚焦問(wèn)題

考慮方程（1）當(dāng)S=0 時(shí)在區(qū)域[ ?1,1]×[?1,1]上滿足如下條件的問(wèn)題：

其中C FL數(shù)為0.6,時(shí) 間t=0.2,g=0.98,采用周期性邊界條件.圖7～9 分別表示二維激波聚焦問(wèn)題的模擬圖、密度等值線圖和速度等值線圖.圖8(a)、9(a)和9(c)是僅由熵穩(wěn)定旋轉(zhuǎn)通量格式得到的數(shù)值結(jié)果，圖8(b)、9(b)和9(d)是由旋轉(zhuǎn)混合格式得到的結(jié)果，可以看出旋轉(zhuǎn)混合格式過(guò)渡帶更窄且無(wú)振蕩，具有更高的分辨率.

圖7 激波聚焦模擬圖Fig.7 The simulation diagram for the shock wave focusing problem

圖8 二維激波聚焦問(wèn)題密度等值線圖：(a) 非混合格式密度結(jié)果； (b) 混合格式密度結(jié)果Fig.8 The density contour for the 2D shock wave focusing problem: (a) the non-mixed scheme density results; (b) the mixed scheme density results

圖9 二維激波聚焦問(wèn)題速度圖：(a) 非混合格式x 方向的速度u； (b) 混合格式x 方向的速度u；(c) 非混合格式y(tǒng) 方向的速度v ； (d) 混合格式y(tǒng) 方向的速度vFig.9 The velocity contour for the 2D shock wave focusing problem: (a) the non-mixed scheme x direction velocity u results; (b) the mixed scheme x direction velocity u results; (c) the non-mixed scheme y direction velocity v results; (d) the mixed scheme y direction velocity v results

例4二維潮汐模擬問(wèn)題

在區(qū)域[ ?2,2]×[?2,2]上考慮方程(1)，底部地勢(shì)函數(shù)和初始條件為

zb(x,y)=1+0.01cos(πx/2)cos(πy/2),

h(x,y,0)=zb(x,y,0)+ζ(x,y,0),u(x,y,0)=v(x,y,0)=0,

其中初始水深表達(dá)式為

取CFL數(shù)為0.45， 時(shí)間t=0.3,g=0.98, 采用周期性邊界條件.圖10 為二維潮汐問(wèn)題的水面高度模擬圖，圖11為潮汐問(wèn)題的密度等值線圖，圖12 為二維潮汐問(wèn)題的速度圖.圖11(a)、12(a)和12(c)是僅由熵穩(wěn)定旋轉(zhuǎn)通量格式得到的數(shù)值結(jié)果，圖11(b)、12(b)和12(d)是由熵穩(wěn)定和HLL 耦合的旋轉(zhuǎn)混合格式得到的結(jié)果，通過(guò)對(duì)比可以發(fā)現(xiàn)兩種格式?jīng)]有明顯差別， 均對(duì)稱性良好，結(jié)構(gòu)清晰.

圖10 潮汐問(wèn)題模擬圖Fig.10 The simulation diagram for the tidal problem

圖11 二維潮汐問(wèn)題密度等值線圖：(a) 非混合格式密度結(jié)果；(b) 混合格式密度結(jié)果Fig.11 The density contour map for the 2D tidal problem: (a) the non-mixed scheme density results; (b) the mixed scheme density results

圖12 二維潮汐問(wèn)題速度圖：(a) 非混合格式x 方向的速度u；(b) 混合格式x 方向的速度u；(c) 非混合格式y(tǒng) 方向的速度v；(d) 混合格式y(tǒng) 方向的速度vFig.12 The velocity contour for the 2D tidal problem: (a) the non-mixed scheme x direction velocity u results; (b) the mixed scheme x direction velocity v results; (c) the non-mixed scheme y direction velocity v results; (d) the mixed scheme y=direction velocity v results

4 總 結(jié)

本文構(gòu)造了一種求解二維淺水波方程的旋轉(zhuǎn)通量混合格式.利用淺水波方程的旋轉(zhuǎn)不變性將熵穩(wěn)定格式和HLL 格式耦合得到混合格式進(jìn)行數(shù)值模擬時(shí)，對(duì)于帶源項(xiàng)的淺水波方程，混合旋轉(zhuǎn)通量法與非混合格式相比沒(méi)有明顯差別，對(duì)稱性良好，結(jié)構(gòu)清晰；對(duì)于源項(xiàng)為零的淺水波方程，混合旋轉(zhuǎn)通量法具有更高的分辨率.驗(yàn)證可得，新的旋轉(zhuǎn)混合格式魯棒性好、分辨率高、具有很好的激波捕捉能力.