华莹，周香莲，，张军

（1.高新船舶与深海开发装备协同创新中心（船海协创中心)，上海200240；2.上海交通大学船舶海洋与建筑工程学院海洋工程国家重点实验室，上海200240）

随机波作用下埋管海床动态响应及液化研究

华莹1，周香莲1，2，张军2

（1.高新船舶与深海开发装备协同创新中心（船海协创中心)，上海200240；2.上海交通大学船舶海洋与建筑工程学院海洋工程国家重点实验室，上海200240）

基于广义Biot动力理论和Longuet-Higgins线性叠加模型，构建波浪-海床-管线动态响应的有限元计算模型，求解随机波作用下，多层砂质海床中管线周围土体孔隙水压力和竖向有效应力的分布。采用基于超静孔隙水压力的液化判断准则，得出液化区的最大深度及横向范围，从而判断海床土体液化情况。考虑海洋波浪的随机性，将海床视为多孔介质，海床动态响应计算模型采用u-p模式，孔隙水压力和位移视为场变量。并考虑孔隙水的可压缩性、海床弹性变形、土体速度、土体加速度以及流体速度的影响，忽略孔隙流体惯性作用。参数研究表明：土体渗透系数、饱和度以及有效波高等参数对海床土体孔隙水压力、竖向有效应力和液化区域分布有显著影响。

随机波；海床；管线；动态响应；液化；参数分析

管线运输是目前最安全快捷和经济可靠的海上运输方式，由于海洋环境的特殊性海底管线受力复杂，因此埋管海床的动态响应研究是管线在建设与运营过程中的关键技术问题之一。最初海床动态响应分析大多采用相对简单的线性波，Madsen（1978）和Yamamoto（1978）分别推导出海床动态响应水力各向同性和水力各向异性的解析解。Jeng（1997）推导出准静态模式下海床动态响应的解析解，并和Seymour（1997）研究了短峰波和变渗透系数对有限厚度海床动态响应的影响。钟佳玉等（2009）和李安龙等（2012）分别对波浪作用下砂质和粉土海床孔压响应做了试验研究。

实际波浪参数具有高度的不规则性。对于随机波的研究，Longuet－Higgins（1957）分析电子管噪声电流，通过叠加线性波模拟出随机波的传递过程。Hasselmann等（1973）在“联合北海波浪计划”中，由实测海浪谱分析拟合得到JONSWAP谱，包含反映能量水平、峰的频率尺度和谱形的参量。

随着波浪理论的发展，随机波开始应用于海床动态响应研究中，但对波浪－海床－管线相互作用的分析仍较少。Sumer（1999）实验模拟海床在随机波作用下的变形及响应，表明孔压形成过程与规则波作用的规律大致相同。王忠涛等（2008）采用复变量解析法，考虑波浪荷载随机特性建立了海床动态响应与液化的解析数值模型。Jeng（2012）阐述了多孔介质理论下的随机波浪和海床交互作用。邓海峰等（2014）研究了随机波浪、管－土接触和惯性效应对海床动态响应的影响。张军等（2015）将随机波与一阶Stokes波和椭圆余弦波进行对比，分析不同波浪理论对海床动态响应的影响。

本文考虑海洋波浪的随机性，采用Longuet－Higgins模型模拟随机波；考虑海床土体分层和管线的作用，使用COSMOL多物理场耦合软件，构建波浪－海床－管线动态响应的有限元模型，确定多孔弹性海床中管线周围的孔压和竖向有效应力大小及分布；基于超静孔隙水压力的液化判断准则，得出直观的液化区分析。本文基于Biot（1941）动力理论，采用u-p模式，海床视为多孔介质，孔压和位移视为场变量，考虑孔隙水的可压缩性、海床弹性变形、土体速度、土体加速度以及流体速度的影响，忽略孔隙流体惯性作用。

1 随机波浪模拟及海床模型分析

1.1 随机波浪的模拟

采用Longuet－Higgins线性叠加模型，海浪谱或者水面高度可表示为：

式中：M为组成波的数量，数值充分大；下标i代表第i个组成波；ai为波浪幅值；为波浪频率，在（fi－1，f）i中随机分布；ζi为初始波相角，在（0，2π）中随机分布。在（1）式中，知道第i个组成波的频率和水深d之后，波数ki可由下式求得：

式中：g为重力加速度。组成波的波浪幅值ai可由已知的随机波频谱公式S（f）求得：

采用改进Jonswap频谱（Gōda，2010）作为S（f）：

其中：

式中：Tp为频谱峰值点对应的波浪周期，fp为频谱峰值点对应的波浪频率，Tp=1/fp。γ为谱峰增强因子，通常情况下γ=1～7（平均值为3.3）。TH1/3和H1/3分别为波浪的有效周期和有效波高。

1.2 海床动态响应有限元模型分析

