2. 高新船舶与深海开发装备协同创新中心, 上海 200240;
3. 上海交通大学 船舶海洋与建筑工程学院 海洋工程国家重点实验室, 上海 200240
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
近年,我国近海工程不断增多,近海海域(相对波高小于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.2 椭圆余弦波荷载的计算
椭圆余弦波,是指有限水深条件下,具有稳定的有限振幅的长周期波。本文使用Xu等提出的改进精细积分算法进行椭圆余弦波的模拟。如图1所示,原点固定于静水面上的一阶椭圆余弦波波面方程为:
$$\eta = H\left( {{1 \over {{m^2}}} - 1 - {E \over {{m^2}K}}} \right) + H{\rm{c}}{{\rm{n}}^2}\left[ {2K\left( {{x \over L} - {t \over T}} \right)} \right]$$
(1)
$$k = \int_0^{{\pi \over 2}} {{1 \over {\sqrt {1 - {m^2}{{\sin }^2}} \varphi }}} {\rm{d}}\varphi $$
(2)
$$E = \int_0^{{\pi \over 2}} {\sqrt {1 - {m^2}{{\sin }^2}\varphi } } {\rm{d}}\varphi $$
(3)
对于模量m的计算已经有较为成熟的方法,在此不再赘述,而对于椭圆余弦函数的计算则较复杂,选用改进后的精细积分算法进行计算,大致的思路如下:
计算${\rm{cn}}\left[ {2K\left( {{x \over L} - {t \over T}} \right)} \right]$时,对于任意的${\left( {{x \over L} - {t \over T}} \right)}$,定义一个ξ(ξ≤2),使得
$${\rm{cn}}\left[ {2K\left( {{x \over L} - {t \over T}} \right)} \right] = {\rm{cn}}(2K\xi ,m)$$
(4)
应该找到满足式(4)的ξ,使得ζ=2Kξ210,之后通过椭圆余弦函数的泰勒展开式计算sn(ζ,m)、cn(ζ,m)、dn(ζ,m),再通过椭圆函数的二倍公式计算得到S1=sn(2ζ,m),C1=cn(2ζ,m),D1=dn(2ζ,m),重复第二步10次可得到S10,C10,D10。最后得到y=cn(2Kξ,m)=1-C10。
浅水地区的海床表面的波压力pb可表示为:
pb=γw·η
(5)
考虑土体为各向同性介质,渗透系数为常数,Biot方程可表示为:
$${K_z}{\nabla ^2}p + {\rho _f}{{{\partial ^2}\varepsilon } \over {\partial {t^2}}} - {\gamma _w}n\beta p{{\partial p} \over {\partial t}} = {{{\gamma _w}} \over {Kz}}{{\partial \varepsilon } \over {\partial t}}$$
(6)
$$\beta = {1 \over {{k_w}}} + {{1 - S} \over {{P_{w0}}}}$$
(7)
$$\varepsilon {\rm{ = }}{{\partial {u_x}} \over {\partial x}} + {{\partial {u_z}} \over {\partial z}}$$
(8)
简化了的u-p形式的控制方程为:
$${\sigma _{ij,j}} + \rho {g_i} = \rho {{\ddot u}_i}$$
(9)
\[-{{P}_{,i}}+{{\rho }_{f}}{{g}_{i}}={{\rho }_{f}}{{{\ddot{u}}}_{i}}+\frac{{{\rho }_{f}}{{g}_{i}}}{{{k}_{i}}}\dot{\bar{w}}\]
(10)
\[{{{\dot{u}}}_{i,i}}+{{{\dot{\bar{w}}}}_{i,i}}=-\frac{n}{{{K}_{f}}}\dot{p}\]
(11)
海床表面忽略水的粘性和摩擦力,即垂直方向应力${{\sigma '}_{zz}}$和剪力σzx忽略不计(${{\sigma '}_{zz}}$=σzx=0)。海床表面超静孔压pf等于波压力pb。视海床底部刚性不透水,海床土骨架不存在水平位移与竖向位移,法向流量$({{\partial p} \over {\partial z}})$为0。在两侧土骨架边界处,水平位移及水平流量$({{\partial p} \over {\partial x}})$为0。
2 椭圆余弦波作用下海床动态响应模型验证及参数分析 2.1 模型验证为了验证本文所采用模型的正确性,将本文所用的椭圆余弦波模型退化到使用Stokes波,将结果与解析解进行比较,如图3,实线为当前模型,点为解析解;输入的数据为北海波浪设计值:T=15 s,d=70 m,L=311.59 m,h=25 m,G=107 N/m3,μs=1/3,n=0.3,S=1.0,Kz=0.01 m/s。可见本文模型与解析解吻合较好,验证了本文结果的有效性。
2.2 参数分析本部分研究了饱和度,渗透系数对于海床孔压,竖向有效应力和竖向位移沿海床深度分布的影响。标准参数见表1、2,为了保证结果的合理性,选取第二周期末的时间点进行分析,使用单位化后的孔压p/P0,竖向有效应力${{\sigma '}_{zz}}$/P0,竖向位移uz2G/P0进行描述。为研究单独某项参数对于动态响应结果的影响,不考虑各参数之间的联动情况。
厚度 s/m | 弹性模量 E/(N· m-2) | 长度 l/m | 孔隙率 n | 密度 ρs/(kg· m-3) | 渗透系数Kz/(m· s-1) | 泊松比 μs | 饱和度S |
60 | 4.8 × 107 | 100 | 0.4 | 2 650 | 0.01 | 0.25 | 0.98 |
饱和度是海床土体的重要参数,图4~6描述了不同的饱和度(S=1,0.98,0.96,0.94)下孔压,竖向有效应力以及竖向位移沿海床深度的分布。图7为孔压分布云图。
如图4所示,孔压随着深度增加逐渐减小,减小过程中梯度逐渐增大,随着饱和度的增大,孔压减小得更慢。
图5表明竖向有效应力呈现先减小后增大的趋势,随着饱和度减小,竖向有效应力峰值增大。图6中可见竖向位移随着海床深度增大而逐渐的减小直至出现反向位移,竖向位移随着饱和度的减小而增大。
图7描述了整体海床中的孔压分布,可见随着饱和度的减小,孔压沿深度方向减小更快。
2.2.2 渗透系数Kz的影响图8~10给出了不同渗透系数下沿海床深度方向的孔压、竖向有效应力、竖向位移。如图8所示,孔压随着海床深度的增加而逐渐减小,渗透系数越大,减小得越慢。由图9,竖向有效应力沿深度方向呈现先减小后增大的趋势,渗透系数越大,竖向有效应力峰值越小。图10给出了位移随海床深度的变化,可见,渗透系数越小,海床表面产生的竖向位移越大。
2.2.3 不同波浪理论的影响
对线性波浪理论和椭圆余弦波理论下孔压的变化做出比较。由图11可见,不同的波浪理论对于海床表面的孔压影响较大,孔压随深度减小趋势基本相同,但由于两种波浪作用下的海床表面所受波动水压力差异较大,导致靠近海床表面部分,椭圆波作用下海床的孔压明显大于线性波作用下的孔压,因此在计算浅海地区波浪作用下的海床动态响应时不能简单地使用线性波理论进行简化。
3 液化判断准则及液化参数分析 3.1 液化判断准则
在波浪荷载持续作用下孔隙水压力上升,有效应力下降,海床的一部分可能变得不稳定而发生液化,本文采用由Zen[10]提出的基于超静孔隙水压力的液化判断准则,认为当土层中某一点处上层土体骨架的重量小于超静孔隙水压力时海床发生液化,即:
-(γs-γw)z+(pb(x,t)-p(x,z,t))≤0
(12)
使用液化准则分析不同饱和度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 。
3.2.2 渗透系数Kz的影响针对不同的渗透系数(Kz= × 10-3 m/s,1 × 10-4 m/s,1 × 10-5 m/s),图13给出了其液化区域。可见随着渗透系数的减小,液化深度和液化横向范围逐渐增大。
Kz=1 × 10-3 m/s时土体不发生液化,从Kz=1 × 10-4 m/s到Kz=1 × 10-5 m/s,液化深度从0.6 m变化到1 m,液化横向范围从20 m增大到27 m。这一现象具有较大的工程意义。表明近海中表明近海中海床土体的渗透系数较大时,液化的可能性将降低。在工程实践中可以采取铺设渗透系数较大的粗粒料上覆层作为保护管道的措施。
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 fluid-saturated 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. |