孤立波作用下埋管斜坡海床及海底管道的响应分析

潘佳禾， 廖晨聪,b,c， 陈锦剑,b,c

(上海交通大学 a. 土木工程系； b. 海洋工程国家重点实验室； c. 高新船舶与深海开发装备协同创新中心， 上海 200240)

随着人类对海洋的开发，近岸海域内工程结构物的数量日益增加.以斜坡为主的近岸海床是海底油气管道“登陆”的重要过渡地带.在近岸浅水环境下，海床及管道的受力状态易受波浪作用的影响，浅水区域的波浪呈现出诸如孤立波的特殊形态.孤立波在斜坡海床上爬升和回落时，会引起孔隙水压力变化，可能导致海床的液化与滑坡，进而危及海底管道等近岸结构物.

孤立波是一种模拟近岸浅水环境的常用波形，也被广泛地用于模拟海啸等极端波浪[1].当孤立波传播到近岸斜坡海底地带时，由于水深变浅，波浪的波幅和非线性均有所增加，而到达海岸线后波浪还会继续向陆地爬升及回落.对于孤立波在斜坡上的破碎、爬升及回落过程已有较多实验观测和数值模拟研究.Synolakis[2]通过水槽实验捕捉斜坡上孤立波的变化波面.宣瑞韬[3]采用考虑造波区波动非稳态特性的造波方法，在水槽实验中模拟孤立波和双孤立波波形的传播及爬高过程.王贺等[4]基于雷诺平均Navier-Stokes(RANS)方程、流体体积(VOF)方法及修正的Goring造波方法，实现对孤立波在斜坡上的破碎、爬升及回落过程的数值模拟.实际上，孤立波在近岸斜坡上的这种破碎、爬升及回落过程将会影响海床受力及孔压响应.

关于波浪、海床与结构物之间的相互作用问题，已有不少学者开展了研究工作.刘博等[5]给出一个近似u-p(固相位移-孔压)动力模型的解析解，考虑在波流共同作用下多孔介质海床的孔隙水压力及有效应力的变化，并研究了潮流流向对海床响应的影响.Jeng等[6]通过Navier-Stokes(N-S)方程控制波浪，提出一个2维波浪-海床-防波堤之间相互作用的单向耦合模型，较精准地模拟了三者之间的相互作用.Zhang等[7]在此基础上进一步提出单桩-波浪-海床之间相互作用的3维模型，研究海床土体的动力响应问题.张军等[8]利用COMSOL软件构建了波浪-海床-管线计算模型，模拟了在1阶Stokes波作用下管线周围土体孔压及其有效应力的分布情况.胡翔等[9]采用准静态3维数值分析方法，研究了波浪载荷作用下海床的孔压变化规律及饱和砂质海床中的单桩响应问题.陈宝清等[10]采用计算流体力学的OpenFOAM软件，建立了波浪-海床-结构物数值模型，模拟了在结构物影响下的海床响应.Jeng等[11]提出一种新的数值模型，分别探讨了波浪作用下各项同性与各项异性、均质与非均质海床中海底管道附近的孔压、有效应力以及管道自身的应力响应，并分析了不同波浪参数、海床土体参数、管道尺寸对计算结果的影响.

上述研究主要集中于水平海床(相当于远海环境)在周期性波浪作用下的响应，对孤立波的研究主要聚焦在模拟水动力学问题的破碎、爬升及回落过程.最近，有学者开展了对于近岸斜坡海床及结构物响应的研究.Zhao等[12]建立了2维数值模型，研究在孤立波作用下1∶6和1∶15两种倾斜度多孔海床的累积孔隙水压对液化范围的影响.Zhao等[13]采用数值方法研究在波浪作用下防波堤附近倾斜海床的土体响应，针对海床渗透系数、波浪周期等参数进行分析，并对比了防波堤附近不同位置的孔压.Young等[14]在48 m长的水槽中模拟了海啸波(破碎孤立波)对近岸沙滩的侵蚀和沉积作用.Xiao等[15]利用物理实验数据，对海啸作用下的沙质海床进行数值模拟，分析不同波高和海床渗透系数对液化深度的影响.Gao等[16]设计管土相互作用设备进行足尺物理试验，研究波浪作用下铺设于倾斜海床表面的海底管道上坡及下坡的稳定性.

