基于海床孔压响应的泥沙起动机制研究

国振

沈侃敏

王立忠

郑东生

1 引言

泥沙起动条件是泥沙输运理论的基础，一般从泥沙颗粒的受力平衡出发，推导泥沙由静止转化为运动状态的临界力学条件，如起动剪应力或流速。

对于铺设在海床上的油气管道而言，由于底床崎岖不平、工作热应力等原因易产生局部悬跨段。海水在流经悬跨段时，往往伴随着漩涡的周期性发放，同时也受到底床摩擦的影响，其流态比较复杂。已有研究大多针对较为简单的水流动力条件，如单向流、振荡流等，没有考虑此类复杂流动条件对泥沙颗粒起动的影响，仍需从力学机理层次对其进一步探索。

首先，应选用合适的湍流模型，以便于准确捕捉悬跨段附近流场的分布及变化特征。湍流场的数值模拟可分为直接模拟（Direct Numerical Simulation，简称DNS），大涡模拟（Large Eddy Simulation，简称LES）和Reynolds时均方程（Reynolds Averaged Navier-Stokes，简称RANS）三类。其中，RANS方法对计算能力需求相对较低，仍是目前工程研究中最现实可行的方法，也为各主流商业软件（Fluent、CFX及Comsol等）所采用。Brørs利用k—ε模型和壁函数模拟了海底管道周围的流场及漩涡脱落情况；Liang和Cheng对比分析了k—ε模型，Wilcox高、低雷诺数k—ω模型和SGS模型，认为采用无滑移壁面边界的低雷诺数k—ω模型可以更好地预测管道周围的漩涡脱落情况。

此外，海床土体往往是包含固-液两相的多孔介质，管道悬跨段附近的流场压力直接影响海床变形与内部孔压场的分布，这对泥沙颗粒起动的影响不可忽视。

因此，本文从计算模拟悬跨段管道附近的湍流场入手，分析获得底床表面上的压强分布规律与时域内的波动特性，并探讨了平均流速U₀和间隙比e/D（e为管道底部与海床距离，D为管道直径）的影响；然后，将海床表面上的波动流场压强作为底床的渗流边界条件，在时域内顺次耦合求解获得海床内部的孔压场变化，深入分析了悬跨段附近底床各区域的泥沙起动机制。本文从多孔介质力学和土力学角度完善了对复杂流动下泥沙颗粒起动机制的认识。

2 数值模型

2.1 流体控制方程

本文的湍流模型采用Menter等提出的SST（剪应力输运）模型。该模型综合了Wilcox k—ω模型在近壁区的求解精度和标准k—ε模型在湍流核心区的数值鲁棒性，其k和ω方程可分别表述为：

式中 ρ_w——水体单位质量；

μ——水的动力黏度系数。

湍流动能产生率P及P_k、漩涡黏度v_t和应变率模量S计算如下：

式中，为水质点的速度张量；参数a₁=0.31，；其他参数采用φ=F₁φ₁+（1-F₁）φ₂形式进行插值，φ表征参数（β，γ，σ_k，σ_ω）。第一组参数取β₁=0.075，γ₁=0.556，σ_k1=0.85，σ_ω1=0.5；第二组参数取β₂=0.0828，γ₂=0.44，σ_k2=1.0，σ_ω2=0.856。

两个混合函数F₁和F₂分别为：

式中 l_w——计算节点距离壁面的最近距离。

2.2 底床控制方程

海水在流经悬跨段时，其雷诺数Re通常处于亚临界范围（150<Re<3×10⁵），此时漩涡将以一个较窄频带周期性释放，漩涡脱落频率f_s=SU₀/D，其中S为Strouhal数。对于大多数实际的油气管道而言，其Strouhal数可近似取为0.2。试验研究发现，在管道下游尾流区内，因受漩涡周期性脱落影响，底床附近的流场速度也呈周期性变化，表现出“波动”特性。

因此，本文借鉴了波浪—海床作用理论的研究思路，结合Biot-Darcy方程和Seed残余孔压源项的二维格式，来模拟在海床表面流场压强周期性作用下其内部孔压的振荡响应和累积增长趋势，这分别表征了在循环动载下海床土体骨架的弹性和塑性体变。

2.2.1 振荡孔压响应

假定海床土体为均匀、各向同性的多孔弹性介质，其平面应变条件下的控制方程为：

