基于无积分节点间断有限元的二维水沙模型：(2)泥沙运动与地形演变

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

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

泥沙运动对于河口海岸演变、航道淤积以及河口治理等问题有着重要的影响，一直是河口海岸地形演变过程的研究重点。近年来，随着计算机技术的发展，数值模型越来越成为研究泥沙运动的重要手段。目前，二维水沙模型在解决工程泥沙问题中的应用已十分普及。过去几十年来，有限差分法、有限元法和有限体积法等经典数值方法均在二维水沙模型中得到大量应用。有限差分法建立在经典数学逼近理论的基础上，简洁直观，但差分法一般不易保证计算的守恒性[1-2]；有限单元法计算精度相对较高，但不适于计算间断问题[3-4]；有限体积法在实际应用中，可以保证不同网格下的积分守恒，能够计算间断问题，但很难在一般非规则网格上取得高精度[5-6]。与上述三种经典方法相比，间断有限元法(Discontinuous Galerkin Method)因在数值精度、间断捕捉能力、网格自适应特性、紧致性、高度可并行性以及适用于非结构化网格等方面的优势，近年来逐渐开始在二维水沙模型的研究中得到应用[7]。

Gourgue[8]较早采用间断有限元方法建立了二维水沙数值模型，并应用于Scheldt河口泥沙运动研究。石宝海和陈焕贞[9]基于迎风间断有限元方法建立了平面二维水沙输运模型，泥沙运动采用挟沙力模型描述，并给出了一个理想算例结果。赵张益和张庆河[10](2014)采用间断有限元法建立了具有二阶精度的非结构化网格二维悬沙数值模型，模型能够有效限制间断处的数值振荡，保证稳定性，但是该模型未考虑推移质输沙和地形演变的影响。

上述研究表明，间断有限元可以用于二维水沙运动的模拟，然而，这些模型中包含的正交项是基于全阶积分的，这需要大量的计算资源[11]。为了减少这种计算量，Atkins和Shu[12]提出了无积分节点间断有限元格式，该格式避免了离散方程时的积分过程，能够有效减少计算时间以及存储空间，并可以保持间断有限元固有的紧致性和稳健性，具有良好的应用前景，但是该方法要求单元的雅克比系数为常数，不能适用于任意四边形网格。李龙翔和张庆河[13]在此基础上建立了适用于任意四边形网格的无积分节点间断有限元模型，用矩阵运算进行数值积分，有效地减小了四边形网格上模型计算所需的存储空间与计算量。李文俊等[14]在二维浅水方程中加入了科氏力、风应力、底摩阻项，基于无积分节点间断有限元建立了充分反映各种实际条件影响的二维水动力模型。

本文将在李文俊等[14]二维水动力模型基础上，进一步基于无积分节点间断有限元建立二维泥沙运动与地形演变模型，同时，本模型采用了李龙翔和张庆河[15]提出的采用权重重构方法的节点型斜率限制器来消除间断处产生的振荡，为二维高效高精度无积分节点间断有限元二维水沙模型的建立奠定基础。

1 二维泥沙运动及地形演变数学模型

1.1 水动力方程

模型的水动力计算采用二维浅水方程作为控制方程，该方程是由三维N-S方程在采用静压假定，忽略粘性以及垂向加速度的情况下，沿水深积分得到的

(1)

式中：U= [h,hu,hv]T，F(U)=[E(U),G(U)]T，S(U)分别代表守恒变量、通量项以及源项，表达式写为

(2)

(3)

式中：h为总水深；g为重力加速度；u和v分别为x方向和y方向的垂线平均流速；z为底部高程；η=h+z为水位；Sfx和Sfy分别为x方向和y方向底摩阻项，可表示为

Sfx=n2u(u2+v2)1/2h-4/3，Sfy=n2v(u2+v2)1/2h-4/3

(4)

式中：n为曼宁系数。

1.2 泥沙运动控制方程

泥沙运动暂时只考虑非粘性泥沙，按照其运动型式可分为推移质和悬移质，下面分别给出描述泥沙运动的推移质运动公式和悬移质运动控制方程。

1.2.1 推移质运动控制方程

推移质运动控制方程一般用经验公式表示，本文采用van Rijn[16](2007)公式

qbx=0.015uh(d50/h)1.2Me1.5qby=0.015vh(d50/h)1.2Me1.5

(5)

式中：qbx和qby分别为x和y方向的推移质泥沙浓度；d50为泥沙中值粒径；Me表达式为

