浅水动边界问题的位移法模拟

吴锋+钟万勰

摘要：为精确模拟浅水波非线性演化过程中的动边界，提出一种基于位移的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.

（编辑 武晓英）