椭圆余弦波作用下海床动态响应及液化分析

张军，周香莲，2，3，王建华

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



椭圆余弦波作用下海床动态响应及液化分析

张军1，周香莲1，2，3，王建华1

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

摘要:对于相对波高小于1/20浅水地区，传统线性理论不再适用，应该选用适用性更好的椭圆余弦波理论，但由于椭圆余弦函数计算的难度较大，椭圆余弦波理论的应用一直有待推广。本文采用精细积分算法模拟椭圆余弦波，结合考虑土骨架加速度及流体速度的u－p动力模式，建立椭圆余弦波－海床模型，得到椭圆余弦波作用下海床的孔压及竖向有效应力分布情况，进行液化分析。对渗透系数，饱和度等参数进行研究可以得出椭圆余弦波对孔隙水压力的影响比线性波大。

关键词:椭圆余弦波;海床;动态响应;液化;参数研究

近年，我国近海工程不断增多，近海海域(相对波高小于1/20)波浪作用下海床动态响应研究有待深入，近海海域应该选用更符合实际的椭圆余弦波理论，而非传统线性波理论。

学者们对于椭圆余弦波的模拟和应用都进行了较多的研究，其中，Fenton［1］系统地研究了椭圆余弦波，提出椭圆余弦波理论在浅水条件下具有非常好的适用性，但由于雅可比椭圆函数和椭圆积分的计算难度较大，使得其发展和应用都受到了局限。沈先荣［2］总结了椭圆余弦波理论的研究工作进展且简评了数值计算的几种方法。姚征等［3］提出了椭圆函数的精细积分改进算法。Chang等［4］通过PIV技术获得了椭圆余弦波的实测数据。Xu等［5］提出了椭圆余弦函数的近似表达式及其可以接受的使用条件。Xu等［6］基于准静态的Biot固结模型继续研究了椭圆余弦波作用下海床孔压及有效应力的分布。

对于海床动态响应的研究，Biot［7］将其固结理论推广到Partly dynamic模型。Hsu等［8］研究了波浪作用下细砂质海床动态响应问题。Ulker［9］得出了波浪作用下饱和多孔海床响应的一组广义解析解。Zen［10］基于超静孔隙水压力提出了液化判断准则。Gao等［11］求解了非线性波浪作用下海床的动态响应问题。Xu［12］研究了海床在波浪作用下的液化进程。Wen等［13］研究了非线性波流共同作用下埋管海床的动态响应问题。针对多孔多层海床，Zhou等［14］研究土壤特性对波浪引起的海床动态响应的影响。Zhou等［15］研究了波流共同作用下各向异性海床的动态响应问题。

本文使用Comsol Multiphysics软件，基于Biot动力方程，将海床考虑为多孔弹性介质，采用u－p模式将孔压位移视为场变量，考虑土体位移加速度，忽略孔隙流体惯性项作用，结合采用精细积分法得到的椭圆余弦波，建立波浪－海床模型，运用隐式欧拉向后差分法(BDF法)进行求解。模型得到验证后进行参数分析，研究了渗透系数、饱和度对于孔压和竖向有效应力及液化区域的影响。

1　椭圆余弦波作用下海床动态响应理论分析

1.1计算模型的建立

计算模型如图1所示，厚度为h的海床之下为固定的不透水基层，水深为d，波高为H，波长为L的椭圆余弦波沿x轴正向传播。

图1　椭圆余弦波－海床计算模型Fig.1 Sketch of the present model

1.2椭圆余弦波荷载的计算

椭圆余弦波，是指有限水深条件下，具有稳定的有限振幅的长周期波。本文使用Xu等提出的改进精细积分算法进行椭圆余弦波的模拟。如图1所示，原点固定于静水面上的一阶椭圆余弦波波面方程为:

式中:H为波高;L为波长;T为周期;t为时间;cn(·)表示椭圆余弦函数;m为模量;E和K分别为第一类和第二类完全椭圆积分，可由下式计算得到:

对于模量m的计算已经有较为成熟的方法，在此不再赘述，而对于椭圆余弦函数的计算则较复杂，选用改进后的精细积分算法进行计算，大致的思路如下:计算时，对于任意的，定义一个，使得

浅水地区的海床表面的波压力pb可表示为:

式中:γw为孔隙流体的重度，k为波数。单位化参数P0取为γwH。如图2所示，本方法计算出来的椭圆余弦波波面结果，与Chang等［4］的实测及其模拟结果均吻合良好，选取的计算参数为:周期T=2 s，波高H=0.036 m，水深d=0.24 m。图中可见，椭圆余弦波的能量集中于波峰，波峰大大高于静水位，而波谷则离静水位很近。

图2　本文椭圆余弦波模拟结果与文献［4］的实测及数值模拟结果对比Fig.2 Comparison of present and previous results between this papers and reference［4］

1.3椭圆余弦波作用下土体响应分析

考虑土体为各向同性介质，渗透系数为常数，Biot方程可表示为:

式中:p为孔隙水压力;γw为孔隙流体的重度;n为海床土壤孔隙率;Kz为土壤渗透系数;ρf为流体的密度;β为孔隙流体的压缩系数;ε为土体体积应变。β、ε分别由以下两式定义:

式中:kw为孔隙水的体积模量，通常取kw=2× 109N/m2;S为土体饱和度;u和w分别为土体位移在水平及竖直方向上的分量;Pw0为孔隙水的绝对压力。

简化了的u－p形式的控制方程为:

式中:σij为土体总应力;ρ和ρf分别是土体和孔隙流体的密度;为土体骨架的加速度;为孔隙流体相对于土体骨架的平均位移;Kf为孔隙流体的体积模量。

2　椭圆余弦波作用下海床动态响应模型验证及参数分析

2.1模型验证

为了验证本文所采用模型的正确性，将本文所用的椭圆余弦波模型退化到使用Stokes波，将结果与解析解进行比较，如图3，实线为当前模型，点为解析解;输入的数据为北海波浪设计值:T=15 s，d=70 m，L=311.59 m，h=25 m，G=107N/m3，μs=1/3，n=0.3，S=1.0，Kz=0.01 m/s。可见本文模型与解析解吻合较好，验证了本文结果的有效性。

图3　当前模型和解析解文献［8］的比较Fig.3 Comparison of previous result between present model and reference［8］

2.2参数分析

本部分研究了饱和度，渗透系数对于海床孔压，竖向有效应力和竖向位移沿海床深度分布的影响。标准参数见表1、2，为了保证结果的合理性，选取第二周期末的时间点进行分析，使用单位化后的孔压p/P0，竖向有效应力σ'zz/P0，竖向位移uz2G/P0进行描述。为研究单独某项参数对于动态响应结果的影响，不考虑各参数之间的联动情况。

表1　海水参数Table 1 Seawater parameters

表2　海床参数Table 2 Seabed parameters

2.2.1饱和度S的影响

饱和度是海床土体的重要参数，图4～6描述了不同的饱和度(S=1，0.98，0.96，0.94)下孔压，竖向有效应力以及竖向位移沿海床深度的分布。图7为孔压分布云图。

图4　不同饱和度下孔压沿海床深度的分布Fig.4 Pore pressure distribution versus seabed depth with different degree of saturations

如图4所示，孔压随着深度增加逐渐减小，减小过程中梯度逐渐增大，随着饱和度的增大，孔压减小得更慢。

图5　不同饱和度下竖向有效应力沿海床深度的分布Fig.5 Vertical stress distribution versus seabed depth with different degree of saturations

图5表明竖向有效应力呈现先减小后增大的趋势，随着饱和度减小，竖向有效应力峰值增大。图6中可见竖向位移随着海床深度增大而逐渐的减小直至出现反向位移，竖向位移随着饱和度的减小而增大。

图6　不同饱和度下竖向位移沿海床深度的分布Fig.6 Vertical displacement distribution versus seabed depth with different degree of saturations

图7　不同饱和度下的孔压分布云图Fig.7 Contour of pore pressure distribution with different degree of saturations

图7描述了整体海床中的孔压分布，可见随着饱和度的减小，孔压沿深度方向减小更快。

2.2.2渗透系数Kz的影响

图8～10给出了不同渗透系数下沿海床深度方向的孔压、竖向有效应力、竖向位移。如图8所示，孔压随着海床深度的增加而逐渐减小，渗透系数越大，减小得越慢。由图9，竖向有效应力沿深度方向呈现先减小后增大的趋势，渗透系数越大，竖向有效应力峰值越小。图10给出了位移随海床深度的变化，可见，渗透系数越小，海床表面产生的竖向位移越大。

图8　不同渗透系数下孔压沿海床深度的分布Fig.8 Pore pressure distribution versus seabed depth with different permeability