式中，——在水平x和竖直z轴方向的土体位移；——流场动压引起的振荡孔压；——土骨架弹性体变；

G_s——土体剪切模量；

μ_s——泊松比；

k——渗透系数；

n_s——孔隙比；

γ_w——孔隙水体单位容重；

β_s——孔隙水的压缩系数；

K_w——纯水的体积模量，K_w=2×10⁹N/m²；

P_w0——计算点静水压力与大气压力之和；

S_r——土体饱和度。

2.2.2 残余孔压计算

基于seed等人提出的残余孔压累积计算的一维模式，本文推广建立了残余孔压累积增长的二维计算模型，其控制方程可表述如下：

式中——流场动压作用引起的残余孔压增长；

c_v2——平面应变条件下的固结系数；

f（x，z，t）为二维的残余孔压累积源项^［13］，表述为：

式中——计算点处某一时刻t的剪应力幅值，可由振荡孔压计算中获得；

T_s——该点处循环荷载的平均作用周期；

α_r，β_r——经验参数；——土体的有效正应力平均值。

对于K₀固结海床，可取为：

式中 K₀——静止土压力系数；

γ′——土体有效容重，等于（γ_s-γ_w），γ_s为海床土体的饱和容重。

2.3 计算域及边界设定

如图1所示，本文的计算模型分为两个计算区域：域1赋予流体控制方程，域2赋予底床控制方程。采用笛卡尔坐标系xoz，设定x轴正方向为水平来流方向，z轴的正方向竖直向上。计算区域水平长度为-10≤x/D≤20；域1竖向高度为0≤x/D≤6，域2竖向高度为-8≤x/D≤0，管道间隙比为e/D。

2.3.1 域1边界条件

边界AB：海床表面，为无滑移壁面边界。为保证湍流场计算精度，边界AB以上第一层网格节点的高度满足z⁺=z·u^*/μ_w<1；

边界BC：流场出口边界，设为p=0；

边界CD：设为自由滑动边界；

边界AD：流场入口边界，设为完全发展的湍流场，其水质点水平流速u（z）和竖向流速w（z）分别为：

式中 u^*——摩阻流速；

κ——卡门常数。

2.3.2 域2边界条件

边界AB：，；

边界BE与AF：，；

边界EF：，。

本数值模型的常规计算参数取值参见图1。

3 数值模拟结果与分析

泥沙颗粒的起动一般涉及两个方面因素的影响：一是近底床的水流条件；二是泥沙力学特性。因此，基于本文所建二维数值模型，首先研究了流场漩涡脱落规律及水质点速度分布，分析获得海床表面的流体压强特征，探讨了平均流速U₀（=0.1m/s，0.5m/s，1.0m/s，2.0m/s）和间隙比e/D（=0.1，0.4，0.6，1.0）的影响；基于此，研究流体压强持续作用下底床内部孔压（振荡孔压与残余孔压）的发展规律，并探讨了在底床不同区域内的泥沙起动机制。

3.1 管道附近流场分布特征

圆柱绕流是一个经典而传统的流体力学问题，对于悬跨段管道而言，还需考虑底床的近壁效应影响。试验表明，当管道间隙比e/D<0.3时，漩涡脱落现象被抑制，故选择e/D=0.4，平均流速U₀=0.5m/s。C_L=F_L/（2ρ_wDU₀）是管道上横向涡激力作用系数，如图2所示，涡激力幅值在最初几个周期后趋于稳定，其作用周期为3.0s，即反映了漩涡脱落的主导频率等于0.33，计算得Strouhal数等于0.2，侧面验证了本湍流模型的可靠性。

图3为一个典型周期内（t₀→t₀+3×T_s/4）的流场水质点速度变化情况。由图中可见，漩涡从管道的顶部和底部交替脱落，水质点与底床存在明显的相互作用。除了管道底部水流速度较大之外，在下游底床附近的水质点速度也存在明显的波动，这在以往研究中被认为是尾流区冲刷坑发展的主要原因^［10］。然而，这些研究仅从近底床水质点切向流速或切向力角度来评估底床的泥沙起动条件，忽视了尾流区的复杂流场（切向、垂向流速）带来的海床表面水压的动态变化，也无法反映海床内部的多孔介质渗流特性，没有从本质上反映底床泥沙的起动机制。

因此，本文在计算中选择了在海床表面均匀布置的8个测点（图4），以此来分析管道附近流场影响区域内表面压强的动态变化特征。

