基于无积分节点间断有限元的二维水沙模型：(1)水动力

李文俊，张庆河，李龙翔，冉国全

(天津大学 水利工程仿真与安全国家重点实验室，天津 300072)

数值模拟方法在二维河流与海岸自由水面流动模拟中得到了广泛应用，它可以合理预测水动力变化规律，并为模拟二维泥沙运动和地形演变、预测污染物输移、指导工程建设[1]等工作提供基础。过去几十年间，二维水动力模型取得了相当多的成果，有限差分法、有限体积法和有限元法等都在求解二维浅水方程中获得了大量应用。最近若干年来，间断有限元方法也逐渐开始应用于二维浅水方程的数值模拟，一系列研究者的工作使得这一模型逐步具有了计算干湿、保证和谐特性、保证水深为正等特性[2-4]。这一模型也逐渐被应用到立波、溃坝以及大范围潮流、海啸、风暴潮等一系列复杂水动力模拟中[5-8]。

间断有限元方法融合了有限元法和有限体积法的诸多优点，不仅可以通过提高单元插值基函数的阶数来获得高精度数值结果，而且可以对数值通量采用有限体积法中的Riemann算子精确计算，确保间断捕捉的准确性和质量的守恒性。模型可采用非结构化网格，能够很好地模拟复杂岸线，同时可应用多种自适应网格算法。另外，在求解给定单元内部的变量时，间断有限元只需要相邻单元的变量，使得处理器之间的信息传递量保持最小，易于实现算法的并行化[9]。当然，间断有限元方法也存在未知数个数较多、计算量较大的缺点。针对这一问题，李龙翔和张庆河[10]建立了任意四边形网格的无积分节点间断有限元模型。本模型相比于已有间断有限元浅水模型，将数值积分转换为矩阵运算，省略传统数值计算中求解积分节点处函数值的过程，减少计算量的同时减少了单元的内存储量大小。同时本模型可采用任意四边形网格计算，相比于三角形网格可以有效地减少未知量个数，从而减少计算量。

为了开发高精度高效间断有限元二维水沙模型，本文将在李龙翔和张庆河[10]基于任意四边形网格无积分节点间断有限元方法求解对流扩散方程的模型基础上，将控制方程扩展为二维浅水方程，并加入相关源项，建立间断有限元二维水动力模型。该模型将考虑底坡源项、科氏力、风应力、底摩阻和地形干湿变化，以充分反映在各种实际条件影响下浅水流动的二维水动力特性。

1 二维水动力数学模型

1.1 控制方程

非恒定二维浅水方程的守恒形式为

(1)

式中：Q=[h,hu,hv]T为守恒变量，其中h为水深，u、v分别为沿x和y方向的水平流速。通量项E和G分别为x，y方向的通量，表达式为

源项S考虑了底坡、科氏力、底摩阻和表面风应力影响，最终表达式为

(2)

式中：η=h+z为水位，z为底坡高程，ρ为水体密度。

科氏力系数f可由式(3)给出[11]

f=2ωsinφ

(3)

式中：ω为地球绕地轴旋转的角速度，大小为ω=7.29×10-5rad/s，φ为地理纬度。

τbx和τby分别为水底床面沿x和y方向的摩阻应力，τwx和τwy分别为自由表面沿x和y方向的风场切应力,分别由经验公式(4)和(5)给出[12]

(4)

(5)

式中：n为糙率，Cw为无因次风应力系数，ρa为空气密度，ωw为水体自由表面以上10 m处的风速，β为风向角，一般取风向与x方向的夹角。

上述方程写成以分量形式表示的方程即

(6)

为了推导方便，下面的数值离散仍以矢量形式为主。

1.2 节点间断有限元数值离散

利用在单元上连续可微、单元间允许存在跳跃的分段光滑函数Qh逼近Q[13]，Qh表达式为

(7)

在节点间断有限元方法中，插值基函数与试验函数ji(x)均采用Lagrange插值函数定义，定义分段光滑函数Qh所在的空间为Vh。利用Qh替换控制方程中的Q，并用试验函数v∈Vh乘以替换后的公式，在单元上积分，可得

(8)

对上式用分部积分和格林公式可得

(9)