Me=max(0，|uc|-ucr,c)/[(s-1)gd50]0.5

(6)

式中：uc为水平方向的合速度；s为泥沙密度与液体密度的比值；ucr,c为泥沙运动的临界流速，表达式为

(7)

式中：d90表示一个样品的累计粒度分布数达到90%时所对应的粒径。

1.2.2 悬移质运动控制方程

悬移质泥沙运动控制方程采用沿水深平均的二维对流扩散方程

(8)

式中：c为沿垂向平均后的泥沙浓度；εx和εy为泥沙的水平扩散系数；Fs为源项。

源项可以采用挟沙力公式或参考浓度法描述，本文采用参考浓度法[17]进行计算

Fs=P-D=(ca-c0)wf

(9)

式中：P为泥沙起悬通量；D为泥沙沉积通量；wf为泥沙沉速。ca和c0是定义在参考高度处a的浓度，ca与底床切应力有关，表示悬沙平衡时参考高度处的泥沙浓度，c0与对流扩散方程的垂向平均浓度有关，表示在水深平均浓度对应于参考高速度处的浓度。当ca>c0时，泥沙起悬，当ca

(10)

式中：τs,max为最大底部应力；τcr为泥沙颗粒临界启动切应力；βd为垂向浓度分布转换系数，由下式计算

(11)

(12)

式中：κ为冯卡门系数，一般取0.4；u*c为底摩阻流速。

1.3 地形演变方程

推移质和悬移质的泥沙浓度会引起底床的冲刷和淤积变形，底床的地形演变可通过以下方程描述

(13)

式中：p为泥沙的孔隙率。

2 控制方程数值离散

水动力模型中浅水方程的数值离散过程已在文献[14]中给出，下面重点给出悬沙输移对流扩散方程与地形演变方程的离散过程。

2.1 对流扩散方程与地形演变方程的离散

式(8)中含有2阶偏导项，为将方程改为一阶偏导方程，引入辅助变量[19]

Q=εhc

(14)

此时二维对流扩散方程可改写为

(15)

为方便离散，将(13)、(15)式写为

(16)

式中：W为守恒通量，W=[hc,z]T，P(W)=[K(W),L(W)]和R(W)分别代表守恒变量通量项以及源项，表达式写为

(17)

(18)

与水动力模型的离散过程类似[14]，将求解域划分为Ne个非重叠单元Ωe[6]，定义Ωe上至多为p阶的多项式空间Vp(Ωe)，取Lagrange插值函数作为试验函数φVP(Ωe)，将试验函数乘以式 (14)和(16)，并在单元Ωe上积分，可得

(19)

(20)

对式(18)进行分部积分，得到离散方程

(21)

式中：n为本单元边界处的外法向向量，由于相邻单元在交界面两侧的值不同，为保证两侧单元流入与流出的通量守恒，引入数值通量(hc)*，采用中心通量计算[20]

(hc)*=0.5(hc-+hc+)

(22)

(23)

2.2 时间递进

对控制方程离散后的强形式进行进一步离散，得到如下的半离散格式

(24)

式中：uh为未知变量的矢量。

式(28)的时间递进可以采取欧拉法、龙格库塔法等求解。为了尽可能减小时间离散误差，本文采用显式四阶龙格库塔方法[22]用于半离散方程的时间递进。

(25)

式中：系数ai，bi，ci的取值参照文献[22]。

3 模型验证与应用

下面分别利用前人进行的泥沙悬扬、沙波演变和溃坝波冲淤水槽试验对建立的泥沙运动与地形演变模型进行验证。最后，还将模型应用于海南省红塘湾附近海域的模拟，与现场实测悬沙资料进行比较。

3.1 泥沙悬扬算例

该算例为van Rijn[23]在图1所示的水槽中所做的泥沙悬扬实验。水槽中分为定床段和铺沙段，铺沙段长30 m、宽0.5 m、高0.7 m。水深为0.25 m，平均流速0.67 m/s，底床泥沙中值粒径为0.23 mm。Lin 和 Falconer[24]曾用数学模型对此实验进行了模拟，本文模拟时的相关系数均参考文献[24]进行取值。悬沙沉降速度取为0.022 m/s，水平扩散系数取为0.002 m2/s。

图2为模型模拟结果与实测垂向平均含沙量结果的对比，结果表明，本文模型的数值解与实测值基本接近，泥沙输移趋势相同，模型可以有效地模拟泥沙悬扬以及输移过程。

图1 实验布置图Fig.1 Layout of experiment图2 垂向平均含沙量模拟值与实测值对比Fig.2 Comparison of simulated and measured vertical average sediment concentration