图1 海床管线相互作用基本模型

如图1，波浪沿x轴方向传播，水深为d。刚性不透水基岩上有厚度为h的多孔饱和砂质海床，海床等分为三层，每层各向同性。第1层海床中埋有管径为D的管线，埋深为e。

本文使用COMSOL多物理场耦合软件PDE模块进行有限元分析，海床模型尺寸为200 m×30 m，离散化处理后最大Lagrange单元尺寸0.3 m，单元总数N=2 335。求解孔压和竖向有效应力数值及分布，管线周围的网格作了局部细化处理，如图2。

图2 模型网格划分

1.2.1 控制方程

每一层海床土体渗透系数为常数，考虑超静孔隙水压力和土骨架变形的Biot固结方程为：

式中：Kz为土体渗透系数；p为孔隙水压力；ρf为流体的密度；ε为土体体积应变；γw为孔隙流体的重度；n为土体孔隙率；β为孔隙流体压缩性系数。

式中：u和w为土体位移在水平和竖直方向上的分量；kw为孔隙水体积模量，一般取kw=2×109N/m2；S为海床土体的饱和度；Pw0为孔隙水的绝对压力。

根据Biot固结方程，Zienkiewicz等（1980）忽略孔隙流体惯性项，得到u-p形式的控制方程：

式中：σij为总应力；ρ为海洋土密度；ρw为孔隙水密度；Kf为孔隙水的压缩模量。

1.2.2 边界条件

（1）海床表面：忽略水粘性和摩擦，超静孔隙水压力等于波压力且竖向有效应力和剪应力满足：

随机波作用下海床表面受到的波压力pb为：

得：

（2）海床底部：视为刚性不透水基岩，海床土体骨架没有位移，且法向流量为零：

（3）海床侧边界：土体骨架水平位移、水平流量为零：

（4）管线表面：视为不透水边界，则超静孔隙水压力梯度为零：

（5）三层海床接触界面：位移、孔压和竖向有效应力相等，第i层和第i+1层之间连续条件为：

1.2.3 模型验证

一方面，不考虑管线将随机波作用下沿海床深度的孔压结果与Liu等（2007）的解析解结果对比，具体参数如图3。本文结果与解析解相比，波面压力值几乎相同，变化趋势基本相同。存在差异主要因为本文使用Biot动态模型，Liu等使用准静态的Biot固结模型。

图3 当前模型和解析解（Liu等（2007））的比较

图4 当前模型和（邓海峰等）的比较

另一方面，对于随机波作用下海床-管线模型也进行了验证，具体参数取值如图4。对比本模型和邓海峰等（2014）的模型中管线圆周的孔压分布，变化趋势基本相同，数值上有差异主要因为邓海峰等假定管线为线弹性材料。综合两方面，表明本文数值模型得到的海床动态响应合理可靠。

2 海床动态响应及液化分析

考虑管线作用，研究渗透系数、饱和度和有效波高等参数，对孔压和竖向有效应力沿管线圆周和过管线中心剖面分布的影响。选取表1中的参数作为标准参数。

表1 计算参数的选取

液化是波浪作用下砂质海床响应的极限状态（赵子丹等，1995），本文采用基于超静孔隙水压力的液化判断准则（Zen et al，1990），即：当土层中某一点处上层土体骨架的重量小于超静孔隙水压力时海床土体发生液化：

式中，γs为海洋土天然重度；pb为海床表面波压力；p为海床内深度为z处的孔压。本文判断液化程度的标准包括液化区的最大深度和横向范围。以波谷传播到管线上方海床表面时作为液化区域分析点，同时数值计算从波谷点前取5个时间步长开始，考虑了波压力向海床内部传递的时滞性。

2.1 渗透系数的影响

渗透系数是衡量海床排水能力的重要参数，选取三组参数Kz=0.1，0.01，0.001（m/s）、Kz=0.01，0.01，0.001（m/s）和Kz=0.001，0.001，0.001（m/s）分别代表三层、双层和单层海床，其余参数均相同。

不同渗透系数下，无量纲化后的孔压及竖向有效应力沿管线圆周的分布如图5。图5（a）表明随着渗透系数降低，管线周围孔压明显减小，峰值均在90°附近（记管线最左侧为0°处，顺时针为正方向）；图5（b）表明随着渗透系数降低，竖向有效应力增大，但变化的幅度较小。

图5 不同渗透系数下孔压及竖向有效应力沿管线圆周分布

过管线中心剖面下孔压和竖向有效应力沿深度的分布如图6。图6（a）表明第1层海床土体渗透系数降低导致孔压衰减更快，到一定深度后差距减小，海床底部孔压几乎为0，且海床土体分层对孔压的分布影响明显；图6（b）表明随着第1层土体渗透系数的减小，竖向有效应力在浅层海床中增长更快，到管线以下差距逐渐缩小直至重合。