本文针对近岸斜坡海床，建立孤立波-斜坡海床-海底管道耦合模型.采用N-S方程和k-ε湍流模型对孤立波进行数值模拟，并将得到的斜坡表面波压力作用于海底斜坡模型.对于斜坡土体，基于Biot多孔弹性固结理论建立海床控制方程；对于海底管道，根据线弹性理论建模，模拟分析埋管斜坡海床和管道在孤立波作用下的响应.通过对比数值模拟结果与试验数据及其解析结果，验证所建模型的准确性.根据计算结果，分析在孤立波作用下斜坡上波压力的水平分布随时间的变化特点及不同位置的孔压随时间的变化规律.此外，对斜坡模型进行参数分析，研究渗透系数对孔压纵向分布的影响.

1 模型与计算方法

孤立波在传播至近岸斜坡时会发生波浪的破碎，并在水平面以上的斜坡继续爬升，至最高点后回落.为了准确地模拟孤立波在斜坡上的破碎、爬升及回落过程中波压力的变化，且以此为基础研究海床的孔压响应，将模型简化为波浪、斜坡海床、海底管道3个部分进行耦合计算.

首先，通过Flow-3D软件中的水槽模拟孤立波在斜坡上的破碎、爬升及回落过程，获得斜坡表面的水压力，以水压力减去初始静止水压得到波浪压力pb；然后，采用COMSOL软件对斜坡海床和海底管道进行建模，以边界条件的方式利用插值函数将pb施加于斜坡模型的上边界(水土交界处)，计算与分析斜坡土体在孤立波作用下的孔压响应.孤立波-斜坡海床-海底管道耦合模型如图1所示，模型坐标原点O设置于水平海床底面及倾斜斜坡底面的交界处，孤立波由左侧产生，x轴为波浪行进方向，z轴为波浪高度方向.

图1 孤立波-斜坡-海底管道耦合模型Fig.1 The solitary wave-sloping seabed submarine-pipeline coupling model

1.1 孤立波模型

在孤立波数值模型中，液体的运动控制方程采用N-S方程，选取k-ε湍流模型以更真实地反映波浪与海底斜坡之间的相互作用.Flow-3D软件中用FAVOR法划分网格会产生流体面积和体积分数项，其连续性方程及动量方程表达式如下：

(1)

(2)

(3)

式中：Ax,Az分别为x和z方向流体面积分数；φF为流体体积分数；vx,vz分别为x和z方向流体速度；ρ为流体密度；p为孔隙水压力；μ为动力黏度；g为重力加速度；ax,az分别为x和z方向黏性加速度.

连续性方程和动量方程的k-ε湍流模型：

PT+fT+DkT-εT

(4)

(5)

式中：kT为湍流动能；εT为湍流动能耗散率；PT为湍流动能产生项；fT为湍流浮力产生项；DkT为湍流动能扩散项；DεT为湍流动能耗散率的扩散项；C1,C2,C3为无量纲计算系数.

孤立波模型的左侧为波浪入射边界，基于孤立波理论造波；孤立波模型的右侧为波浪流出边界(事实上，由于斜坡长度足够长，保证波浪爬坡后会完全回落，并不会从右侧流出)；孤立波模型的上部为水和空气的交界面，气压等于一个标准大气压(101.3 kPa)；孤立波模型的底部为水土交界面，采用墙面边界，边界上的流体法向速度为0.

1.2 埋管海床模型

1.2.1海床模型 采用COMSOL软件对斜坡海床建模，基于Biot固结方程，利用偏微分方程模块对海底斜坡设置本构关系方程.以边界条件的方式用插值函数将pb施加于斜坡模型的上边界(水土交界处)，计算斜坡土体在孤立波作用下的孔压响应.基于Biot多孔弹性固结理论的海床控制方程为

(6)

(7)

假设海床为均匀、饱和多孔、各向同性的弹性土体，其几何非线性可以忽略，并且土体的渗透系数k在各方向上皆相同，则土体孔隙中的液体压力可以表述为

(8)

式中：ns为土体孔隙率；γw为水的重度；βs为孔隙中液体的压缩系数，

(9)

Kw为水的体积模量(此处取2 GPa)；pw0为绝对水压力；S为饱和度.

(10)

海床底部可以视为不透水的刚性底面.该底面海床土体的水平位移及纵向位移为0，海水的法向流量也为0，因此边界条件可以表示为

us=ws=0

(11)

在纵向边界上，x方向的位移为零.

1.2.2管道模型 采用COMSOL软件对海底管道进行建模，该模型基于线弹性理论，采用与海床控制方程类似的偏微分方程控制海底管道响应，

(12)

(13)

式中：up,wp分别为海底管道在x，z方向的位移；Gp为管道材料的切变模量；υp为管道材料的泊松比.

