淺水動邊界問題的位移法模擬

吳鋒+鐘萬勰

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

關鍵詞：淺水波； 位移法； 動邊界； 保辛算法； 有限元； Hamilton變分原理； Well-balanced算法； 不平水底

中圖分類號： O353.2

文獻標志碼： 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 引 言

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

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

2 基于位移的淺水方程

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

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

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

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

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

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

4.2 考慮摩阻的拋物線河床

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

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

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

4.4 考慮摩阻潰壩

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

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

5 結束語

本文在Lagrange坐標下研究淺水的動邊界問題，通過Hamilton變分原理導出淺水波方程，數值計算時，可以方便地采用有限元進行空間離散和保辛算法進行時間積分，因此在數值上具有優勢.通過上文的分析過程以及幾個具體算例表明，位移法在分析淺水問題上的優勢可以歸為以下幾點.

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

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

3）對SVE的離散時，如果底坡項的離散不恰當，會導致靜水變動現象的出現，而SWE-D可以保持靜水的靜止狀態，是一種Well-balanced模型.采用SWE-D計算淺水流，不會出現負水深問題，可以精確捕捉動邊界的位置和移動速度.

參考文獻：

[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] 柏祿海， 金生. 帶源項淺水方程的高階格式研究[J]. 水動力學研究與進展： 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]. 水動力學研究與進展： 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]. 大連理工大學學報， 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] 鐘萬勰. 應用力學的辛數學方法[M]. 北京： 高等教育出版社， 2006.

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

[10] 吳鋒， 鐘萬勰. 淺水問題的約束Hamilton變分原理及祖沖之類保辛算法[J]. 應用數學和力學， 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] 鐘萬勰， 高強， 彭海軍. 經典力學-辛講[M]. 大連： 大連理工大學出版社， 2013.

[12] 鐘萬勰， 高強. 約束動力系統的分析結構力學積分[J]. 動力學與控制學報. 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]. 計算力學學報， 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] 邢譽峰， 楊蓉. 動力學平衡方程的Euler中點辛差分求解格式[J]. 力學學報， 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]. 大連： 大連理工大學， 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] 宋利祥， 周建中， 鄒強， 等. 一維淺水方程的強和諧Riemann求解器[J]. 水動力學研究與進展： 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] 郭彥. 基于特征思想的高分辨率格式的研究和應用[D]. 合肥： 中國科學技術大學， 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.

（編輯 武曉英）