图6 不同渗透系数下孔压及竖向有效应力沿深度分布

在第1层海床中，选取不同的渗透系数Kz=0.005（m/s）、Kz=0.003（m/s）和Kz=0.001（m/s），图7给出了相应的液化区域。其中图中云线含义为超静孔隙水压力与土层某点处上层土体骨架的重量之差由图7得Kz=0.005（m/s）时不发生液化，Kz=0.001（m/s）时液化最大深度达2.0 m，横向宽度达14.1 m。表明渗透系数减小，液化区的最大深度和横向宽度均增大。同时，本文的计算结果对比Jeng（1997）的结论：渗透系数Kz＞0.001 m/s时海洋土体不发生液化，验证了其合理性。

图7 不同渗透系数下的液化区

2.2 饱和度的影响

选取三组不同的饱和度参数S=0.99，0.97，0.95、S=0.97，0.97，0.95和S=0.95，0.95，0.95，分别代表三层、双层和单层海床土体，其余参数均相同。

无量纲化的孔压和竖向有效应力沿管线圆周的分布如图8。图8（a）表明饱和度降低，管线周围孔压均降低，且降低的幅度沿管线圆周较均匀，峰值均在90°附近；图8（b）表明饱和度减小导致管线周围竖向有效应力增大，但变化的幅度较小。

过管线中心剖面的孔压和竖向有效应力沿深度的分布如图9。图9（a）反映出第1层土体的饱和度越小孔压衰减越快，到管线以下差距逐渐减小直至重合；图9（b）反映出竖向有效应力主要在第1层海床中有不同，饱和度越大，管线以上竖向有效应力峰值越小，管线以下差距逐渐减小直至重合。

在第1层海床中，选取不同的饱和度参数S=0.99、S=0.98和S=0.97，土体液化情况如图10。S=0.99时不发生液化，S=0.97时液化最大深度达1.9 m，横向宽度达12.6 m。结果表明：随着饱和度的减小，液化区的最大深度和横向宽度均增大。本文计算结果与王忠涛等（2008）的结论“海床最大液化深度在96.5%位置处出现峰值，随着饱和度的逐渐增加，液化深度逐渐减小，当土体完全饱和时海床土体不再发生液化”相一致。

图9 不同饱和度下孔压及竖向有效应力沿深度分布

图10 不同饱和度下的液化区

2.3 有效波高的影响

本部分通过选取不同的有效波高H1/3=1.0 m、H1/3=1.5 m和H1/3=2.0 m，第1、2、3层海床土体参数均相同，分析波高对孔压和竖向有效应力的影响。

无量纲化的孔压及竖向有效应力沿管线圆周分布如图11。图11（a）表明有效波高增大，管线周围孔压明显增大；管线90°附近为孔压峰值，且增长最为显著。图11（b）表明随着有效波高的增大，管线周围竖向有效应力变小，且分布更均匀。

过管线中心剖面的孔压及竖向有效应力沿深度分布如图12。图12（a）表明有效波高越大，海床表面孔压越大衰减越快，随着深度的增加差距逐渐缩小；图12（b）表明竖向有效应力的差距主要在第1层海床，在管线附近的竖向有效应力数值均随着有效波高的增大而增大，分别向海床两端收敛。

在第1层海床中，选取不同的有效波高参数H1/3=1.0 m、H1/3=1.5 m和H1/3=2.0 m，土体液化情况如图13。H1/3=1.0 m时液化最大深度为0.6 m，横向宽度为11.5 m；H1/3=2.0 m时液化最大深度达1.45 m，横向宽度达15.5 m。表明有效波高增大，液化区的最大深度和横向宽度均增大

图11 不同有效波高孔压及竖向有效应力沿管线圆周分布

图12 不同有效波高孔压及竖向有效应力沿深度分布

图13 不同有效波高下的液化区

3 结论

考虑随机特性的波浪接近实际情况，在波压力作用下海床土体液化，引起海底管线失稳甚至断裂破坏，因此海底管线的设计需考虑诸多因素，并对铺设管线的海洋环境作出综合的考量。

（1）渗透系数大，孔隙水压力易消散，故增大渗透系数可有效减小海床液化的可能性，为实际工程应用提供指导：尽量选择渗透系数较大的海床路径，当无法避免时可以采用渗透系数较大的粗粒料作为上覆层再铺设管线。

（2）一般情况下砂土饱和度越高越容易液化，但对于接近饱和的海床土，饱和度越高反而不容易液化。原因是：当饱和度很高时，随机波浪荷载作用下波压力沿深度方向传递时衰减慢，形成的压差小，导致产生超孔压的吸力减小。