图5为截取计算时间150s内8个测点处的水流压强时程变化曲线。如图所示，所有测点的压强波动幅值和周期在约60s后逐渐趋于稳定；由上游经管道底部至下游，各测点处的压强波动幅值先增大，再减小；各测点的压强波动周期几乎一致，大约等于漩涡的脱落周期3.0s，表现出很好的周期性和规律性。

为了更好地总结规律，本文绘制了流场稳定之后200s时间内的海床表面压强分布，如图6所示：①上游床面压强始终为正压，波动性较小；②管道底部一定范围内床面的压强为负值，并表现出强烈的波动性；③在1D～4D范围内的床面压强变为正压，波动较明显；④至下游6D之后，床面压强的波动可忽略不计。本文的湍流模型很好地描述了底床上平均压强的分布规律及其波动情况。可以发现，悬跨段管道漩涡脱落现象的主要影响区域大约在管道底部及下游尾流区4D范围内；底床不同区域的泥沙经历了压强水平和波动幅度各异的动水压力作用。

针对漩涡脱落的主要影响区0≤x≤4D，本文计算获得了不同平均流速U₀和间隙比e/D对x/D=0.0，2.0，4.0三个测点处压强的影响。图7、图8分别为平均流速U₀=0.1m/s，0.5m/s，1.0m/s，2.0m/s和间隙比e/D=0.1，0.4，0.6，1.0的数值计算结果。

由图7可见，对于间隙比e/D=0.4的情况，随着平均流速U₀的逐渐增加，三个测点处的压强波动幅值显著增大，其波动周期逐渐减小。如图8所示，对于平均流速U₀=0.5m/s的情况，当间隙比e/D=0.1时，底床的近壁效应使得漩涡脱落现象被完全抑制，各测点压强几乎保持为定值；当e/D达到0.4后，各测点压强波动明显；随着间隙比的进一步增大，因管道远离海床，各测点受到的流场压强波动幅值逐渐降低，但是作用周期几乎不变。

3.2 海床内部孔压发展

传统的泥沙起动研究针对于突出床面上的单个颗粒，考虑其受到流体拖曳力和颗粒上下压差造成的上举力的平衡。对于平坦海床而言，大多数颗粒受其周围颗粒的嵌固作用，拖曳力的影响往往较小，主要是上举力使其脱离床面^［1］，而上举力更多地来自于海床内部渗流作用。因此，本文依次对床面流体压强循环作用下海床内部的振荡孔压分布和残余孔压发展情况进行分析。在本算例中，平均流速U₀=2.0m/s，间隙比e/D=0.4，其他计算参数见图1。

3.2.1 振荡孔压分布

海床中的振荡孔压产生与消散和土体骨架的弹性体积变形密切相关。因此，仅针对计算时程中一个周期内的振荡孔压场分布情况展开分析，如图9所示。伴随着上部流场中漩涡的产生与脱落，海床内部振荡孔压由管道底部向其下游区域逐渐传播，振荡幅值递减，这与图6所示海床表面流体压强分布特征相一致。

这里采用Jeng提出的瞬时液化准则（20）来判定海床瞬时液化深度的发展，即满足：

图10表征了在海床三个位置处，即x/D=0.0，1.0，2.0，振荡孔压幅值沿深度方向的分布。结合图中所示土体有效正应力均值曲线，可以发现：管道底部区域的瞬时液化深度约0.1D=0.03m；在x=1D处其液化深度显著减小，仅为0.05D=0.005m；而在x=2D之后的底床不发生瞬时液化。振荡孔压的影响区域基本集中在管道底部至下游1D范围以内的海床。

3.2.2 残余孔压累积

作为砂性海床液化的另一种机制，残余液化是流体压强循环作用下土体骨架逐渐产生塑性体变的结果，这导致孔隙水压力随时间不断累积，直至达到该点处的土体有效正应力平均值，此时残余液化发生。如图11所示，在200s计算时间内海床内残余孔压累积主要发生在管道下游1D～2D范围内。

同样，针对在海床三个位置（x/D=0.0，1.0，2.0）残余孔压沿深度方向随时间的发展进行分析，如图12所示：管道底部完全不会发生残余液化；在x=1D处海床的残余液化发展相对较快，200s时其液化深度已超过0.3D=0.09m，远大于该处的瞬时液化深度；x=2D处在大约100s后才开始产生残余液化，200s时其液化深度接近0.15D=0.045m。