图9　不同渗透系数下竖向有效应力沿海床深度的分布Fig.9 Vertical stress distribution versus seabed depth with different permeability

图10　不同渗透系数下竖向位移沿海床深度的分布Fig.10 Vertical displacement distribution versus seabed depth with different permeability

2.2.3不同波浪理论的影响

对线性波浪理论和椭圆余弦波理论下孔压的变化做出比较。由图11可见，不同的波浪理论对于海床表面的孔压影响较大，孔压随深度减小趋势基本相同，但由于两种波浪作用下的海床表面所受波动水压力差异较大，导致靠近海床表面部分，椭圆波作用下海床的孔压明显大于线性波作用下的孔压，因此在计算浅海地区波浪作用下的海床动态响应时不能简单地使用线性波理论进行简化。

图11　不同波浪理论下孔压沿深度分布Fig.11 Pore pressure distribution versus seabed depth with different wave theories

3　液化判断准则及液化参数分析

3.1液化判断准则

在波浪荷载持续作用下孔隙水压力上升，有效应力下降，海床的一部分可能变得不稳定而发生液化，本文采用由Zen［10］提出的基于超静孔隙水压力的液化判断准则，认为当土层中某一点处上层土体骨架的重量小于超静孔隙水压力时海床发生液化，即:

3.2液化参数分析

使用液化准则分析不同饱和度S，渗透系数Kz下的液化区域，图中深色区域表示发生液化的土体。

3.2.1饱和度S的影响

饱和度对于海床是否液化有着很重要的影响。选取S=1、0.97、0.94这3种不同参数进行分析，如图12所示，随着饱和度的减小，液化深度和液化横向范围都逐渐增大。S=1时不发生液化。从S=0.96到S=0.94，液化深度从1.2 m变化到1.5 m，液化横向范围从25 m增大到28 m。

图12　不同饱和度下液化区域的变化Fig.12 The liquefaction location with different degree of saturations

3.2.2渗透系数Kz的影响

针对不同的渗透系数(Kz=1×10－3m/s，1× 10－4m/s，1×10－5m/s)，图13给出了其液化区域。可见随着渗透系数的减小，液化深度和液化横向范围逐渐增大。

Kz=1×10－3m/s时土体不发生液化，从Kz=1× 10－4m/s到Kz=1×10－5m/s，液化深度从0.6 m变化到1 m，液化横向范围从20 m增大到27 m。这一现象具有较大的工程意义。表明近海中表明近海中海床土体的渗透系数较大时，液化的可能性将降低。在工程实践中可以采取铺设渗透系数较大的粗粒料上覆层作为保护管道的措施。

图13　不同渗透系数下液化区域的变化Fig.13 The liquefaction location with different permeability

4　结论

本文使用精细积分法计算椭圆余弦函数，模拟一阶椭圆余弦波，结合u－p动力模式建立波浪－海床动态响应模型。研究了渗透系数、饱和度对孔压、竖向有效应力、竖向位移以及液化区域的影响。可得到以下结论:

1)由于椭圆余弦波和线性波作用下的海床表面的所受波动水压力差异很大，椭圆波作用下孔压明显大于线性波作用下的孔压。

2)随着饱和度的增大，孔压减小得更慢，竖向有效应力峰值减小，竖向位移减小，液化区减小。

3)渗透系数越大，孔压减小得越慢，竖向有效应力峰值越小，竖向位移越小，液化区越小。

参考文献:

［1］FENTON J D.The cnoidal theory of water waves［M］.Houston:Gulf Professional Publishing，1998:1－34.

［2］沈先荣.椭圆余弦波理论研究工作的进展［J］.海洋工程，1990，8(4):77－87.SHEN Xianrong.Advances of studies on cnoidal wave theory ［J］.The Ocean Engineering，1990，8(4):77－87.

［3］姚征，钟万勰.椭圆函数的精细积分改进算法［J］.数值计算与计算机应用，2008，29(4):251－260.YAO Zheng，ZHONG Wanxie.The improved precise integration method for elliptic functions［J］.Journal on Numerical Methods and Computer Applications，2008，29(4):251－260.

［4］CHANG K A，HSU T J，LIU P L F.Vortex generation and evolution in water waves propagating over a submerged rectangular obstacle:Part II:Cnoidal waves［J］.Coastal Engineering，2005，52(3):257－283.

