淺水動(dòng)邊界問題的位移法模擬

吳鋒+鐘萬勰

摘要：為精確模擬淺水波非線性演化過程中的動(dòng)邊界，提出一種基于位移的Hamilton變分原理，并進(jìn)而導(dǎo)出一種基于位移的淺水方程（Shallow Water Equation based on Displacement，SWE-D）.SWE-D以位移為基本未知量，可以精確滿足動(dòng)邊界處的零水深要求并精確捕捉動(dòng)態(tài)邊界位置，且解具有協(xié)調(diào)性.在Hamilton變分原理的框架下，分別采用有限元和保辛積分算法對(duì)該淺水方程進(jìn)行空間離散和時(shí)間積分，可有效地處理不平水底情況，保證對(duì)非線性演化進(jìn)行長(zhǎng)時(shí)間仿真的精度.數(shù)值算例表明該方法適用于淺水動(dòng)邊界問題的數(shù)值模擬.

關(guān)鍵詞：淺水波； 位移法； 動(dòng)邊界； 保辛算法； 有限元； Hamilton變分原理； Well-balanced算法； 不平水底

中圖分類號(hào)： O353.2

文獻(xiàn)標(biāo)志碼： A

Abstract：To exactly simulate the moving boundaries of shallow water flow in the procedure of nonlinear evolution， a Hamilton variational principle based on displacement is proposed. Furthermore， a Shallow Water Equation based on Displacement （SWE-D） is developed in terms of the Hamilton variational principle. Taking the displacement as a basic unknown variable， SWE-D can exactly satisfy the requirement of zero water depth at moving boundary and exactly capture the location of the moving boundary， and the solutions are well-balanced. In the frame of Hamilton variational principle， the finite element method and symplectic method are respectively used for the spatical discretization and time integral， which can effectively deal with the uneven water bottom and keep the accuracy in simulating the long time nonlinear evolution. The numerical examples show that the method is suitable to the simulation on the shallow water flow with moving boundaries.

Key words：shallow water wave； displacement method； moving boundary； symplectic algorithm； finite element； Hamilton variational principle； Well-balanced algorithm； uneven water bottom

0 引 言

淺水模型在河道流動(dòng)、近海波浪等實(shí)際工程中有十分重要的應(yīng)用.目前最常使用的淺水方程是圣維南方程（de Saint-Venant system of Equations，SVE）[1]，該方程是在Euler坐標(biāo)下建立的基于流速的非線性一階偏微分方程，其中主要由對(duì)流項(xiàng)和源項(xiàng)組成.對(duì)流項(xiàng)的存在使得其空間離散比較講究，目前常用的方法是有限體積法[2]，原因在于有限體積法的離散格式具有物理意義，可以滿足質(zhì)量守恒和動(dòng)量守恒.當(dāng)考慮水底不平順和水底存在摩阻時(shí)，SVE中會(huì)出現(xiàn)源項(xiàng).源項(xiàng)一般分為兩項(xiàng)：一項(xiàng)是底坡項(xiàng)，由水底不平順帶來的；另一項(xiàng)為摩阻項(xiàng)，可以視為淺水流動(dòng)的阻尼.但是，當(dāng)考慮水底不平順?biāo)鶐淼牡灼略错?xiàng)時(shí)，方程的守恒性被破壞[3].底坡項(xiàng)的不適當(dāng)離散常常導(dǎo)致靜水變動(dòng)現(xiàn)象，即所謂解的協(xié)調(diào)性問題，而此時(shí)質(zhì)量守恒也不能保證.另一個(gè)會(huì)破壞質(zhì)量守恒的因素是干濕交界面的處理.在干河床或海灘上，由于水的流動(dòng)，導(dǎo)致干、濕2個(gè)區(qū)間之間的交界面隨時(shí)間不斷變化，也就是所謂的動(dòng)邊界問題.目前基于流速的SVE的數(shù)值方法在處理動(dòng)邊界問題中存在困難[4]，主要有：1）動(dòng)邊界的判定問題，在Euler坐標(biāo)下，以流速為基本變量，計(jì)算區(qū)域必須包含干區(qū)域和濕區(qū)域，而往往計(jì)算網(wǎng)格是固定的，因此就必須判定干濕交界面的位置和判定單元是干還是濕，對(duì)干、濕2種狀態(tài)的判定常常導(dǎo)致計(jì)算數(shù)值的震蕩；2）由于在干濕交界面上的理論水深是0，而數(shù)值計(jì)算時(shí)總會(huì)帶有誤差，因此會(huì)出現(xiàn)負(fù)水深問題；3）當(dāng)考慮摩阻項(xiàng)時(shí)，水深會(huì)出現(xiàn)在摩阻項(xiàng)的分母中，而在干濕交界面處的水深為0，因此干濕交界處是摩阻項(xiàng)的一個(gè)奇點(diǎn)，這也導(dǎo)致數(shù)值計(jì)算上的困難，若出現(xiàn)負(fù)水深，則還可能導(dǎo)致負(fù)阻尼的出現(xiàn).