式中：F表示通过边界的法向通量，n为边界∂Ωk处的单位外法向量，式(9)便是弱形式的空间离散方程组。由于间断有限元方法中物理通量F在边界两侧的值不相同，故将通过边界的法向通量用数值通量F*代替[14]。利用二维浅水方程的旋转不变性，将整体坐标系绕原点逆时针旋转角度φ(φ为边界外法线向量与x轴正方向的夹角)，使坐标轴x向与边界外法线向量平行，并将原点平移至对应边界中点处，从而建立局部坐标系。由此将数值通量的求解转化为局部坐标系下关于守恒变量向量Q的一维Riemann问题的求解。本文采用HLL格式[15]的数值通量，其形式为

(10)

其中，QL和QR、EL和ER分别为局部坐标系下交界面左右的守恒变量、通量向量值

(11)

SL和SR为波速，分别由式(12)和式(13)计算得到[16]

(12)

(13)

其中，

(14)

(15)

(16)

对式(9)左端第二、第三项再运用一次分部积分，得到方程组

(17)

式(17)为强形式的空间离散方程组。相对于弱形式而言，强形式无积分离散格式的存储量更小。

本文采用4阶5步显式Runge-Kutta(RK)格式对半离散方程(17)进行时间积分，最终可得到不同时刻具有高阶时间与空间精度的数值解。

1.3 无积分格式

传统间断有限元方法在求解式(17)时采用数值积分，这种方法要求被积函数在积分节点处的值已知，因此需要先计算出积分节点处的守恒变量值与通量值，从而造成计算量较大而影响计算效率。为了提高计算效率，Atkins和Shu[17]曾提出一种无积分格式，通过假定单元内雅克比行列式为常数，将数值积分运算转化为标准单元内矩阵计算。这一方法省去了数值积分过程，但其只适用于结构化平行四边形网格，对于任意四边形单元并不适用。李龙翔和张庆河[10]提出的无积分格式无需假定单元内雅克比行列式为常数，可适用于任意四边形网格。该格式将积分节点和插值节点均选定为同一套边界点的Legendre-Gauss-Lobatto(LGL)节点，在计算式(17)时将任意四边形单元通过雅克比变换映射到标准四边形单元，并将体积分中的单元刚度矩阵表示为质量矩阵与导数矩阵乘积形式，最终可得到标准单元坐标系下包含NP个方程的方程组(18)

(18)

式中：NP为单元插值节点个数，1≤j≤NP，r、s为标准单元坐标变量，J为雅克比变换矩阵行列式，由式(19)给出

(19)

无积分格式将数值积分转化为上述矩阵的运算，避免了求解积分节点处的函数值。同时，确定标准单元后，其导数矩阵、体积分质量矩阵和面积分质量矩阵都为定值，可以预先调用矩阵运算库计算，从而有效减少计算量并节省大量储存空间。特别需要指出的是，当采用高阶精度格式时，无积分运算求解数值积分产生的误差与采用高阶基函数近似产生的截断误差同阶[18]。因此，无积分方法同样对于格式计算误差影响不大。

1.4 干湿处理及和谐特性

1.4.1 干湿处理

近岸二维水动力模型模拟实际问题的难点在于处理复杂的干湿边界条件[19]。对于水位变化引起的岸线变化的处理，最主要的两种方法是动边界法和干湿算法。其中动边界往往需要持续地对网格进行重新划分，所需计算量偏大。因此，二维水动力模型大多采用干湿算法。干湿处理需要解决三个主要问题：(1)保持水深为正值；(2)消除非物理压力梯度造成的虚假振荡；(3)保持干湿交界处的数值稳定[20]。

1.4.2 和谐性

在守恒型浅水方程模型中，水位梯度项被分为压力项和底坡源项之和

(20)

为了使控制方程离散后该等式始终成立，即满足静水条件下水位为常数、流速为0的和谐(Well-balanced)特性[24]，需要引入一些特殊处理手段，如：改变单元水力模型，消去底坡项；或者对底坡源项进行改造，使底坡源项与压力项在静水条件下能够平衡[25]。本模型采用改造底坡源项的方法，将底坡源项作如下拆分

(21)

并采用Bollermann[26]的方法，将底坡方程用一个连续分段线性函数近似，以保证模型的和谐特性。

2 模型验证

按照上述数值模式建立了二维水动力模型，下面将分别利用不同算例对模型进行验证，其中包括：(1)利用赤道孤立Rossby波验证模型描述浅水长波和科氏力的能力；(2)利用平底坡Stommel问题验证模型描述风应力作用以及科氏力与底摩阻联合作用的能力；(3)利用平缓岸滩波浪爬坡算例验证模型的干湿处理能力；(4)利用斜坡静水算例验证模型的和谐特性。

2.1 赤道孤立Rossby波