海底管道外壁视为不透水边界，故管道外壁表面法向流量为0.假设管道外壁与海床土体之间无相对位移，即管道外壁在x、z方向的位移与海床土体在x、z方向的位移相等，则管道外壁的边界条件可以表示为

∂p/∂n=0，us=up，ws=wp

(14)

式中：n为垂直于海底管道壁切线的方向.

2 模型验证

2.1 孤立波模型验证

Synolakis[2]针对孤立波在斜坡上的破碎、爬升及回落过程，在水槽和固定斜坡上开展试验，并采用该试验数据验证孤立波模型的可靠性.设定模型参数：最大水深为0.29 m，波高水深比为0.28，斜坡倾斜度为1∶20.

孤立波模型与Synolakis试验数据[2](TS1)的对比结果如图2所示，其中t为波浪运动时间.图2(a)～(c)所示为孤立波在斜坡上的破碎及爬升过程，图2(d)所示为波浪从最高点回落的阶段.

图2 孤立波模型与TS1数值结果的对比Fig.2 Comparison of numerical results between solitary wave numerical model and TS1

由图2可见， TS1基本分布在孤立波模型模拟的波浪表面，虽然孤立波模型波峰高度的模拟结果比TS1值略大，但其波浪的爬升高度及波浪形状与TS1的结果吻合得较好.因此，孤立波模型能较好地模拟孤立波在斜坡上的破碎、爬升及回落过程.

2.2 埋管海床模型验证

为了验证埋管海床模型的准确性，将埋管海床模型的数值计算结果与文献[17]的解析结果(TS2)和文献[18]的试验数据(TS3)进行对比.设波浪参数：波高为 5.24 mm；水深d为0.533 m；周期为0.9 s；波长为1.25 m.海底管道参数：弹性模量Ep为 68 GPa；υp为0.32；管道外径R为 0.042 m；埋深e为0.659 m.斜坡海床参数：海床厚度H为0.826 m；ns为0.42；υs为0.33；侧向土压力系数为0.41；k为1.1 mm/s；S为1；土体密度ρs为 2 053 kg/m3；Gp为0.64 MPa.

在波浪作用下，管道周围土体孔压沿管道圆周分布，埋管海床模型的计算结果与TS2及TS3的数据对比结果如图3所示，其中p0为初始孔隙水压力.由图3可知：埋管海床模型与TS2在数值计算结果及变化趋势上均吻合得较好；与TS3在数值计算结果上略有出入，但变化趋势基本一致.因此，可以认为埋管海床模型能较好地模拟实际工程的应用情况.

图3 埋管海床模型与TS2及TS3的数据结果对比图Fig.3 Comparison of numerical results among seabed model with a buried pipeline, TS2 and TS3

3 计算结果与分析

基于上述已验证的孤立波模型和埋管海床模型，建立孤立波-斜坡海床-海底管道耦合模型；模拟孤立波在斜坡上一个完整的破碎、爬升及回落过程；讨论管道埋置于水平海床、斜坡坡脚及斜坡海岸线附近3种典型情况下，海床和海底管道响应；讨论埋深、土体参数等对海床和海底管道响应的影响.孤立波参数：波高为3 m；d为10 m；ρ为 1 000 kg/m3；Kw为2 GPa.海底管道参数：Ep为68 GPa；υp为0.32；管道密度ρp为 2 700 kg/m3；R为1 m；e为2 m；管壁厚度为0.1 m.斜坡海床参数：海床厚度为50 m；斜坡高度为20 m；水平海床长度为100 m；斜坡长度为200 m；斜坡坡度为1∶10；υs为0.35；ns为0.425；ρs为 2 053 kg/m3；侧向土压力系数为0.41；k为0.1 mm/s；S为1；Gs为5 GPa.

3.1 埋管斜坡海床的响应特征

3.1.1海底管道位置及埋深的影响 海底管道位于3个典型位置，当波峰、波谷经过不同埋深的海底管道圆心正上方时，p沿海底管道圆周分布的情况如图4所示.随着埋深的增加，海底管道上下方孔压峰值的差距逐渐减小，而峰谷值之间的差距却在增加，p的谷值均出现在海底管道两侧.

当海底管道位于水平海床(见图4(a))与斜坡坡脚(见图4(b))时，p的峰谷值之间的差距以及两个峰值之间的差距均较小，不同埋深对p分布的影响也较小；当海底管道位于斜坡海岸线附近时(见图4(c))，p的峰谷值间的差距以及两峰值间的差距均较大，且不同埋深对p分布的影响非常明显；在孤立波回落阶段(见图4(d))，海底管道周围孔压出现负值，管道正下方出现明显的谷值，且随着埋深的增加，管道上下方孔压差距逐渐减小.