DURAN等[5]指出，用于非線性淺水問題的有效算法，必須滿足以下3個(gè)條件：1）能很好地處理動(dòng)邊界；2）針對(duì)不同拓?fù)湫螤畹牡灼马?xiàng)，必須具有方便且統(tǒng)一的離散格式；3）能保持靜水的靜止?fàn)顟B(tài).這3個(gè)要求是針對(duì)SVE而言的.然而，從基于流速的SVE本身來看，要同時(shí)滿足這3個(gè)要求是很困難的，根本原因在于SVE是在Euler坐標(biāo)下建立的，其基本未知量是流速，而動(dòng)邊界問題以及水深實(shí)際上是位置問題，都需要用位移來描述，因此如果從位移法入手，問題就好辦得多.近年來，對(duì)于淺水問題的位移法研究逐漸增多.文獻(xiàn)[6]在Lagrange坐標(biāo)下研究淺水波，以位移為基本變量，并假定水平位移與豎向坐標(biāo)z無關(guān)，給出基于位移的淺水波方程，并將該方程用于研究三峽升船機(jī)水箱中水的晃動(dòng)問題.文獻(xiàn)[7]和[8]中均采用位移法研究淺水波中孤立波問題[9]，得到孤立波的一個(gè)解析解，并與傳統(tǒng)Euler坐標(biāo)下的孤立波解進(jìn)行數(shù)值比較.文獻(xiàn)[10]基于Lagrange坐標(biāo)，以位移和壓強(qiáng)為基本變量，將不可壓縮條件視為約束，提出一種微分-代數(shù)方程形式的淺水波方程以及相應(yīng)的約束Hamilton變分原理，并利用保辛的祖沖之類算法[11-13]研究該淺水波方程的求解.已有的研究表明位移法十分適用于淺水問題的數(shù)值仿真.本文將位移法用于淺水動(dòng)邊界問題的仿真計(jì)算，通過數(shù)值算例發(fā)現(xiàn)位移法十分適用于淺水動(dòng)邊界問題的研究.本文方法可以精確計(jì)算動(dòng)邊界處的零水深，可以精確捕捉動(dòng)態(tài)邊界，且沒有靜水變動(dòng)問題.一般將能保持靜水靜止?fàn)顟B(tài)的算法稱為Well-balanced算法[5]，因此本文方法是一種已經(jīng)實(shí)現(xiàn)Well-balanced的算法.數(shù)值算例表明本文方法適用于淺水動(dòng)邊界問題的數(shù)值模擬.

2 基于位移的淺水方程

2.1 無摩阻項(xiàng)的位移淺水方程（Shallow Water Equation based on Displacement， SWE-D）

考慮圖1所示淺水域，定義u（x，z，t）為初始時(shí)刻在坐標(biāo)（x，z）處的質(zhì)點(diǎn)在t時(shí)刻的水平位移；w（x，z，t）為初始時(shí)刻在坐標(biāo)（x，z）處的質(zhì)點(diǎn)在t時(shí)刻的豎向位移.水底形狀為z=-h（x），其中0≤x≤L.初始時(shí)刻的水面z=α（x）.

當(dāng)不考慮摩阻項(xiàng)時(shí)，淺水波是保守系統(tǒng)，其時(shí)間積分最好采用保辛算法.當(dāng)考慮摩阻項(xiàng)時(shí)，摩阻項(xiàng)相當(dāng)于動(dòng)力學(xué)中的阻尼.有文獻(xiàn)研究表明，對(duì)于帶有阻尼的動(dòng)力非保守系統(tǒng)，保辛算法的性能仍然比傳統(tǒng)非保辛算法要好[14]，因此無論考不考慮摩阻項(xiàng)，都建議采用保辛算法.關(guān)于保辛算法，有許多文獻(xiàn)（如文獻(xiàn)[11]）可以參考，這里不再闡述.

需要注意的是FREI[17]是在Euler坐標(biāo)下基于SVE對(duì)此問題展開研究的，沒有考慮豎向加速度的影響.現(xiàn)在采用位移法對(duì)該問題進(jìn)行仿真計(jì)算，也不考慮豎向加速的影響，將仿真結(jié)果與解析解進(jìn)行比較.1和2 s時(shí)的水滴輪廓以及速度分布分別見圖2和3.