（3）随机波作用下海床-管线动态响应中，有效波高的影响非常显著。有效波高增大会使管线周围的孔压和有效应力均明显变化，液化的深度和横向范围也明显扩大，对管线安全极为不利。

Biot M A,1941.General theory of three-dimensional consolidation.Journal of Applied Physics,12(2):155-164.

Gōda Y,2010.Random seas and design of maritime structures.World Scientific.

Hasselmann K,Barnett T P,Bouws E,et al,1973.Measurements of windwave growth and swell decay during the Joint North Sea Wave Project(JONSWAP).Deutches Hydrographisches Institut.

Jeng D S,1997.Wave-induced seabed response in front of a breakwater.University of Western Australia.

Jeng D S,2012.Porous models for wave-seabed interactions.Springer Science&Business Media.

Jeng D S,Seymour B R,1997.Response in seabed of finite depth with variable permeability.Journal of Geotechnical and Geoenviron mental Engineering,123(10):902-911.

Liu H,Jeng D S,2007.A semi-analytical solution for random waveinduced soil response and seabed liquefaction in marine sediments.Ocean Engineering,34(8):1 211-1 224.

Longuet-Higgins M S,1957.The statistical analysis of a random,moving surface.Philosophical Transactions of the Royal Society of London A:Mathematical,Physical and Engineering Sciences,249(966):321-387.

Madsen O S,1978.Wave-induced pore pressures and effective stresses in a porous bed.Geotechnique,28(4):377-393.

Sumer B M,Fredsøe J,Christensen S,et al,1999.Sinking/floatation of pipelines and other objects in liquefied soil under waves.Coastal Engineering,38(2):53-90.

Yamamoto T,Koning H L,Sellmeijer H,et al,1978.On the response of a poro-elastic bed to water waves.Journal of Fluid Mechanics,87(1):193-206.

Zen K,Yamazaki H,1990.Mechanism of wave-induced liquefaction and densification in seabed.Soils and Foundations,30(4):90-104.

Zienkiewicz O C,Chang C T,Bettess P,1980.Drained,undrained,consolidatinganddynamicbehaviourassumptionsinsoils.Geotechnique,30(4):385-395.

邓海峰，王忠涛，刘鹏，2014.随机波浪作用下海底管线与海床的相互作用研究.水利与建筑工程学报，12(4)：43-49.

李安龙，李广雪，林霖，等，2012.波浪作用下粉土海床中的孔压响应试验研究.海洋通报，31(1)：15-20.

王忠涛，栾茂田，Jeng D S，等，2008.随机波浪作用下海床动力响应及液化的理论分析.岩土力学，29(8)：2 051-2 076.

张军，江俊达，华莹，等，2015.随机波作用下海床动态响应分析.海洋工程，33(3)：80-85.

赵子丹，别社安，1995.沙质海床对波浪的响应及其稳定性研究的回顾.海洋通报，14(4)：85-104.

钟佳玉，郑永来，倪寅，2009.波浪作用下砂质海床孔隙水压力的响应规律实验研究.岩土力学，30(10)：3 188-3 193.

Numerical study of random wave induced seabed-pipeline response and liquefaction

HUA Ying1,ZHOU Xiang-lian1,2,ZHANG Jun2

（1.Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration（CISSE),Shanghai 200240,China;2.State Key Laboratory of Ocean Engineering,Shanghai Jiao Tong University,Shanghai 200240,China）

Based on general Biot's dynamic consolidation theory and Longuet-Higgins random wave theory,a finite element calculating model was established for random wave induced seabed-pipeline dynamic response and liquefaction.The model simulated the distribution of pore pressure,vertical effective stress surrounding the pipeline under random wave.The liquefaction criteria based on excess pore water pressure was adopted to obtain the maximum depth and crosswise range of liquefaction zone.The stochastic feature of ocean wave loading had been taken into account and the seabed was treated as porous medium.In the calculation model u-p mode formulations were adopted,in which both pore pressure and displacement were defined as field variables.Besides,the influences from compressibility of pore water,elastic deformation of seabed,soil velocity,soil acceleration and fluid velocity were considered.However the inertia effect of pore fluid was ignored.The results show that the parameters of seabed,including the permeability,degree of saturation and effective wave height,influence the distribution of pore pressure,vertical effective stress and liquefaction zone significantly.

random wave;seabed;pipeline;dynamic response;liquefaction;parameter analysis

P756.2

A

1001－6932（2017）06－0644－08

10.11840/j.issn.1001-6392.2017.06.006

2016-05-23；

2016-10-11

国家自然科学基金（41372286；41572243）；海洋工程国家重点实验室（GKZD010059）

华莹（1994-），硕士研究生。主要从事海洋+土力学方面的研究。电子邮箱：mochahy@sjtu.edu.cn。

周香莲（1972-），副教授。电子邮箱：zhouxl@sjtu.edu.cn。

袁泽轶）