［5］XU Yunfeng，XIA Xiaohe，WANG Jianhua.Calculation and approximation of the cnoidal function in cnoidal wave theory ［J］.Computers＆Fluids，2012，68:244－247.

［6］XU Yunfeng，XIA Xiaohe，WANG Jianhua，et al.Numerical analysis on cnoidal wave induced response of porous seabed with definite thickness［J］.Journal of Shanghai Jiaotong University(Science)，2013，18(6):650－654.

［7］BIOT M A.Theory of propagation of elastic waves in a fluidsaturated porous solid.I.low－frequency range［J］.The Journal of the Acoustical Society of America，1956，28(2):168－178.

［8］HSU J R C，JENG D S.Wave－induced soil response in an unsaturated anisotropic seabed of finite thickness［J］.International Journal for Numerical and Analytical Methods in Geomechanics，1994，18(11):785－807.

［9］ULKER M B C，RAHMAN M S，JENG D S.Wave－induced response of seabed:various formulations and their applicability［J］.Applied Ocean Research，2009，31(1):12－24.

［10］ZEN K，YAMAZAKI H.Mechanism of wave－induced liquefaction and densification in seabed［J］.Soils and Foundations，1990，30(4):90－104.

［11］GAO F P，JENG D S，SEKIGUCHI H.Numerical study on the interaction between non－linear wave，buried pipeline and non－homogenous porous seabed［J］.Computers and Geotechnics，2003，30(6):535－547.

［12］XU Haixia.Wave－induced liquefaction processes in marine sediments［D］.Dundee:The University of Dundee，2012:39－62.

［13］WEN F，JENG D S，WANG J H，et al.Numerical modeling of response of a saturated porous seabed around an offshore pipeline considering non－linear wave and current interaction［J］.Applied Ocean Research，2012，35:25－37.

［14］ZHOU Xianglian，XU Bin，WANG Jianhua，et al.An analytical solution for wave－induced seabed response in a multi－layered poro－elastic seabed［J］.Ocean Engineering，2011，38(1):119－129.

［15］ZHOU Xianglian，WANG Jianhua，ZHANG Jun，et al.Wave and current induced seabed response around a submarine pipeline in an anisotropic seabed［J］.Ocean Engineering，2014，75:112－127.

［16］JENG Dongsheng.Porous models for wave－seabed interactions［M］.Berlin:Springer－Verlag，2013:7－32.

Numerical study of the cnoidal wave-induced dynamic response of the seabed and liquefaction analysis

ZHANG Jun1，ZHOU Xianglian1，2，3，WANG Jianhua1
(1.Center for Marine Geotechnical Engineering，Shanghai Jiao Tong University，Shanghai 200240，China;2.Collaborative Innovation Center for Advanced Ship and Deep－Sea Exploration，Shanghai 200240，China;3.State Key Laboratory of Ocean Engineering，Shanghai JiaoTong University，Shanghai 200240，China)

Abstract:Cnoidal wave theory can be applied to periodic wave progression in water when the wave height is less than 1/20 of the wavelength and for which traditional linear theory is no longer appropriate.However，cnoidal wave theory has not been widely applied in practical engineering applications because the formula for the wave profile involves a Jacobian elliptic function.In this study，Taylor's expression and a precise integration method were used to estimate the Jacobian elliptic function.The seabed was treated as a porous medium and characterized by Biot's partly dynamic equations(u－p model).A cnoidal wave－seabed system was modeled and analyzed.Based on the numerical results，the effects of seabed characteristics such as permeability and degree of saturation on the cnoidal waveinduced excess pore pressure and liquefaction phenomenon were studied.It was concluded that the effect of a cnoidal wave on the pore water pressure is larger than that of a linear wave.

Keywords:cnoidal wave;seabed;dynamic response;liquefaction analysis;parametric study

通信作者:周香莲，E－mail:zhouxl@ sjtu.edu.cn.

作者简介:张军(1990－)，男，硕士;周香莲(1972－)，女，副教授，博士.

基金项目:国家自然科学基金资助项目(41372286，41572243);海洋工程国家重点实验室青年创新基金资助项目(GKZD010059).

收稿日期:2014－08－13.网络出版时间:2015－12－21.

中图分类号:P75;TU4

文献标志码:A

文章编号:1006－7043(2016)01－0053－06

doi:10.11990/jheu.201408014

网络出版地址:http://www.cnki.net/kcms/detail/23.1390.u.20151221.1555.028.html