由圖2和3可以看出，本文的位移法解與解析解吻合很好.圖3中亦給出高分辨率算法結(jié)果，數(shù)據(jù)是通過GetDataW_cn.exe軟件從文獻(xiàn)[16]抓取的.文獻(xiàn)[16]基于SVE采用高分辨率算法研究此問題，其計(jì)算結(jié)果流速的計(jì)算誤差較大，尤其是在干濕交界（動(dòng)邊界）處誤差更加明顯.這說明本文方法在處理動(dòng)邊界問題上具有顯著優(yōu)勢(shì).

由圖4和5可知：當(dāng)考慮豎向加速度時(shí)，水滴坍塌更慢，水滴形狀也不再能保持拋物線形狀，而速度場(chǎng)的分布也不再是線性的.考慮豎向加速度與不考慮豎向加速度兩者計(jì)算結(jié)果差異較大，顯然考慮豎向加速度更合乎實(shí)際的物理情況.

4.2 考慮摩阻的拋物線河床

SAMPSON等[18]曾解析地分析過拋物線型河床的淺水問題，可用于檢驗(yàn)本文方法處理干濕界面和不平水底問題的能力.初始時(shí)刻的水域剖面見圖6.

圖7給出不同時(shí)刻數(shù)值水面與解析解的比較，其中，x=-a處的位移u（-a，t）隨時(shí)間的周期變化見圖8.由圖7和8可知，本文的數(shù)值解與解析解吻合，表明本文提出的SWE-D可準(zhǔn)確處理底坡項(xiàng)和干濕界面問題[19]，同時(shí)也表明保辛算法在長(zhǎng)時(shí)間仿真帶有阻尼的淺水波動(dòng)問題中具有優(yōu)勢(shì).

由圖10和11可以看出本文方法的計(jì)算結(jié)果與解析解吻合很好.圖11同時(shí)給出文獻(xiàn)[2]和[21]計(jì)算該問題得到的速度分布，其中速度的數(shù)據(jù)是通過GetDataW_cn.exe軟件從文獻(xiàn)[2]和[21]中抓取的.文獻(xiàn)[2]和[21]中采用的模型是SVE，以流速為基本未知量，采用基于特征思想的高分辨率格式數(shù)值算法，計(jì)算得到的數(shù)值水面與解析水面吻合很好，這里不再給出.其數(shù)值流速與解析流速偏差較大，主要體現(xiàn)在潰壩后涌出的水與下游干河面的交界處（也就是動(dòng)邊界處）的流速偏差較大.實(shí)際上動(dòng)邊界問題是Euler坐標(biāo)下天然存在的困難，有許多其他文獻(xiàn)針對(duì)這一問題展開研究，然而都是采用SVE，因此在計(jì)算本算例時(shí)效果不佳，如文獻(xiàn)[22]中的守恒格式和文獻(xiàn)[23]中的CWENO-type中心逆風(fēng)格式.

4.4 考慮摩阻潰壩

SCHOKLITSH[24]曾經(jīng)做過一個(gè)潰壩試驗(yàn).試驗(yàn)中水槽寬0.096 m，高0.080 m，長(zhǎng)20.000 m.水槽由光滑木材制成，在10.000 m處有一壩，壩的左邊蓄水，蓄水深0.074 m，壩右邊水槽不蓄水以模擬干河床，其Manning糙率為0.009 s/m1/3.[25]壩在瞬間被抽走，水開始向右涌出.在文獻(xiàn)[25]和[26]中也可以查到該試驗(yàn)的數(shù)據(jù).這里采用位移法對(duì)該問題進(jìn)行仿真計(jì)算，在[0， 10] m內(nèi)采用200個(gè)線性單元空間離散，采用Euler中點(diǎn)辛差分格式進(jìn)行時(shí)間積分，時(shí)間步長(zhǎng)取0.01 s，將仿真數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)進(jìn)行比較.潰壩后3.75和9.40 s的實(shí)測(cè)水面與仿真水面見圖12，其中實(shí)測(cè)水面是通過GetDataW_cn.exe軟件從文獻(xiàn)[25]中抓取的.

由圖12可知：實(shí)測(cè)的水面輪廓與仿真的水面輪廓吻合很好，這說明本文方法適用于分析帶有摩阻的干底河床的潰壩問題.

5 結(jié)束語