综上所述，海底管道的埋深越大，其圆周附近的p越小，管道上下方两个p峰值之间的差距越小，管道正上方p峰值与其侧面谷值之间的差距越大.相比于水平海床和斜坡坡脚位置时的情况，斜坡海岸线附近的海底管道周围p分布更不均匀，且受管道埋深的影响更为显著.

图4 波峰、波谷经过不同埋深管道时p沿管道圆周的分布情况Fig.4 Distribution of p around pipe for different pipe buried depth and location when peak or trough passes

图5 不同k下p及沿纵向及圆周的分布情况Fig.5 The vertical distribution and around pipe distribution of p and for different k

图6 不同S下p及沿纵向及圆周的分布情况Fig.6 The vertical distribution and around pipe distribution of p and for different S

3.2 管道的纵向受力与位移特点

3.2.1管道的纵向受力 在孤立波作用下，当海底管道位于水平海床(x1=-50 m)，斜坡坡脚(x2=0 m)以及斜坡海岸线附近(x3=86 m)3个位置，e=2 m时，管道的纵向受力F随时间的变化趋势如图7所示.其中，F向上为正，向下为负.在波峰到达海底管道前数秒，海底管道受到较小的浮力F1；当波峰经过海底管道时，管道受到向下的压力F2.|F1| 与 |F2| 的值在x3(斜波海岸线附近)处大于x1(水平海床)和x2(斜坡坡脚)处.与x1和x2位置不同，当海底管道埋置于x3处时，在孤立波回落阶段会受到持续约20 s的较大上浮力.这可能是由于该时间段内在海底管道上方位置，即斜坡海岸线附近，出现“水洼”，导致海底管道上方p大幅下降，从而使海底管道受到的上浮力增加.

在孤立波作用下，埋置于x3位置，e=2，4，6，8 m时的F随时间的变化趋势如图8所示.由图8可知，e越小的管道受到的|F|越大.在孤立波回落阶段，e=2 m处的管道受到的最大上浮力远大于其他位置，约为e=4 m处的2倍，e=6 m处的3倍，e=4 m处的4倍.因此，当海底管道埋置在斜坡较浅的位置时，通过增加埋深可以大幅降低在孤立波作用下受到的上浮力.

图7 不同位置海底管道的F随t的变化趋势Fig.7 The variation of F with t for different buried location

图8 x3处不同埋深海底管道的F随t的变化趋势Fig.8 The variation of F near coastline with t for different pipeline depth on x3

3.2.2海底管道的位移 选取海底管道顶端的位移代表海底管道在海床中的整体位移.当海底管道位于3种典型位置，e=2 m时，wp随时间的变化曲线如图9所示.当波峰经过海底管道时，海底管道出现最大沉降，且x3处的海底管道的沉降大于x2及x1处的海底管道的沉降，其值分别为0.032, 0.026与0.021 m.在x3处的波浪回落阶段(t=40～60 s)，受到海底管道上方出现波谷(水洼)的影响，海底管道出现的上浮位移达到0.013 m，海底管道上下浮动的位移差值达到0.045 m，明显大于其他两个位置的位移.

当海底管道处于不同位置且有波峰(谷)经过时，e与wp的关系分布如图10所示.其中：x3(爬升，波峰)对应斜坡海岸线附近位置，孤立波爬升过程中的波峰出现时刻；而x3(回落，波谷)对应斜坡海岸线附近位置，孤立波回落过程中的波谷(水洼)出现时刻.由图10可知，e的大小对wp的极值有一定的影响，并且在3个位置作用相近.e越大，|wp|越小，且基本呈线性变化.

图9 不同位置，wp随t的变化曲线Fig.9 The variation of wp with t for different buried locations

图10 不同位置有波峰(谷)经过时wp随e的变化Fig.10 The variation of wp with e for different buried locations when peak or trough passes

4 结论

(1)e越大，管道圆周附近|p|越小.p分布越不均匀，管道周围p的峰值与谷值之间的差距越大.相比于水平海床和斜坡坡脚位置，斜坡海岸线附近管道周围p的分布更不均匀，且受到e的影响更为显著.

(5) 海岸线附近位置并且e较浅的海底管道在孤立波回落阶段会受到较大的、持续的上浮力，可以通过增加e，大幅降低浮力的最大值，从而缓解不利情况.

(6)wp主要受管道埋置位置的影响.其中，海底管道在斜坡海岸线附近出现最大上浮和下沉，最大位移差为0.045 m，该数值远大于水平海床与海床斜坡坡脚位置的最大位移差.此外，e越大，|wp|越小.