NUMERICAL STUDY ON SUSPENDED SEDIMENT TRANSPORT UNDER STRONGLY NONLINEAR WAVE INDUCED BY SHALLOW WATER
-
摘要:
波浪作用下的泥沙运动是海岸工程的研究热点；波形效应和表面波效应是公认的影响近岸泥沙输运的原因。由于构建大尺度波浪水槽和相应造波的困难，波浪作用下的泥沙输运实验和数值研究大多基于一维往复流水槽开展，波浪水槽的结果较少。与传统往复流水槽不同的是，波浪水槽可以同时考虑波形效应和表面波效应对泥沙输运的影响，在尺度上更加接近真实情况。该文构建了适用于浅水强非线性波输沙且计算成本低的数值波浪水槽，采用边界造波法造波和阻尼法消波，由VOF模型追踪自由表面。将体现了相位差、质量守恒、加速度效应和边界层非对称发展的泥沙近底边界条件拓展至二维波浪水槽，避免了两相模型对近底高含沙区域颗粒碰撞、摩擦以及动床面捕捉的复杂计算。此外，利用已有成果在泥沙运动中考虑了颗粒尾流涡效应和非静水高含沙的相对速度理论修正。该波浪水槽成功模拟了强非线性二阶椭圆余弦波作用下的泥沙输运，可靠性得到了大尺度波浪水槽实验及包含自由表面波的两相模型结果验证。数值波浪水槽同往复流水槽的结果对比显示，在自由表面波影响下存在额外的向岸流使向岸输沙显著增强。
Abstract:Sediment transport under nonlinear waves is a hot research topic in coastal engineering. Wave profile and free surface wave effects are recognized as important factors in sediment transport. Due to the difficulty of constructing large-scale wave flume and corresponding wave generation, most experiments and numerical studies on sediment transport are carried out upon an oscillating water tunnel, and few results are obtained from wave flume. In comparison with a traditional oscillating water tunnel, the wave flume performed both wave profile and free surface wave effects on sediment transport in nearly real scale situation. A numerical wave-sediment flume is constructed for strongly nonlinear propagation of wave over sediment bed with low computational work. The nonlinear wave is generated by boundary wave maker and absorbed by a damping function, and classical VOF method is used to handle the free surface. The near-bed sediment condition for asymmetric wave is extended to the two-dimensional flume, which introduces phase lag, mass conservation, acceleration effect and asymmetric boundary layer development. Therefore, complex calculations for near-bed high sediment concentration in a two-phase model are avoided, i.e., the computation of particle collision, friction, and detection of mobile bed surface. In addition, the particle wake vortex effect and the relative velocity modification of non-static high sediment concentration are considered in sediment transport. The numerical wave flume successfully simulates the sediment transport under strong nonlinear second-order conical waves, and its reliability is verified by the large-scale wave flume experiment and the two-phase model including free surface waves. Comparison of the results of the numerical wave flume with the oscillating water tunnel shows that the presence of additional onshore streaming under the influence of free surface waves significantly enhances the onshore sediment transport.
-
在沿海地区，波浪及其净流动对地形变化有很大的影响，涉及了岸滩演变、港口航道冲刷淤积、海岸建筑物的局部冲刷、海洋生态环境变化等问题。尤其在大风暴来临或强非线性波浪情况下，泥沙颗粒快速悬浮并随水流在海床底部大量输移，可能在短时间内造成巨大的海底形态演变[1]。波浪自深水传播至浅水区时，近底水流对泥沙颗粒剪切应力增加，水中泥沙颗粒悬浮增加并产生向岸或离岸泥沙净输运[2-4]。由于近岸复杂的水动力环境以及水沙相互作用机制，近岸泥沙运动一直是具有挑战性的热点问题。其中，波浪非线性和表面波效应是公认的影响泥沙净输运的原因[3]。
波浪在近岸区的传播会发生明显的浅水变形，波浪形态逐渐变为具有高而窄的波峰和低而平的波谷。对中等或粗泥沙颗粒观察到了向岸的泥沙输运。然而对细颗粒泥沙，许多往复流水槽(Oscillating Water Tunnel, OWT)实验观察到了离岸净泥沙输运，这可归因于相位滞后效应[4]。向岸半周期卷起的泥沙在流速转向时刻并未完全落下，并在离岸半周期被继续输运，由于离岸周期较长进而导致离岸泥沙净输运[5-6]。在波浪破碎附近，破碎变形会导致波浪具有加速度偏度，通常会出现锯齿波。由于加速度偏斜，近底边界层厚度非对称导致向岸净输运。OWT实验表明在这种流动下，泥沙净输移是向岸的[7-9]。非线性波形引起泥沙净输运的这一机制又称为波形效应[3]。
OWT因其可以获得较大的底床切应力，在往复流驱动泥沙输运方面已经取得了很多成果。但理想化的OWT和真实波浪存在很大差异。真实波浪条件下近底边界层内水流的水平和垂向速度分量的相位差并非π/2，垂向速度分量中存在与水平速度分量同相位的成分必然引起一个二阶的向岸边界层净流动。LONGUET-HIGGINS[10]最早在层流边界层内提出这一现象。KRANENBURG等[3]将这一边界层向岸净流动归因于表面波效应。基于OWT的研究可以很好的考虑波形效应对泥沙输运的影响。然而其忽略了垂向流速，因此并不能完全表述真实波浪作用下的泥沙输运。在过去十几年里，大尺度波浪水槽中进行了几次强非线性规则波驱动的泥沙输运实验[4, 11]，DOHMEN-JANSSEN等[4]以及SCHRETLEN[11]发现在真实波浪作用下，近岸泥沙的输沙率约为同等条件下OWT内观测结果的2.5倍；并且二阶斯托克斯波作用下的细沙输运并不是OWT中测得的离岸输运，而是向岸输运。这说明在近岸区域真实波浪作用下，表面波效应引起的边界层净流动不可忽略且同波形效应相互竞争，共同决定泥沙输移。
传统的OWT数值模型忽略了垂向流速，将水平方向均匀近似并把外边界层的自由流速度转换为水平压力梯度来驱动数值模型[12]。考虑到求解包含自由表面流动问题的复杂性和数值不稳定性，部分学者将水平速度梯度转换为时间导数进而将表面波效应包含进OWT数值模型中[13-16]。KRANENBURG等[3, 15]证明了在细颗粒泥沙情况下表面波效应产生的边界层向岸净流动能够将净输沙从离岸(OWT)逆转为向岸(波浪水槽)。TAN等[16]利用上述方式考虑表面波效应，提出了一种精度与两相模型相似的两层泥沙输移数值模型。目前已知这种将水平速度梯度转换成时间导数的方式在小振幅波情况下是适用的，但对非线性波的适用性还不清楚，因为非线性波常伴随着波形的快速变化[17]。KIM等[17]基于OpenFOAM框架开发出一种名为SedWaveFoam的模型，该模型集成了泥沙输运的欧拉两相模型SedFoam和表面波求解器wave2Foam，通过构建数值波浪水槽模拟包含自由表面流动的泥沙输运问题。经与大尺度波浪水槽实验对比，发现其可以很好的描述包含自由表面效应的泥沙输运过程。
随着计算机快速发展，物理模型实验、理论解析和数值模拟相结合的方式逐渐成为了研究波浪运动及泥沙输运问题的重要手段[18-19]。然而数值波浪水槽涉及造波、消波和自由面追踪，需要比较大的计算域；同时两相模型涉及复杂的颗粒应力、两相作用力、紊流计算和相间迭代，在高含沙区域的网格尺度一般需与泥沙粒径一致。基于数值波浪水槽的两相模型往往计算量很大，目前并不适合大尺度的工程运用[20-22]。波浪运动下泥沙模拟成本的降低绕不开精确的泥沙近底边界条件，它涉及由波浪的非恒定流属性引起的诸多复杂因素[23-24]；而目前广泛使用的泥沙近底边界条件仍然主要依靠恒定流推导的结果。因此，CHEN等[23-24]考虑了泥沙运动相对波浪的相位差(相位漂移、相位残留)、质量守恒、波浪破碎带的加速度效应、波浪浅水变形导致的边界层非对称性，获得了适用于非对称歪斜波输沙的泥沙近底边界条件。进一步地，LIANG等[21]通过CHEN等[23-24]的边界条件建立混相模型应用于层移输沙问题，将计算量减少了一个数量级，同时保持了两相模型的高精度。
从进一步减少数值波浪输沙水槽计算量的角度考虑，发展合适的泥沙输移扩散模型是一个较为理想的选择，具有较为广阔的工程应用前景。它相比混相模型可以进一步减少混合物相对速度的迭代、两相紊动及相关应力计算。因此，本文基于泥沙的输移扩散模型构建适用于非线性波的数值波浪输沙水槽，吸收LIANG等[21]混相模型对悬沙输移扩散性能的考虑以及CHEN等[23-24]发展的适用于非对称歪斜波输沙的泥沙近底边界条件，进而探讨包含真实表面的浅水强非线性波作用下的泥沙输运。
1 数值模型
本文讨论的是二阶椭圆余弦波代表的浅水强非线性波作用下的泥沙输运，如图1所示。其中：D为水深；L为波长；H为波高；c为悬沙浓度；x为笛卡尔坐标，下标i=1、2。
![]()
图 1 浅水强非线性波作用下泥沙输运示意图Figure 1. Sketch of sediment transport in strongly nonlinear waves1.1 基本方程
雷诺平均的连续性方程如式(1)：
∂ui∂xi=0 (1) 动量方程如式(2)：
∂ui∂t+∂(uiuj)∂xj=∂∂xj(−¯u′iu′j)−∂pρ∂xi−S(xi)+∂∂xj[v(∂ui∂xj+∂uj∂xi−23δij∂ul∂xl)]+gi (2) 式中：ui为xi方向的水流速度；t为时间；ρ为水的密度；p为压强；v为运动黏度；δij为克罗内克尔符号；下标l为求和指标；−¯u′iu′j为雷诺应力项；gi为重力加速度；S(xi)为阻尼耗散项。
悬沙运动输运方程如式(3)：
∂c∂t+∂∂xi[(ui−δi2ws)c]=∂∂xi(vtSc∂c∂xi) (3) 式中：ws为泥沙沉速；vt为水流涡粘系数；Sc为Schmidt数，泥沙扩散系数γ定义为vt/Sc。颗粒尾流涡效应主要体现在Sc的给定方法上，前人的研究大多将Sc与w0/u*建立关系，认为Sc随着w0/u*在一定范围内增加而减小[25-27]，w0为单颗粒泥沙静水沉速，u*为摩阻流速。SHI等[20]认为随着粒径增长，泥沙颗粒惯性变大，水沙相对运动越加明显，颗粒与绕流流体的相互作用变得更加复杂。颗粒局部绕流形成的尾流区内紊流增强，相应的颗粒脉动也得到增强。对于粗粒径的泥沙，尾流促进颗粒扩散效应更明显，泥沙紊动扩散会更加符合实际[20]。本文采用LIANG等[21]基于SHI等[20]的工作进一步得到的结果，如式(4)：
Sc=(1+1.5C4/3Dψ1/31w20/k)−1 (4) 式中：k为水相紊动能；CD为泥沙拖曳力系数，定义为CD=24(1+0.15Re0.687s)/Res，Res为泥沙颗粒雷诺数，定义为Res=w0ds/v，ds为泥沙中值粒径；ψ1为修正系数，如式(5)：
ψ1={1−exp[−0.007(Res−2)],Res⩾ (5) 考虑到细沙(Res<2)可以较好的响应紊流流动，忽略颗粒绕流产生的局部紊动能增强(ψ1=0)是合理的[20-21]。高含沙时，泥沙颗粒群对颗粒沉降产生阻碍，因此泥沙沉速为：
{w_{\rm{s}}} = {(1 - c)^n}{w_0} (6) 式中，n为理论推导的与泥沙颗粒雷诺数有关的参数[21]，n=2.2+2.65/(1+0.06{{Re}}_{\rm{s}}^{0.7})。式(4)~式(6)[21]利用颗粒惯性及颗粒绕流的相对重要性决定扩散系数大小，并且通过沉速修正有效表示了非静水条件下高浓度泥沙与水的相对运动速度。
紊流闭合使用标准k-ε模型，方程如式(7)~式(9)：
\frac{{\partial k}}{{\partial t}} + \frac{\partial }{{\partial {x_i}}}\left( {k{u_i}} \right) = \frac{\partial }{{\partial {x_i}}}\left[ {\left( {v + \frac{{{v _{\rm{t}}}}}{{{\sigma _k}}}} \right)\frac{{\partial k}}{{\partial {x_i}}}} \right] + {G_k} - \varepsilon (7) \frac{{\partial \varepsilon }}{{\partial t}} + \frac{\partial }{{\partial {x_i}}}\left( {\varepsilon {u_i}} \right) = \frac{\partial }{{\partial {x_i}}}\left[ {\left( {v + \frac{{{v _{\rm{t}}}}}{{{\sigma _\varepsilon }}}} \right)\frac{{\partial \varepsilon }}{{\partial {x_i}}}} \right] + {C_{1\varepsilon }}\frac{\varepsilon }{k}{G_k} - {C_{2\varepsilon }}\frac{{{\varepsilon ^2}}}{k} (8) {v _{\rm{t}}} = {C_\mu }\frac{{{k^2}}}{\varepsilon } (9) 式中：ε为紊动能耗散率；Gk为梯度项产生的紊动能；模型常数为C1ε=1.44，C2ε=1.92，Cμ=0.09，σk=1.0，σε=1.3。
水相自由表面压强等于大气压。自由面变化属于水气两相流动，本文选用VOF方法，求解单一的动量方程并跟踪区域内流体的体积分数F(单元中空气占网格体积比例)[28]。F=0表示该单元为水所占据；F=1表示该单元被空气所占据；若0<F<1，则表示该单元为交界面。F满足的方程如式(10)：
\frac{{\partial F}}{{\partial t}} + \frac{\partial }{{\partial {x_i}}}\left( {F{u_i}} \right) = 0 (10) 1.2 水相边界条件
计算域底部为无滑移壁面，顶部压力入口。入口边界使用在数值波浪研究工作中广泛使用的边界造波法来构建数值波浪水槽[29-32]。入口给定椭圆余弦波的波面、水平流速、垂向流速，F在波面以下为0，波面之上为1。LAITONE[33]基于摄动法推导出了经典的椭圆余弦波二阶解，将H/Dt看作扰动参数，Dt为波谷处水深。相较于一阶解，其考虑了高阶项对结果的影响，主要表现在压力变化不再是静压，垂向速度和加速度在量级上不再是小值。而ISOBE等[34]参考了LAITONE[33]，选定H/D作为扰动参数推导出了二阶解。在已知波高、水深、周期等参数的情况下，相对一阶解得到的波形非线性更强，在一定程度上适用性更广且节省了更多的计算工作。二阶解的波面位置为：
N = D({A_0} + {A_1}{c_1} + {A_2}{c_2}) (11) 水平流速为：
\begin{split} {u_1} = & - \sqrt {gD} \left\{ {\frac{1}{2}{{\left( {\frac{{{x_2}}}{D}} \right)}^2}({B_{01}} + {B_{11}}{c_1} + {B_{21}}{c_2})} \right\} + \\& \sqrt {gD} ({B_{00}} + {B_{10}}{c_1} + {B_{20}}{c_2}) \end{split} (12) 垂向流速：
\begin{split} {u_2} = &- \sqrt {gD} \frac{{4KD}}{L}{\rm{csd}}\frac{1}{6}{\left( {\frac{{{x_2}}}{D}} \right)^3}({B_{11}} + 2{B_{21}}{c_1}) + \\& \sqrt {gD} \frac{{4KD}}{L}{\rm{csd}}\left( {\frac{{{x_2}}}{D}} \right)({B_{10}} + 2{B_{20}}{c_1}) \end{split} (13) 椭圆余弦波理论涉及到椭圆函数使其求解过程十分复杂。为了简化计算需要根据式(14)迭代计算出椭圆模数m：
\frac{{16{m^2}{K^2}}}{3} = \frac{{gH{T^2}}}{{{D^2}}}\left( {1 + \frac{H}{D}\left( {\frac{{ - 1 - 2\lambda }}{4}} \right)} \right) (14) K = \int_0^{\pi /2} {\frac{{{\rm{d}}\theta }}{{\sqrt {1 - {m^2}{{\sin }^2}\theta } }}} (15) E = \int_0^{\pi /2} {\sqrt {1 - {m^2}{{\sin }^2}\theta } {\rm{d}}} \theta (16) 式中：N为波面；T为周期；m为椭圆模数；K为第一类完全椭圆积分；E为第二类完全椭圆积分；A0、A1、A2、c1、c2、csd、λ、A0、A1、A2、B00、B10、B20、B01、B11、B21为系数，具体可参考文献[34]。在已知波高、水深、周期等参数下，一旦求得椭圆模数m就可以计算第一类完全椭圆积分K、第二类完全椭圆积分E以及椭圆函数，进而得到波面和流速表达式。
为了避免出口处波浪反射对前端计算域的影响，在计算域末端设置1倍~2倍波长的消波区。本文将消波边界设置为压力出口，选用阻尼消波法消波。其原理在于在动量方程中添加阻尼源项使ui按规律衰减。x1方向阻尼源项表达式为：
S\left( {{x_1}} \right) = \alpha \frac{{{x_1} - {X_0}}}{{{X_1} - {X_0}}}{u_1},\;\;\;\; {X_0} \leqslant {x_1} \leqslant {X_1} (17) 式中：α的取值与ui有关，本文取8；X0和X1为消波起始点和终点的水平坐标位置。x2方向阻尼源项形式同上，α同取8。
1.3 泥沙相边界条件
参照LIANG等[21]，自由表面处给定不可穿透条件，进口和出口泥沙浓度均设置成无梯度条件。表面波效应带来的额外泥沙向岸净输移主要发生初始床面以上区域[17]，因此本文重点关注初始床面以上的泥沙输运。近底泥沙浓度边界条件的给法通常有两种：一种是在某一参考高度处给定参考浓度；另一种是给定扬沙率，本文使用常见的参考浓度法。考虑到非线性波切应力的加速度效应、边界层紊流非对称发展的影响以及泥沙运动的相位差现象和质量守恒，参考浓度采用对于非对称歪斜波具有很好适用性的CHEN等[23-24]公式：
c({x_1},t) = {c_{\rm{m}}}\exp \left[ { - \left(1 + \frac{{2{d_{\rm{s}}}}}{Z}\right)} \right] (18) \frac{{Z({x_1},\sigma t + \psi )}}{{{d_{\rm{s}}}}} = ({\beta _1}{\theta _{\rm{m}}} + {\beta _2}\theta )\left(1 + 0.011\frac{{{u_{\rm{m}}}}}{{{w_0}}}\right) (19) \theta ({x_1},t) = \frac{{{\tau _{\rm{b}}}}}{{\rho g(s - 1){d_{\rm{s}}}}} (20) 式中：cm为泥沙体积浓度最大值[6]，取0.6；Z为侵蚀深度；τb为床面切应力；σ为频率，σ=2π/T；s为泥沙比重，s=2.65；ψ为相位滞后参数，ψ=2σZm/w0，Zm为最大侵蚀深度；θ为希尔兹数；θm为最大希尔兹数；β1、β2与相位滞后参数ψ有关，β1θm为流速转向时的泥沙残留，定义为β1=3.2exp(−0.2/ψ)，β2=3.2−β1；um为水质点速度极值。将本文的近底泥沙参考浓度同基于恒定流动提出的ZYSERMAN和FREDSØE[35](ZF94)参考浓度进行对比。
往复流水槽[23-24]只需确定一个位置的参数。数值波浪水槽受沿程切应力和消波等因素的影响，不同位置的波能以及周期切应力变化不同，这就导致不同的相位漂移和相位残留，进而参考浓度也存在一定差异。为了拓展CHEN等[23-24]近底泥沙边界条件到波浪水槽，计算中记录了周期内近底沿程的床面切应力。每个位置根据其最大切应力求得最大希尔兹数和最大侵蚀深度，在计算中调用瞬时切应力得到瞬时希尔兹数，进而得到不同位置处的相位滞后、相位残留以及侵蚀深度，根据式(18)计算出准确的参考浓度，之后周期迭代达到收敛。
2 结果与讨论
2.1 计算模型与参数
模型验证采用DOHMEN-JANSSEN等[4]在德国GWK大型波浪水槽进行的强非线性规则波实验，工况对应于实验中mh组次。基本参数为H=1.6 m、T=6.5 s、D=3.5 m，厄塞尔数Ur=54.17，表征波浪非线性强弱。该波浪水槽长300 m，宽度为5 m，高度7 m。实验在计算域85 m~130 m段设置了0.75 m厚的沙层，ds为0.24 mm。实验中水相流速使用声学多普勒测速仪在x1=109 m断面处、x2=0.109 m处测得，实验波浪形态及流动特征形似椭圆余弦波。为节约计算成本，本文参考了KIM等[17]的设定，计算域参数取长280 m、高7 m、水深3.5 m，沙层顶部x2=0处为初始床面。对于不同波浪条件，KIM等[17]出于数据的可用性选用H=1.55 m、T=6.5 s为基本参数，采用十阶流函数波匹配工况mh的流速分布。本模型选用ISOBE等[34]推导的二阶椭圆余弦波理论以波高H=1.6 m、T=6.5 s为基本参数匹配实验的水平流速。
计算域网格均采用长方形网格，横向网格宽度为均匀的0.1 m，纵向网格划分以静水面上下1 m为界分为三个区域。泥沙运动主要集中在近底，所以在近底区域局部加密。下部区域近底第一层网格高度为4ds，以匹配参考高度2ds。为节约计算量向上采用比例为1.06的等比网格。静水面上下1 m区域及上部空气域采用0.04 m的均匀纵向网格，保证网格长宽比小于5∶1。
2.2 模型验证
本文采用一致性系数IA来衡量计算值与实验值之间的偏差[17, 36-37]。IA为无量纲数，当IA=0时意味着完全不一致，IA=1时意味着完全一致。IA表达式定义如下：
IA = 1 - \frac{{\displaystyle\sum\limits_{j = 1}^n {{{( {{P_j} - {O_j}} )}^2}} }}{{\displaystyle\sum\limits_{j = 1}^n {{{( {| {{P_j} - \overline O } | + | {{O_j} - \overline O } |} )}^2}} }} (21) 式中：j为数据点；Oj为实验值；Pj为计算值；\overline O 为Oj的平均值；
波浪传播稳定后从x1=80 m~280 m某一时刻的两相体积分数如图2所示。本文的数值波浪符合强非线性椭圆余弦波的形状特征，即波面在波峰附近变得很陡，两个波峰之间相隔很长且存在较为平坦的波谷。同时可以看到在选定合适的阻尼消波参数后，波浪在计算域出口端有明显的衰减，出口处波面平稳，说明末端波能消散良好。DOHMEN-JANSSEN等[4]文中水相相关结果仅给出了周期流速分布，并未详细给出其余变量的结果。KIM等[17]介绍了包含表面波效应之后的流速、紊动能和泥沙浓度分布等结果，因此后面将本模型计算的其余结果与KIM等[17]的两相模型结果对比。
为减少网格密度对近底泥沙浓度结果的影响，计算之前进行了网格无关性验证。对三套疏密不同的网格进行了计算，其中网格变动主要体现在更改了近底区域第一层网格之上的网格节点数量，其余区域网格密度保持一致，三套网格数量分别为170 000、232 500、268 750。针对图3中论述的水平流速周期分布、波峰时刻浓度剖面、波谷时刻浓度剖面结果，三套网格对应的IA值分别为：第一套网格IA=0.947、0.901、0.925；第二套网格IA=0.986、0.943、0.962；第三套网格IA=0.989、0.950、0.971。综合考虑计算效率和计算精度，最终选定计算域网格总数为232 500的第二套网格。数值计算采用的计算机CPU配置为32核。数值计算的时间步长使用T/1000，本模型计算花费约为6.7 h，使用两相模型计算花费约为21.9 h。
x1=109 m、x2=0.109 m处的水平流速周期分布如图3(a)所示，入口处的理论值“ISOBE”也放入比较。本模型很好地捕捉了椭圆余弦波水平流速的特点，即正半周期历时短、波峰陡、流速大，负半周期历时长、波谷平缓、流速小。模型计算结果相较于实验存在微小的差别(IA=0.986)，即正半周期历时更短、负半周期历时更长。出现差别的原因可能是本文用非线性更强的二阶椭圆余弦波，相比一阶解正半周期更短、负半周期更长。模拟的椭圆余弦波的非线性很强，入口处的理论值相比实验值已经有一些相位偏差。ISOBE等[34]在提出二阶椭圆余弦波理论时也探讨了不同阶数椭圆余弦波水平流速分布规律，不同阶数椭圆余弦波在同等参数下波面差异不大，但是流速存在较大差异。波浪流速相位超前会使加速阶段流速稍微大于减速阶段。此外在波浪传播过程中存在能量损耗和一定的数值振荡，容易给流速分布带来大小、相位上的微小差异。
图3(b)~图3(c)展示了波峰时刻和波谷时刻的泥沙浓度剖面。实验数据只有三个点，x2=0.9 mm、2.2 mm、4.2 mm；点与点之间用直线连接。x2=0 mm处的数据超出图标未显示。对于三个点的数据，本文数值计算结果更契合实验值(波峰时刻IA=0.943、波谷时刻IA=0.962)。本文数值计算结果和KIM等[17]两相模型结果在近底和实验值的线性连接均有差距。根据传统的动床面泥沙浓度分布理论[23]，在靠近床面的地方合理的浓度变化不应该是线性的，而是接近于负指数或幂函数，如KIM等[17]和本文的计算结果。x2=0.9 mm和x2=0 mm处实验数据的线性连接与负指数或幂函数变化产生了很大的误差。模型结果表明：泥沙量几乎与波浪近底边界层(WBBL)以上的水平方向流速大小直接相关，与物理实验中观察到的一致[4]。
CHEN等[23-24]提出的近底泥沙参考浓度和ZF94公式的对比如图4(a)所示。CHEN等[23-24]的公式受质量守恒的限制不会无限增大，所以参考浓度最大值落后于ZF94公式。由于考虑了泥沙运动的相位残留，CHEN等[23-24]公式在流速转向时刻不为0。ZF94公式仅和希尔兹数有关，增长不受质量守恒限制；与外边界层流速同相位，θ=0时浓度减为0，对于波浪输沙是不合理的。因为泥沙运动存在相位残留，即在流动转向时仍有大量泥沙运动。CHEN等[23-24]的参考浓度公式成功捕捉到了上述现象。需要额外指出的是传统OWT数值模型计算参考浓度时，切应力模化与外边界层流速相关，参数化床面剪应力的准稳态方法无法捕捉相对于波浪近底边界层上方流速的相位超前。WBBL中的这种不稳定影响可能进一步导致中粗沙向岸运输预测不足[6, 38]。本模型使用数值计算近底切应力，近底流速相较于外边界层流速存在相位超前[38-39]。因此，利用ZF94公式会超前外边界层流速一个相位；CHEN等[23-24]的公式因考虑了浓度相位漂移，相位会合理滞后。但本文模拟工况的相位超前大小和相位滞后相差不大，相位滞后参数大约为8.6°，相位超前值大约为7.9°，这导致结果与图3(a)展示的外边界层水平流速相位仍然相近。
![]()
图 3 水平流速周期分布及浓度剖面Figure 3. Horizontal velocity periodic distribution and concentration profile不同高度处的泥沙浓度周期分布结果如图4(b)~图4(d)所示。对应于图3(a)波形的强非线性，向岸阶段(t/T=0~0.4)的浓度值大于离岸阶段(t/T=0.4~1)。随着位置的升高，浓度幅值逐渐减小，KIM等[17]两相模型结果减小得比本模型更快。本模型很好的捕捉了初始床面(x2>0)以上悬沙层中的泥沙浓度分布，模型结果与实验结果在正半周期匹配较好，在负半周期本模型计算幅值偏大。实验捕捉到了泥沙浓度的相位滞后，而本模型考虑的相位滞后却和近底切应力存在的相位超前存在一定抵消，因此模型浓度分布存在一定相位超前。模型结果在x2=1 mm、2.3 mm、4.2 mm处的IA值分别为0.784、0.762、0.837，KIM等[17]两相模型结果在对应位置处的IA值分别为0.745、0.666、0.614。整体来看，本模型在悬沙层的浓度分布结果优于两相模型。
沙床上方0.2 m内的周期平均浓度剖面如图5所示。本模型计算结果在上层区域存在一定的浓度垂向梯度，与传统的数值结果一致[21]，而实验结果却几乎没有浓度梯度。DOHMEN-JANSSEN等[4]表示上层区域浓度垂向梯度很小的现象确实是令人疑惑。其观点是可能在测量过程中存在细颗粒沉积物或者是波浪边界层上方出现了额外的紊流。KIM等[17] 两相模型结果与实验相比，仅在近底层移输沙层底部区域较好，悬沙区域较差。其采用AMOUDRY等[40]的研究更改Sc重新模拟，发现依旧不能匹配用超声回馈系统测得的悬沙数据。由于粗粒径泥沙颗粒的尾流扰动更强，导致粗粒径颗粒速度脉动强度的增量大于细粒径颗粒脉动强度的增量，尾流促进颗粒扩散的效应会更明显且更符合实际情况[20]。因此对于本文计算的中沙工况来说，泥沙扩散系数中考虑了颗粒尾流涡效应[20-21]，本模型的结果相较于实验结果偏大。
2.3 孤立表面波效应
本模型计算的流速剖面同两相模型结果对比如图6所示，< >代表周期平均。本模型的波峰时刻流速剖面在近底更偏向向岸，波谷时刻流速剖面几乎吻合，整体匹配度是非常理想的。本模型的周期平均流速在近底更偏向离岸，这是因为波形效应产生了更大的离岸净流速，抵消掉更多由表面波效应引起的向岸净流动，最终导致周期平均流速结果更偏向离岸。表面波效应引起的净流动究其本质是因为波浪水槽中的垂向速度分量与水平速度分量存在同相位的成分，这必然会引起一个二阶的水平向岸恒定水流，进而增大向岸的泥沙净输运[15-16]。而往复流水槽(OWT)是一种边界层泥沙输送模型，除不能包括表面波效应外，往复流水槽在研究泥沙输运问题上具有和波浪水槽相同的能力。因此预计近底水平流速剖面和泥沙输运特性的差异仅是由于表面波效应的存在。本文参考KIM等[17]的思路，使用相同模型参数建立往复流水槽，以此来分离表面波效应产生的净流动。图6(a)~图6(b)显示在波峰、波谷时刻波浪水槽计算的流速结果相较于OWT均会更偏向向岸，很明显波浪水槽存在额外的向岸流动。图6(c)显示波浪水槽中的离岸净流动明显弱于OWT。两种水槽的差异“计算值-OWT”表示表面波效应引起的净流动，该流动是正值，会带来明显的向岸净输沙[4]。
当水流强度很大时，泥沙大量处于悬浮状态，尽管表面波效应引起的净流动值相较于最大流速而言很小，但在近底高浓度输沙层内依旧会对泥沙净输运造成很大影响。波浪水槽和往复流水槽在波峰、波谷时刻的瞬时泥沙通量剖面如图7(a)~图7(b)所示，可以明显看到波浪水槽结果受表面波效应的向岸净流动影响，结果在近底更偏向向岸输运。周期平均泥沙通量剖面如图7(c)所示，两种水槽计算的周期平均泥沙通量都是向岸的，但是波浪水槽下的结果明显呈现出了更大的向岸输移。由此可见在初始床面以上，表面波效应增强了泥沙的向岸净输移，这和KIM等[17]得到的结论一致。
3 结论
本文基于泥沙的输移扩散模型构建适用于浅水强非线性波的低计算成本数值输沙水槽，探讨包含自由表面的浅水强非线性波作用下的泥沙输运。主要结论如下：
(1) 采用边界造波法模拟了浅水强非线性二阶椭圆余弦波，得到符合预期的波面和流速。波面在波峰附近很陡，两个波峰之间相隔很长且存在较为平坦的波谷；水平流速正半周期历时短、波峰陡、流速大，负半周期历时长、波谷平缓、流速小。
(2) 数值波浪水槽同时考虑了波形效应和表面波效应，水相结果的可靠性得到两相流模型结果的验证。通过往复流水槽计算的流速结果，分离出表面波效应对泥沙输运的影响。结果显示无论波峰波谷时刻的瞬时泥沙通量还是周期平均泥沙通量，波浪水槽均具有比往复流水槽更大的近底向岸输移。
(3) 计算的泥沙浓度优于两相模型的结果，整体分布趋势也较为理想，可反映近底泥沙浓度相对流速的相位滞后。考虑的颗粒绕流涡作用促进了泥沙悬浮，计算的周期平均泥沙浓度大于实验值，在上层区域存在合理的浓度梯度。相比两相模型，本文采用更少的网格和方程达到了较好的结果，在保证精度的同时大大的节省了计算量，整体结果是可接受的。
-
![]()
图 1 浅水强非线性波作用下泥沙输运示意图
Figure 1. Sketch of sediment transport in strongly nonlinear waves
![]()
图 3 水平流速周期分布及浓度剖面
Figure 3. Horizontal velocity periodic distribution and concentration profile
-
[1] TAN W K, YUAN J. Net sheet-flow sediment transport rate: Additivity of wave propagation and nonlinear waveshape effects [J]. Continental Shelf Research, 2022, 240: 104724. doi: 10.1016/j.csr.2022.104724
[2] STACHURSKA B, STAROSZCZYK R. Laboratory study of suspended sediment dynamics over a mildly sloping sandy seabed [J]. Oceanologia, 2019, 61(3): 350 − 367. doi: 10.1016/j.oceano.2019.01.006
[3] KRANENBURG W M, RIBBERINK J S, SCHRETLEN J J L M, et al. Sand transport beneath waves: The role of progressive wave streaming and other free surface effects [J]. Journal of Geophysical Research: Earth Surface, 2013, 118(1): 122 − 139. doi: 10.1029/2012JF002427
[4] DOHMEN-JANSSEN C M, HANES D M. Sheet flow dynamics under monochromatic nonbreaking waves [J]. Journal of Geophysical Research:Oceans, 2002, 107(C10): 1 − 21. doi: 10.1029/2001JC001045
[5] RIBBERINK J S, AL-SALEM A A. Sheet flow and suspension of sand in oscillatory boundary layers [J]. Coastal Engineering, 1995, 25(3/4): 205 − 225.
[6] O'DONOGHUE T, WRIGHT S. Concentrations in oscillatory sheet flow for well sorted and graded sands [J]. Coastal Engineering, 2004, 50(3): 117 − 138. doi: 10.1016/j.coastaleng.2003.09.004
[7] WATANABE A, SATO S. A sheet-flow transport rate formula for asymmetric, forward-leaning waves and currents [M]// SMITH J M. Coastal Engineering 2004. Singapore: World Scientific Publishing, 2005: 1703 − 1714.
[8] VAN DER A D A, O’DONOGHUE T, RIBBERINK J S. Measurements of sheet flow transport in acceleration-skewed oscillatory flow and comparison with practical formulations [J]. Coastal Engineering, 2010, 57(3): 331 − 342. doi: 10.1016/j.coastaleng.2009.11.006
[9] RUESSINK B G, MICHALLET H, ABREU T, et al. Observations of velocities, sand concentrations, and fluxes under velocity-asymmetric oscillatory flows [J]. Journal of Geophysical Research: Oceans, 2011, 116(C3): C03004.
[10] LONGUET-HIGGINS M S. Mass transport in water waves [J]. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 1953, 245(903): 535 − 581.
[11] SCHRETLEN J L M. Sand transport under full-scale progressive surface waves [D]. Enschede: University of Twente, 2012.
[12] HASSAN W N M, RIBBERINK J S. Modelling of sand transport under wave-generated sheet flows with a RANS diffusion model [J]. Coastal Engineering, 2010, 57(1): 19 − 29. doi: 10.1016/j.coastaleng.2009.08.009
[13] HOLMEDAL L E, MYRHAUG D. Wave-induced steady streaming, mass transport and net sediment transport in rough turbulent ocean bottom boundary layers [J]. Continental Shelf Research, 2009, 29(7): 911 − 926. doi: 10.1016/j.csr.2009.01.012
[14] FUHRMAN D R, SCHLØER S, STERNER J. RANS-based simulation of turbulent wave boundary layer and sheet-flow sediment transport processes [J]. Coastal Engineering, 2013, 73: 151 − 166. doi: 10.1016/j.coastaleng.2012.11.001
[15] KRANENBURG W M, RIBBERINK J S, UITTENBOGAARD R E, et al. Net currents in the wave bottom boundary layer: On waveshape streaming and progressive wave streaming [J]. Journal of Geophysical Research:Earth Surface, 2012, 117(F3): F03005.
[16] TAN W K, YUAN J. A two-layer numerical model for coastal sheet-flow sediment transport [J]. Journal of Geophysical Research:Oceans, 2021, 126(6): e2021JC017241.
[17] KIM Y, CHENG Z, HSU T J, et al. A numerical study of sheet flow under monochromatic nonbreaking waves using a free surface resolving Eulerian two-phase flow model [J]. Journal of Geophysical Research:Oceans, 2018, 123(7): 4693 − 4719. doi: 10.1029/2018JC013930
[18] 刘传林, 韩新宇, 董胜. 半椭圆型防波堤水动力特性的数值模拟[J]. 工程力学, 2023, 40(10): 247 − 256. doi: 10.6052/j.issn.1000-4750.2022.01.0094 LIU Chuanlin, HAN Xinyu, DONG Sheng. Numerical simulation of hydrodynamic characteristics of semi-elliptic breakwater [J]. Engineering Mechanics, 2023, 40(10): 247 − 256. (in Chinese) doi: 10.6052/j.issn.1000-4750.2022.01.0094
[19] 张治国, 叶铜, 张成平, 等. 基于Pasternak海床模型的椭圆余弦波浪荷载作用下埋置管线动力响应解析解[J]. 工程力学, 2024, 41(1): 76 − 89. doi: 10.6052/j.issn.1000-4750.2022.02.0190 ZHANG Zhiguo, YE Tong, ZHANG Chengping, et al. Analytical solution for dynamic response of buried pipeline under cnoidal wave load based on Pasternak seabed model [J]. Engineering Mechanics, 2024, 41(1): 76 − 89. (in Chinese) doi: 10.6052/j.issn.1000-4750.2022.02.0190
[20] SHI H B, YU X P. An effective Euler-Lagrange model for suspended sediment transport by open channel flows [J]. International Journal of Sediment Research, 2015, 30(4): 361 − 370. doi: 10.1016/j.ijsrc.2015.03.012
[21] LIANG L X, YU X P, BOMBARDELLI F. A general mixture model for sediment laden flows [J]. Advances in Water Resources, 2017, 107: 108 − 125. doi: 10.1016/j.advwatres.2017.06.012
[22] 何康, 施华斌, 余锡平. 基于两相流理论的稀疏和致密颗粒流统一模型[J]. 工程力学, 2023, 40(8): 24 − 35. doi: 10.6052/j.issn.1000-4750.2021.12.0978 HE Kang, SHI Huabin, YU Xiping. The unified model for dilute and dense granular flows based on the two-phase fow theory [J]. Engineering Mechanics, 2023, 40(8): 24 − 35. (in Chinese) doi: 10.6052/j.issn.1000-4750.2021.12.0978
[23] CHEN X, NIU X J, YU X P. Near-bed sediment condition in oscillatory sheet flows [J]. Journal of Waterway, Port, Coastal, and Ocean Engineering, 2013, 139(5): 393 − 403. doi: 10.1061/(ASCE)WW.1943-5460.0000193
[24] CHEN X, LI Y, WANG F J. Mobile bed thickness in skewed asymmetric oscillatory sheet flows [J]. Acta Mechanica Sinica, 2018, 34(2): 257 − 265. doi: 10.1007/s10409-017-0686-3
[25] RIJN L C V. Sediment transport, part II: Suspended load transport [J]. Journal of Hydraulic Engineering, 1984, 110(11): 1613 − 1641. doi: 10.1061/(ASCE)0733-9429(1984)110:11(1613)
[26] LIU H J, SATO S. Modeling sediment movement under sheetflow conditions using a two-phase flow approach [J]. Coastal Engineering Journal, 2005, 47(4): 255 − 284. doi: 10.1142/S0578563405001276
[27] BOMBARDELLI F A, JHA S K. Hierarchical modeling of the dilute transport of suspended sediment in open channels [J]. Environmental Fluid Mechanics, 2009, 9(2): 207 − 235. doi: 10.1007/s10652-008-9091-6
[28] 程友良, 许强, 焦慎俐. 速度边界造波法与推板造波法的差异化分析[J]. 海洋技术学报, 2021, 40(5): 105 − 116. doi: 10.3969/j.issn.1003-2029.2021.05.013 CHENG Youliang, XU Qiang, JIAO Shenli. Differentiation analysis of velocity boundary wave making method and push plate wave making method [J]. Journal of Ocean Technology, 2021, 40(5): 105 − 116. (in Chinese) doi: 10.3969/j.issn.1003-2029.2021.05.013
[29] 陈林雅, 郑东生, 王盼娣, 等. 基于波浪-多孔介质海床-结构物耦合模型的单桩基础波浪力分析[J]. 工程力学, 2019, 36(11): 72 − 82. doi: 10.6052/j.issn.1000-4750.2018.08.0455 CHEN Linya, ZHENG Dongsheng, WANG Pandi, et al. Analysis of wave force acting on the monopile based on wave-porous seabed-structure coupled model [J]. Engineering Mechanics, 2019, 36(11): 72 − 82. (in Chinese) doi: 10.6052/j.issn.1000-4750.2018.08.0455
[30] 李邦华, 郑向远, 李炜, 等. 波浪水槽中Stokes五阶波的数值生成[J]. 武汉理工大学学报(交通科学与工程版), 2016, 40(2): 238 − 244, 250. LI Banghua, ZHENG Xiangyuan, LI Wei, et al. Simulation of fifth-order stokes waves in the numerical wave flume [J]. Journal of Wuhan University of Technology (Transportation Science & Engineering), 2016, 40(2): 238 − 244, 250. (in Chinese)
[31] 郑艳娜, 刘卓, 陈昌平, 等. 基于Fluent软件二维数值波浪水槽的研究[J]. 中国海洋平台, 2015, 30(6): 67 − 71. doi: 10.3969/j.issn.1001-4500.2015.06.012 ZHENG Yanna, LIU Zhuo, CHEN Changping, et al. Study on the two-dimension numerical wave tank based on fluent [J]. China Offshore Platform, 2015, 30(6): 67 − 71. (in Chinese) doi: 10.3969/j.issn.1001-4500.2015.06.012
[32] 刘莎莎, 顾煜炯, 惠万馨, 等. 基于边界造波法的波浪数值模拟[J]. 可再生能源, 2013, 31(2): 100 − 103. doi: 10.13941/j.cnki.21-1469/tk.2013.02.007 LIU Shasha, GU Yujiong, HUI Wanxin, et al. Wave numerical simulation based on wave-generation method of defining inlet boundary conditions [J]. Renewable Energy Resources, 2013, 31(2): 100 − 103. (in Chinese) doi: 10.13941/j.cnki.21-1469/tk.2013.02.007
[33] LAITONE E V. The second approximation to cnoidal and solitary waves [J]. Journal of Fluid Mechanics, 1960, 9(3): 430 − 444. doi: 10.1017/S0022112060001201
[34] ISOBE M, KRAUS N C. Derivation of a second-order cnoidal wave theory [M]. Yokohama, Japan: Yokohama National University, Department of Civil Engineering, Hydraulic Laboratory, 1983: 23 − 27. [35] ZYSERMAN J A, FREDSØE J. Data analysis of bed concentration of suspended sediment [J]. Journal of Hydraulic Engineering, 1994, 120(9): 1021 − 1042. doi: 10.1061/(ASCE)0733-9429(1994)120:9(1021)
[36] WILLMOTT C J, WICKS D E. An empirical method for the spatial interpolation of monthly precipitation within California [J]. Physical Geography, 1980, 1(1): 59 − 73. doi: 10.1080/02723646.1980.10642189
[37] WILSON K C. Analysis of bed-load motion at high shear stress [J]. Journal of Hydraulic Engineering, 1987, 113(1): 97 − 103. doi: 10.1061/(ASCE)0733-9429(1987)113:1(97)
[38] NIELSEN P. Sheet flow sediment transport under waves with acceleration skewness and boundary layer streaming [J]. Coastal Engineering, 2006, 53(9): 749 − 758. doi: 10.1016/j.coastaleng.2006.03.006
[39] NIELSEN P. Shear stress and sediment transport calculations for swash zone modelling [J]. Coastal Engineering, 2002, 45(1): 53 − 60. doi: 10.1016/S0378-3839(01)00036-9
[40] AMOUDRY L, HSU T J, LIU P L F. Schmidt number and near-bed boundary condition effects on a two-phase dilute sediment transport model [J]. Journal of Geophysical Research:Oceans, 2005, 110(C9): C09003.
下载:
计量
- 文章访问数: 173
- HTML全文浏览量: 34
- PDF下载量: 53