本文在Lagrange坐標(biāo)下研究淺水的動(dòng)邊界問題，通過Hamilton變分原理導(dǎo)出淺水波方程，數(shù)值計(jì)算時(shí)，可以方便地采用有限元進(jìn)行空間離散和保辛算法進(jìn)行時(shí)間積分，因此在數(shù)值上具有優(yōu)勢(shì).通過上文的分析過程以及幾個(gè)具體算例表明，位移法在分析淺水問題上的優(yōu)勢(shì)可以歸為以下幾點(diǎn).

1）SWE-D可以通過對(duì)作用量進(jìn)行變分得到.作用量中包含動(dòng)能和勢(shì)能，均具有鮮明的物理意義.通過空間有限元離散和變分原理所建立的離散格式，其質(zhì)量矩陣和剛度矩陣均為對(duì)稱矩陣，具有保辛的性質(zhì).在SVE的有限元離散中，常常通過加權(quán)殘數(shù)法離散，離散過程不具備物理意義.

2）SVE是一階偏微分方程組，以水平流速和水深為基本變量，而SWE-D是二階偏微分方程，僅僅以水平位移為基本變量，且不存在對(duì)流項(xiàng)，因此空間離散方便，離散后的自由度較少.空間離散后，可以采用保辛算法進(jìn)行時(shí)間積分，因此具有長(zhǎng)時(shí)間仿真的優(yōu)勢(shì).

3）對(duì)SVE的離散時(shí)，如果底坡項(xiàng)的離散不恰當(dāng)，會(huì)導(dǎo)致靜水變動(dòng)現(xiàn)象的出現(xiàn)，而SWE-D可以保持靜水的靜止?fàn)顟B(tài)，是一種Well-balanced模型.采用SWE-D計(jì)算淺水流，不會(huì)出現(xiàn)負(fù)水深問題，可以精確捕捉動(dòng)邊界的位置和移動(dòng)速度.

參考文獻(xiàn)：

[1] CHALFEN M， NIEMIEC A. Analytical and numerical-solution of saint-venant equations[J]. Journal of Hydrology， 1986， 86（1-2）： 1-13. DOI： 10.1016/0022-1694（86）90002-8.

[2] GUO Y， LIU R， DUAN Y， et al. A characteristic-based finite volume scheme for shallow water equations[J]. Journal of Hydrodynamics， 2009， 21（4）： 531-540. DOI： 10.1016/S1001-6058（08）60181-X.

[3] 柏祿海， 金生. 帶源項(xiàng)淺水方程的高階格式研究[J]. 水動(dòng)力學(xué)研究與進(jìn)展： A輯. 2009， 24（1）： 22-28.

BAI L H， JIN S. Study on high resolutio scheme for shallow water equation with source terms[J]. Journal of Hydrodinamics： Ser. A， 2009， 24（1）： 22-28.

[4] BRUFAU P， GARCIA-NAVARRO P， VZQUEZ-CENDN M E. Zero mass error using unsteady wetting-drying conditions in shallow flows over dry irregular topography[J]. International Journal for Numerical Method in Fluids， 2004， 45： 1047-1082. DOI： 10.1002/fld.729.

[5] DURAN A， LIANG Q， MARCHE F. On the well-balanced numerical discretization of shallow water equations on unstructured meshes[J]. Journal of Computational Physics， 2013（235）： 565-586. DOI： 10.1016/j.jcp.2012.10.033.

[6] 鐘萬勰， 陳曉輝. 淺水波的位移法求解[J]. 水動(dòng)力學(xué)研究與進(jìn)展： A輯， 2006， 21（4）： 486-493.

ZHONG W X， CHEN X H. Solving shallow water waves with the displacement method[J]. Journal of Hydrodynamics： Ser. A， 2006， 21（4）： 486-493.

[7] 鐘萬勰， 姚征. 位移法淺水孤立波[J]. 大連理工大學(xué)學(xué)報(bào)， 2006， 46（1）： 151-156. DOI： 10.3321/j.issn：1000-8608.2006.01.028.

ZHONG W X， YAO Z. Shallow water solitary waves based on displacement method[J]. Journal of Dalian University of Technology， 2006， 46（1）： 151-156. DOI： 10.3321/j.issn：1000-8608.2006.01.028.

[8] 鐘萬勰. 應(yīng)用力學(xué)的辛數(shù)學(xué)方法[M]. 北京： 高等教育出版社， 2006.

[9] REMOISSENET M. Waves called solitons： Concepts and experiments[M]. Berline： Springer， 1996.

[10] 吳鋒， 鐘萬勰. 淺水問題的約束Hamilton變分原理及祖沖之類保辛算法[J]. 應(yīng)用數(shù)學(xué)和力學(xué)， 2016， 37（1）： 3-15. DOI： 10.3879/j.issn.1000-0887.2016.01.001.