综上所述，由于悬跨段管道附近流场的分布特征与波动特性，海床在不同区域的泥沙起动机制也将有所区别：在管道底部至下游1D范围以内，底床泥沙起动的上举力主要来自于振荡孔压；在管道下游1D～2D范围内，由压强波动持续作用产生的残余孔压对泥沙起动起主要贡献，这从多孔介质力学角度解释了Sumer等试验观测到的尾流区冲刷坑现象。

4 结论与建议

本文结合SST湍流模型，Biot-Darcy方程和Seed孔压累积源项构建了考虑复杂流场作用下的泥沙起动机制研究多场耦合模型，深入研究了悬跨段管道周围复杂流场的分布特征与波动特性，分析了管道底床不同区域泥沙的起动机制及规律，主要得到以下结论：

（1）管道底部及下游海床表面水流压强受漩涡脱落影响，呈较为规律的波动性；管道底部海床受向上的负压作用，压强波动强烈，而远离管道的下游方向波动性逐渐减弱。

（2）随平均流速U₀增加，海床表面的压强波动幅值显著增大，其周期逐渐降低；间隙比e/D=0.1时，漩涡脱落现象被完全抑制；当e/D超过0.4并逐渐增大时，海床表面的压强波动幅值明显降低，但是波动周期保持不变。

（3）在悬跨段附近流场作用下，管道底部1D范围内海床的孔压振荡明显，此区域泥沙起动主要受振荡孔压影响；管道下游1D～2D范围内海床内部的残余孔压随时间累积显著，此部分孔压对该区域内的泥沙起动起主要贡献。

参考文献

［1］钱宁，万兆惠.泥沙运动力学［M］.北京：科学出版社，2003.

[2]Brørs B.Numerical modeling of flow and scour at pipelines[J].Journal of Hydraulic Engineering，ASCE，1999，125（5）：p.511-523.

［3］Liang D F，Cheng L.Numerical modeling of flow and scour below a pipeline in currents Part I.Flow simulation［J］.Coastal Engineering，2005，52：25-42.

［4］Cheng Nian-sheng，Chiew Yee-meng.Incipient sediment motion with upward seepage［J］.Journal of Hydraulic Research，1999，37（5）：665-681.

［5］Zhou Chun yan，Li Guang xue，DONG Ping，SHI Jing-hao，XU Jishang.An experimental study of seabed responses around a marine pipeline under wave and current conditions［J］.Ocean Engineering，2011，38：226-234.

［6］Menter F R.Two-equation eddy-viscosity turbulence models for engineering applications［J］.AIAA Journal，1994，32（8）：1598-1605.

［7］Menter F R，Kuntz M，Langtry R.Ten years of industrial experience with the SST turbulence model［J］.Turbulence Heat and Mass Transfer 4，2003.

［8］吴钰骅.海底管道—流体—海床相互作用机理和监测技术研究［D］.杭州：浙江大学，2007.

［9］Lei C，Cheng L，Kavanagh K.Reexamination of the effect of a plane boundary on force and vortex shedding of a circular cylinder［J］.Journal of Wind Engineering and Industrial Aerodynamics，1999，80（3）：263-286.

[10]Sumer B M，Jensen H R，Mao Y，Fredsøe J.Effect of lee-wake on scour below pipelines in current[J].Journal of Waterway，Port，Coastal and Ocean Engineering，1988，114：599-614.

［11］Biot M A.Theory of deformation of a porous viscoelastic anisotropic solid［J］.Journal of Applied Physics，1956，27：459-467.

［12］Seed H B，Rahman M S.Wave-induced pore pressure in relation to ocean floor stability of cohesionless soils［J］.Marine Geotechnology，1978，3（2）：123-150.

[13]Sumer B M，Fredsøe J.The mechanics of scour in the marine environment[M].Singapore：World Scientific，2002.

［14］Jeng D S.Porous Models for Wave-seabed Interactions［M］.German：Springer，2013.

［15］Bearman P W，Zdravkovich M M.Flow around a circular cylinder near a plane boundary［J］.Journal of Fluid Mechanics，1978，89（1）：33-47.

［16］Grass A J，Raven P W J，Stuart R J，Bray J A.The influence of boundary layer velocity gradients and bed proximity on vortex shedding from free spanning pipelines［J］.Journal of Energy Resources Technology，1984，106：70-78.