本算例模拟的是赤道处孤立Rossby波的运动情况，算例的计算域取为[0,48]×[0,16]的矩形区域，初始的孤立Rossby波位于计算域中心，即坐标[24,8]处，其0阶初始波形由式(22)给出[27]

(22)

式中：Γ(x)=0.771a2sech2(ax)，a为确定孤立波振幅的参数，取为0.395。

算例求解的是无量纲化的二维浅水方程，方程中源项除了底坡源项仅包括科氏力项。由于弱非线性和弱色散之间的平衡，该问题的摄动解预测赤道孤立Rossby波会保持其形态向西行进[28]。采用计算网格为4 181个任意四边形单元。

1-a t=0 s1-b t=10 s

1-c t=20 s1-d t=40 s图1 赤道孤立Rossby波算例模拟结果示意图Fig.1 Simulation result of equatorial solitary Rossby wave

算例共计算了40个单位时长，多项式阶数P=2，四周边界均视为墙面，采用可滑移边界条件，模拟结果如图1所示。由图可知，在孤立Rossby波向西传播过程中，有一部分水体分离出来以孤立Kelvin波的形式向东传播，这一现象是由于给定的初始条件不是精确解造成的。摄动解预测孤立Rossby波波峰相位速度为-0.77 m/s，最终会运动到x=-15.68处，波峰高度为0.153 4；Giraldo和Warburton运用间断有限元方法模拟的孤立Rossby波波峰相位速度为-0.78 m/s，最终会运动到x=-15.56处，波峰高度为0.148 4；本模型计算结果孤立Rossby波波峰相位速度为-0.78 m/s，最终会运动到x=-15.64处，波峰高度为0.158 4。本文模型模拟结果与摄动解和已有数值模型的结果均接近。

为了说明提高多项式阶数对计算精度的作用，分别计算了孤立Rossby波在多项式阶数P分别为1、2、3时的结果，三个算例的计算网格数均为1 600个。由于本算例没有精确解，本文选取25 600个网格下P=3的计算结果作为精确解，分别计算了不同阶数下模拟所得水位的二范数误差，如表1所示。由表可知，在相同的网格划分条件下，随着多项式阶数的提高，模拟的误差逐渐减小，模型的计算精度逐渐提高。

表1 不同阶数下孤立Rossby波模拟结果误差比较Tab.1 L2 errors of simulation results in equatorial solitary Rossby wave problem

2.2 平底坡Stommel问题

Stommel问题是科氏力、风应力和底摩阻达到平衡稳定状态时的解，它涉及到浅水方程实际应用中会遇到的多种物理过程，是海洋环流问题的一个很好的原型[29]。参考文献[29]的算例设置，算例中计算域大小为106×106m的矩形区域，计算域内水深均为1 000 m，四周边界视为壁面，采用可滑移边界条件。计算域共划分为502个任意四边形网格，算例选取多项式阶数P=2。本算例中科氏力系数f由式(23)给出，风场切应力由式(24)给出，底摩阻应力由式(25)给出

(23)

(24)

τbx=ργuh，τby=ργvh

(25)

2-a 模拟结果水位等值线2-b 模拟结果流速等值线图2 稳定状态水深等值线图及流速等值线图Fig.2 The free surface elevation and the velocity magnitude at the steady state of the nonlinear Stommel problem

式中:f0=10-4，β=10-111/m，L=106m；τ0=0.2 N/m2；ρ=1 000 kg/m3，γ=2×10-6，各参数选取均与文献[29]相同。

在科氏力、风应力和底摩阻的共同作用下，计算域内初始会在中心形成一个凸起的波，之后波峰慢慢向左侧偏移，当偏移一定距离后达到稳定状态，波的形状和位置均不再发生变化。结果达到稳定时计算域内水位等值线和流速等值线如图2所示，与文献[29]中的结果相近。

2.3 平缓岸滩波浪爬坡

为了说明本模型的干湿处理能力，考虑一个由海底滑坡引起的波浪在平缓岸滩上爬坡的算例。本算例长度为[-500 m,50 000 m]，底坡坡度为1/10，初始水位由式(26)给出

η=a1e-k1(x-x1)2-a2e-k2(x-x2)2
u=0,t=0

(26)

(27)

3-a 初始波面3-b t=160 s、175 s、220 s时爬坡和回落过程波面图3 波浪初始波面及模拟波面图Fig.3 The initial wave profile and the simulation result of wave profile at t= 160 s, 175 s and 220 s