WU F， ZHONG W X. The constrained Hamilton variational principle for shallow water problems and the Zu-type symplectic algorithm[J]. Applied Mathematics and Mechanics， 2016， 37（1）： 3-15. DOI： 10.3879/j.issn.1000-0887.2016.01.001.

[11] 鐘萬勰， 高強(qiáng)， 彭海軍. 經(jīng)典力學(xué)-辛講[M]. 大連： 大連理工大學(xué)出版社， 2013.

[12] 鐘萬勰， 高強(qiáng). 約束動(dòng)力系統(tǒng)的分析結(jié)構(gòu)力學(xué)積分[J]. 動(dòng)力學(xué)與控制學(xué)報(bào). 2006， 4（3）： 193-200.

ZHONG W X， GAO Q. Integration of constrained dynamical system via analytical structural mechanics[J]. Journal of Dynamics and Control， 2006， 4（3）： 193-200.

[13] 吳鋒， 鐘萬勰. 基于祖沖之類方法具有保辛性[J]. 計(jì)算力學(xué)學(xué)報(bào)， 2015， 32（4）： 447-450.

WU F， ZHONG W X. The Zu-type method is symplectic[J]. Chinese Journal of Computational Mechanics， 2015， 32（4）： 447-450.

[14] 邢譽(yù)峰， 楊蓉. 動(dòng)力學(xué)平衡方程的Euler中點(diǎn)辛差分求解格式[J]. 力學(xué)學(xué)報(bào)， 2007， 39（1）： 100-105.

XING Y F， YANG R. Application of Euler midpoint symplectic integration method for the solution of dynamic equilibrium equations[J]. Chinese Journal of Theoretical and Applied Mechanics， 2007， 39（1）： 100-105.

[15] SCH R C， SMOLARKIEWICZ P K. A synchronous and iterative Flux-Correction formalism for coupled transport equations[J]. Journal of Computational Physics， 1996， 128（1）： 101-120. DOI： 10.1006/jcph.1996.0198.

[16] 柏祿海. 淺水方程高分辨率算法的研究[D]. 大連： 大連理工大學(xué)， 2013.

[17] FREI C. Dynamics of a two-dimensional ribbon of shallow water on a f-plane[J]. Tellus， 1993， 45A（1）： 44-53.

[18] SAMPSON J， EASTON A， SINGH M. Moving boundary shallow water flow above parabolic bottom topography[J]. ANZIAM Journal， 2005（47）： 373-387.

[19] 宋利祥， 周建中， 鄒強(qiáng)， 等. 一維淺水方程的強(qiáng)和諧Riemann求解器[J]. 水動(dòng)力學(xué)研究與進(jìn)展： A輯， 2010， 25（2）： 231-238. DOI： 10.3969/j.issn.1000-4874.2010.02.013.

SONG L X， ZHOU J Z， ZOU Q， et al. A well-balanced Riemann solver for one-dimensional shallow water equations[J]. Journal of Hydrodynamics： Ser. A， 2010， 25（2）： 231-238. DOI： 10.3969/j.issn.1000-4874.2010.02.013.

[20] IZEM N， SEAID M， WAKRIM M. A discontinuous Galerkin method for two-layer shallow water equations[J]. Mathematics and Computers in Simulation， 2016， 120： 12-23. DOI： 10.1016/j.matcom.2015.04.009.

[21] 郭彥. 基于特征思想的高分辨率格式的研究和應(yīng)用[D]. 合肥： 中國科學(xué)技術(shù)大學(xué)， 2009.

[22] MOHAMMADIAN A， le ROUX D Y， TAJRISHI M. A conservative extension of the method of characteristics for 1-D shallow flows[J]. Applied Mathematical Modelling， 2007， 31（2）： 332-348. DOI： 10.1016/j.apm.2005.11.018.

[23] FENG J， CAI L， XIE W. CWENO-type central-upwind schemes for multidimensional Saint-Venant system of shallow water equations[J]. Applied Numerical Mathematics， 2006， 56（7）： 1001-1017. DOI： 10.1016/j.apnum.2005.09.002.

[24] SCHOKLITSCH A. Ueber dammbruchwellen[J]. Sitzungsberichte Der Kaiserlichen Akademie Wissenschaften， Viennal， 1917， 126： 1489-1514.

[25] KHAN A A， LAI W. Modeling shallow water flows using the discontinuous Galerkin method[M]. Boca Raton： CRC Press， 2014.

[26] KHAN A A. Modeling flow over an initially dry bed[J]. Journal of Hydraulic Research， 2000， 38（5）： 383-388.

（編輯 武曉英）