3.2 沙波演变算例

该算例在不考虑悬移质泥沙条件下模拟沙波地形演变过程[25]，用以检验模型模拟地形演变数值模型的合理性。初始沙波为正弦形，源项取为0，推移质表达式取为

(26)

式中：a和b为常数；u=(ux，0)为沿水深平均的流速；D=(Dx，0)为流量；Δy=1.2 m为宽度；a=0.001 s2·m-1；b=3；Dx=1 m3·s-1，初始地形的函数可参考一维算例[26]。

图3为y=0.75 m处剖面在t=500 s时模型模拟结果与文献[25]给出的解析解的对比，结果表明，本文模型的模拟值与解析解吻合较好，模型可以精确地模拟地形演变过程。图4为本文模型计算结果与文献[25]基于有限体积方法给出考虑人工扩散项的计算结果与解析解误差的对比图，可以看到，本文的模拟值误差更小，精度更高。

图3 沙波演变模型模拟值与解析解的对比Fig.3 Comparisons of simulated values and analytical solutions图4 间断有限元与文献[25]考虑扩散项有限体积法误差对比Fig.4 Error comparisons of discontinuous Galerkin method and finite volume method with additional diffusion in reference[25]

3.3 溃坝动床算例

图5 溃坝算例的实验布置图 (m)Fig.5 Layout of dam break experiment

该算例出自Leal的溃坝实验[27]，实验在一个长19.2 m、宽0.5 m、高0.7 m的水槽中进行，图5为初始时实验布置的示意图，水槽分为左右两端，左段在固定床面之上铺沙高度为0.19 m，水深0.4 m，右侧铺沙高度为0.071 m，水深0.075 m。泥沙中值粒径为1.22 mm，泥沙沉速为9.92 cm/s。

图6和图7分别显示了t=1 s和t=4 s时模型模拟与实验测量水位值和地形值的比较情况。虽然水位、地形模拟值与实测值在局部有一定差异，二者变化趋势总体上吻合较好，可以认为模型能够合理模拟复杂水流及其引起的泥沙冲淤导致的快速地形演变问题。

图6 t=1 s时模拟值与实测值的对比Fig.6 Comparisons of simulated and measured values at t=1 s图7 t=4 s时模拟值与实测值的对比Fig.7 Comparisons of simulated and measured values at t=4 s

3.4 模型在三亚红塘湾泥沙运动模拟中的应用

图8 三亚红塘湾泥沙测站布置图Fig.8 Layout of sediment station in Hongtang Bay of Sanya

为了进一步考察二维水沙模型模拟实际海域水沙运动的合理性，本节选取三亚市红塘湾附近海域，模拟了2016年4月26日～27日一个完整大潮期间的泥沙运动。本算例计算域与文献[9]相同，潮流模拟结果与实测资料的比较见文献[14]，这里主要给出悬沙运动模拟结果。根据现场实测资料分析，悬沙中值粒径取为0.3 mm，悬沙沉降速度取为0.035 m/s。

现场全潮水文观测包含潮流与悬沙测量，这里选取S1～S6共六个测站的悬浮泥沙资料与模拟结果进行比较。泥沙测站位置如图8所示。模拟结果与实测结果比较见图9。由图可知，模型模拟的泥沙运动与实际观测情况吻合良好。可以认为，本文建立的水沙模型能够合理模拟比较复杂的现场实际情况。

图9 2016年4月26日～27日大潮期悬沙浓度模拟值与实测值对比Fig.9 Comparison of simulated and measured suspended sediment concentration during high tide period (2016-04-26～27)

4 结论

基于无积分节点间断有限元方法建立了泥沙输移二维数学模型，悬沙运动通过对流扩散方程进行离散求解，方程的源项采用参考浓度法，推移质输沙采用经验公式进行计算，地形演变过程则通过推移质和悬移质输沙控制的床面变形方程进行求解。建立的泥沙输移与地形演变数值模型获得了已有水槽冲刷试验、沙波演变解析解和溃坝动床实验的验证。模型还应用于三亚红塘湾大潮过程悬沙运动的模拟，模拟结果与全潮水文观测多点含沙量过程吻合较好。模型验证和应用结果表明，基于无积分节点间断有限元建立的水沙模型可以合理模拟泥沙运动及其导致的地形演变问题。

为了更全面描述河口海岸泥沙二维运动问题，模型将进一步考虑波流共同作用下二维水动力及泥沙运动。