初始波形如图3-a所示，初始速度为0。计算域划分为10 000个网格，计算时间为280 s，多项式阶数P=2。图3-b给出了t=160 s、175 s和220 s三个时刻的波面模拟结果。由图可知，本模型模拟的波浪爬坡(t=220 s)、回落(t=160 s、t=175 s)波面和Carrier等[30]给出的半解析解吻合良好，说明模型具有较为可靠的干湿处理能力。

2.4 斜坡静水算例

为了说明模型的和谐特性，本文选取了Bonev等[20]采用的斜坡静水算例。计算域为[-1.5 m,1.5 m]×[-1.5 m,1.5 m]的矩形区域，底坡b仅为x的函数，表达式如式(28)所示

(28)

4-a t=1 s时斜坡静水水位4-b t=1 s时斜坡静水x向通量图4 斜坡静水算例t=1 s时计算结果Fig.4 The free surface elevation and the flux magnitude at t=1 s of the water at rest on a sloping beach

静水面设置在z=0 m处，g=9.8，最小水深H0=1×10-4m，四周边界均设置为可滑移固壁边界，多项式阶数P=2。这里将未对底坡源项进行拆分的间断有限元浅水模型称为传统方法，图4-a分别绘制了模拟至1 s时本模型和传统方法计算的水位，由图可知，本模型计算的水位基本保持为静水位而传统方法计算的水位则在干湿交接附近有微弱波动。图4-b分别绘制了模拟至1 s时本模型和传统二维浅水方程计算的x向通量，此图可以更直观地看出传统方法存在虚假流动的问题，而本模型通过对底坡源项进行拆分，保证了各项阶数的平衡，从而很好地保持了斜坡上的静水特性。

3 模型应用

5-a 网格划分5-b 测点位置图5 三亚红塘湾算例网格划分及测点布置示意图Fig.5 Computing grid and locations of measuring stations of Hongtang Bay in Sanya

为了验证加入源项后的二维水动力模型模拟实际潮流的合理性，选取三亚市红塘湾附近海域，模拟了2016年4月26日～27日一个完整大潮期间的潮流运动。本算例计算域共划分为2 448个任意四边形网格，如图5-a所示，网格划分由外海至近岸逐渐加密。由于该海域受北部湾复杂潮流影响，为了充分反映模拟海域的潮流运动，本文首先选取北部湾区域建立大模型，运用Mike21的FM模块模拟了大潮期间北部湾大范围的潮流运动，本算例开边界使用的潮位数据由大模型计算结果插值得到。本算例中科氏力系数根据网格点的实际纬度得到，糙率n根据当地泥沙粒径分布及地形特点确定取值为0.017，计算阶数P=1。

图6 2016年4月26日～27日大潮期潮位模拟值与实测值对比Fig.6 Comparison of surface elevation between simulation results and observed data

从2016年4月26日～27日的水文实测资料与计算结果进行验证，选取大潮T1、T2两个测站的潮位资料，V1～V6共6个测站的流速流向测站资料，潮位及流速测站位置如图5-b所示。潮位验证结果如图6所示，流速验证结果如图7所示，流向验证结果如图8所示。由图可知，模型模拟的潮流运动与实际观测的情况吻合良好，模型计算结果较为可靠。

图7 2016年4月26日～27日大潮期流速模拟值与实测值对比Fig.7 Comparison of current speed between simulation results and observed data

图8 2016年4月26日～27日大潮期流向模拟值与实测值对比Fig.8 Comparison of current direction between simulation results and observed data

4 结论与展望

基于节点无积分间断有限元和任意四边形网格建立二维浅水方程模型，在初始控制方程中加入了科氏力、风应力和底摩阻项，使得模型可以模拟风吹流及实际潮流运动。采用赤道孤立Rossby波算例验证了科氏力项的合理性，通过求解Stommel问题的稳定解对科氏力、风应力和底摩阻共同作用下的模型进行了验证。在源项验证完毕后，利用平缓岸滩波浪爬坡算例和斜坡静水算例分别说明了模型的干湿处理能力和静水模拟的和谐特性。模型最后应用于三亚市红塘湾海域，对实际大潮过程进行了模拟，模拟结果与全潮水文观测的水位和流速、流向结果吻合良好。本文建立的模型已经具备了模拟河口海岸二维水动力的基本功能，今后将进一步在模型中加入粘性应力项及辐射应力项，完善波流耦合功能进而实现实际情况下河口海岸地区风、浪、潮相互作用的风暴潮过程模拟。