波浪诱发松散海床渐进式液化的数值分析

栾一晓

(1. 武汉理工大学，湖北 武汉 430070)

波浪诱发松散海床渐进式液化的数值分析

栾一晓1

(1. 武汉理工大学，湖北 武汉 430070)

近海区域广泛分布着第四纪新沉积的松散海洋土，波浪荷载作用下松散海床会发生液化进而对近海结构物的稳定性存在巨大威胁。本文采用中国科学院流体-结构-海床相互作用数值计算模型FSSI-CAS 2D，选用Pastor-Zienkiewicz-Mark Ⅲ(PZⅢ)弹塑性本构研究了波浪诱发的松散海床液化问题。分析了波浪荷载引起的松散海床内超孔隙水压力、有效应力以及应力角的时程变化特性，并预测了松散海床的渐进液化过程。计算结果表明，波浪荷载作用下松散海床内残余孔压会累积增长，海床表面最先发生液化，然后逐渐向下发展至液化最大深度。同时指出海床内超孔隙水压力的竖向分布特征和应力角的变化时程均可以作为判断海床液化的间接参数。最后，通过应力状态分析，讨论了海床渐进式液化的发展过程和趋势。

渐进式液化；松散海床；波浪荷载；FSSI-CAS 2D；PZⅢ本构模型

1 引言

波浪荷载作用下海洋结构物的稳定性是工程设计中需要重点考虑的问题之一。近海广泛存在着松散的第四纪新沉积海洋软土，当遭遇较大波浪荷载时，土体颗粒会伴随着孔隙水的排出而重新排列，海床变得更为密实，过程中超孔隙水压力会逐步累积，海床土体的塑性变形会累积增大，引起海床液化进而导致地基承载力不足造成海洋结构物失稳破坏。因此海床在波浪荷载作用下的动态响应是评价服役期内离岸结构物稳定性的重要因素，而目前对波浪荷载作用下海床内超孔隙水压力的增长模式以及海床的液化机理尚揭示不足，迫切需要针对该问题开展研究工作。

针对波浪-海床相互作用问题已有大量研究成果，最先是基于Biot理论提出的解析解分析模型[1]，但该分析方法的局限在于假定海床是十分密实的土体，在波浪荷载作用下主要是发生弹性变形。Lee借助耦合求解方法进一步考虑海床表面海水和海床内孔隙水的物质交换及压力传递连续条件分析了波浪和海床的相互作用[2]，但仅限于自由海床的分析，尚无法用于分析上覆有结构物工况时海床的动态响应。

除解析解外，数值计算分析在研究波浪作用下海床的动力响应方面取得了重要突破。早期的数值模拟仍采用弹性本构方程，同时假定海床处于密实状态[3]，实际上海洋工程结构物的设计中应更多考虑新沉积的松散海洋土，采用更为先进的本构模型来描述海洋土的非线性和弹塑性。在已有的数值研究成果中，主要分为两类研究方法。第一类方法主要是通过循环荷载周次来表征海床发生液化时土所承受的剪应力，如Seed等[4-5]通过土工实验数据建立了剪切应力比和破坏周数的关系。但该类方法中的剪应力大都是基于多孔弹性理论求解，所求得的剪应力幅值是恒定不变的，无法反映松散海床中随着孔隙水压力的累积剪应力逐渐减小直至伴随着土体的液化最终为零的物理现象，从而严重高估了波浪的动荷载效应。第二类方法是应用弹塑性本构模型来考察松散海床的非线性动力响应。Sassa等[6]基于两层流和移动边界理论，提出孔压预测模型的主要思路是认为循环荷载作用下塑性体应变的变化率依赖于剪切应力比与最终应变状态的乘积，其预测结果与试验结果吻合很好，但尚无法得到松散海床内不同位置的有效应力。2006年，Dunn等[7]应用Pastor等[8]和Zienkiewicz等[9]提出的Pastor-Zienkiewicz Mark-Ⅲ (PZⅢ)模型，研究了渐进波作用下埋置管线周围的松散海洋土的液化情况。Jeng和Ou[10-11]将PZⅢ模型从二维扩展到了三维，但已有研究成果还不足以深入地揭示松散海床土体在波浪荷载作用下的液化规律。

本研究采用FSSI-CAS 2D[12]数值计算模型对新沉积的松散海床土体在波浪荷载作用下的液化规律做进一步的研究和探讨。数值模型中采用PZⅢ弹塑性本构模型，计算过程中孔隙率及相应的渗透系数随着波浪荷载的动态作用而发生的实时变化，通过实时更新刚度矩阵来反映海洋土在循环荷载作用下的非线性动力特性。研究结果表明，基于FSSI-CAS 2D模型的数值分析可以相对有效地反映松散海床在波浪荷载作用下的动力响应特征，通过对松散海床内超孔隙水压力、有效应力以及应力角的时程变化特性的分析，提出了判断海床液化的间接依据，并预测了波浪诱发海床渐进式液化的发展过程。

2 数值模型和土体本构

应用Zienkiewicz等[9]、Biot[1]和Chan[13]提出的u-p格式动力Biot方程作为控制方程，模拟波浪荷载作用下多孔介质海床的动态响应：

(1)

(2)

k2ps-γwnβ+kρf=γw，

(3)

(4)

式中，Sr为海床土体饱和度；pw0为静水压力的幅值；Kf为孔隙水的体积模量，通常取值为2.24×109N/m2；孔隙液体压缩性β反映了海床土体的不饱和程度。实际上，在近海岸区域海床土体的饱和度一般大于90%，接近于完全饱和。

应用广义Newmark[13]法的隐式求解公式(1)至式(3)。对于波浪-结构-海床之间的相互作用问题，Ye[14]建立了耦合的数值模型FSSI-CAS 2D，模型中使用Hsu等[15]提出的Volume Average Reynold Average Navier Stokes(VARANS)方程控制波浪运动和孔隙介质中流体的流动方式，而结构物和海床的动态响应由广义Biot动力方程控制，Ye等[14,16]和Zienkiewicz等[17]的研究成果中提供了将Biot动力方程和VARANS方程进行耦合求解的算法。

以往研究中通常基于变形较小的假设而在计算中保持孔隙比e和渗透系数k为常数。本研究中依据大变形理论提出的方程en+1=(1+en)exp(ΔpQ+Δεvs)-1考虑了海床土体的孔隙比及相应的渗透系数随土体颗粒的变形和重组的变化，其中n表示计算的时间步，Δp为孔隙水压力的增量，Δεvs为海洋土的体应变增量，Q=1/β为孔隙水的压缩性。相应的，海洋土的渗透系数k可以表述为k=Cfe31+e，其中，Cf为经验系数，取决于Cf=k01+e0e30[18]，式中e0为初始孔隙比。波浪荷载作用下松散海床会发生较大变形，海床表面的几何形状会发生显著变化，因此计算中海床表面的静水压力和动水压力将随着海床表面的变形而动态改变，保证了不同时间步内计算域的力边界条件始终实时更新，从而实现耦合求解。

在FSSI-CAS 2D数值模型中采用PZⅢ本构模型[8]描述海洋土的动力特性。PZⅢ模型作为一种弹塑性本构模型，其可靠性已经被一系列单调加载试验、循环荷载试验、离心模型试验及相关工况的数值模拟结果[19]所证实。

3 验证

2012年Ye[14]通过对比已有的数模和物模研究成果，验证了半耦合数值计算模型FSSI-CAS 2D在预测自由海床和防波堤地基在波浪荷载作用下的动力响应方面的有效性和可靠性，如Hsu和Jeng的解析解[20]，规则波和椭圆余弦波的水槽试验[21]，驻波[22]，半潜式防波堤[23]以及组合式防波堤[24]。随后Ye等[25-26]进一步通过与波浪槽试验[27]和离心模型试验[19]的模拟，验证了FSSI-CAS 2D计算模型在求解波浪-松散海床相互作用方面的有效性。

4 模型边界条件和水动力荷载

计算模型为一个400 m长、20 m厚的平坦海床，其中，设置该区域的左下角点为原点。水平方向的网格尺寸为1 m，竖直方向网格尺寸为0.5 m，总共生成12 000个4节点单元。边界条件设置如下：第一，底部边界不允许排水；第二，两个水平边界的水平方向不允许有变形；第三，在海床表面施加静水压力。在每一个时间步，激活施加在海床的静水压力作为施加在海床表面的边界数值，每一步更新公式为ps=ρgd0+ρgsv，式中，d0为初始水深，sv为波浪荷载作用下海床表面节点的竖向残余位移和瞬间位移之和；第四，波流引起的动水压力也作用在海床上，其三阶表达式为[28]：

cos3(λx-ωt)，

(5)

式中，H为波浪高度；λ=L/2π为波数；L为波长；ω=T/2π为角速度；U0为当前速度；d=d0+sv为即时水深。当不考虑流U0=0 m/s时，上述方程退化为经典的三阶非线性波。当在海床上施加静水压力和动水压力时，海床表面土体的有效应力确保为0。

PZⅢ模型参数可以根据近海土体的室内试验获得，本文中采用的是相对密度为60%的Nevada砂的参数[17]，列于表1。其中初始孔隙比e=0.333，饱和度为98%。相应地，海床土体初始的渗透系数为1.0×10-5m/s，初始水深为10 m，波高为1.5 m，波浪周期为8 s，流速U0=0.5 m/s。

表1 PZⅢ模型中松散海床土体参数[18]

5 计算结果

在近海环境中，静水压力会引起海床土的长期固结，但不引起超孔隙水压力。海床土的初始应力状态可由Ye[29]获得，以此开始进一步的动力分析。

5.1 有效应力和超孔隙水压力

图1给出了在松散海床3个典型深度处(z=18 m为在海床土体表面位置，z=10 m为海床模型的中间位置，以及z=2 m为海床底部附近位置)由波浪荷载引起的超孔压和有效应力的时程。超孔隙水压力由震荡孔压和累积孔压两部分组成，震荡孔压的幅值从海床表面向下逐渐减小，反之海床较深位置处残余孔压的幅值要远大于海床表面位置，达到峰值的时间也明显增多。波浪和海流荷载引起的累积孔压不会无限增大，达到峰值后会保持不变。

与累积孔压的趋势一致，松散海床土的有效应力从初始状态会逐渐减小，当有效应力达到0时，海床会失去抗剪强度而发生液化造成上部结构的倒塌。如图1所示，z=18 m处土体的平均有效应力在约150 s时达到0，发生部分液化。在300 s时，有效应力变为0，意味着z=18 m处的海床表面的土体全部发生液化。在有效应力下降的过程中，剪应力的幅值也减小。当z=18 m处海床土体发生全部液化时，由波浪和海流引起的剪应力也变化为0，原因是发生完全液化的土体的性质和流体近似，不能传递剪应力。FSSI-CAS 2D模型可以清楚地捕捉到这样的现象。同时也观察到z=10 m海床中部位置土的平均有效应力没有达到完全液化状态。然而，与初始应力的数值相比下降超过了90%，在约300 s时发生了部分液化。z=2 m接近海床底部位置直到800 s土体也没有发生液化。

图1 海床3个典型位置处波流荷载引发的超孔隙水压力(a)和平均主应力(b)时程Fig.1 Time history of wave & current-induced excess pore pressure (a) and mean principle stress (b) at 3 typical depth

5.2 应力角

应力角是研究松散海床在波流荷载作用下动态响应的另一个重要参数。应力角是基于摩尔库伦准测定义的：

(6)

图2为波流荷载作用下松散海床3个典型深度处的土体应力角变化时程。作为对比，也绘制出相同荷载条件下密实弹性海床土的应力角变化时程。对于密实的海床，应力角只是在初始值附近阶段性地变动。然而，松散海床的应力角在初始阶段逐渐减小，然后逐渐增大直到该位置土体发生部分或者完全液化，之后残余应力角基本保持稳定。如图2所示，在z=18 m处和在z=10 m处，两个位置残余应力角从增加到保持稳定的时间分别为150 s和300 s，时间与松散海床土体这两个位置发生部分或者全部液化的时间十分一致。这说明，残余应力角从增加到保持稳定的时间可以用来作为判断土体发生液化的一个间接依据。

5.3 渐进式液化

室内试验[19]和现场记录[30]都证明松散海床在波浪荷载作用下的液化存在两种机制。一种是瞬时液化，只发生在十分密实的土体中，另外一种是松散土体中由孔隙水压力累积引起的残余液化。本文研究的是发生在近海松散海床土体中的液化即为残余液化，定义Lpotential为液化势：

(7)

了解波浪荷载作用下松散海床中残余孔隙水压力沿竖向深度方向的分布特性是十分必要的。在图3中，给出了在x=200 m处不同时间土中累积孔压和震荡孔压沿竖向的分布。注意到在液化深度范围内震荡孔压的竖向分布是不规则的，超出液化深度后则保持同一趋势。因此震荡孔压的竖向分布可以作为液化判断的另一个间接依据。松散海床中累积孔压随时间持续增长，但不会超过界限值，即液化抵抗线(LRL)，一旦到达时该位置的土体即发生液化。根据式(7)，可以估计松散海床中液化可能深度Lpotential在不同时间(如t/T=10, 20, 45)的数值，如图4所示。图4显示，上层海床Lpotential比下层海床大一些，而且随着波浪的循环加载海床的Lpotential逐渐增加达到或超过0.86而发生液化。图3和图4都说明，累积超孔压达到液化抵抗线的时间随着海床深度增加而增加，这意味着波浪荷载引起松散海床液化是一个渐进的过程，最开始是在海床表面，然后逐渐向下扩展，如图4所示。在t=45T之后，液化深度达到11 m。

图2 波浪荷载诱发松散海床3个典型深度处土体应力角变化时程Fig.2 Time history of wave & current-induced stress angle at 3 typical depth

图3 在x=200 m处不同时间土中累积孔压(a)和震荡孔压(b)沿竖向的分布Fig.3 Vertical distribution of wave & current-induced residual pore pressure(a) and oscillatory pore pressure (b) on symetrical line x=200 m

图4 Lpotential在不同时间沿竖向的分布Fig.4 Vertical distribution of liquefaction potential Lpotential at typical time t/T=10, 20, 45

图5 在x=200 m处不同深度位置点的应力路径(单位：kPa)Fig.5 Stress paths of soil at a series of positions on x=200 m (unit: kPa)

图6 典型时间步下波浪荷载诱发松散海床液化区域分布Fig.6 Wave & current-induced residual liquefaction zone in loose seabed at typical time t/T=10, 25, 50, 100

应力路径是表征松散海床在波浪荷载作用下发生液化的另外一个重要参数。图5给出了在x=200 m处不同深度位置点的应力路径。可以发现，所有应力状态在初始时刻都在原始K0线上。在波浪荷载的循环作用下，松散海床中的孔隙水压力上升，有效应力下降。结果是，应力状态逐渐向零应力状态靠近，即为液化状态。在计算的结尾，在x=200 m处的一系列位置的应力状态达到了零应力状态，发生了部分或者完全液化。

如图6所示，通过Lpotential对整个海床的液化区域分布进行预测。可以看到，松散海床的液化区域随时间逐渐增大，受波浪荷载的影响，起初液化区域是波浪形状的，直到t/T=100，大部分松散海床发生液化。

6 结论

在近海区域内广泛分布着第四纪新沉积的松散海洋土，波浪荷载引起的松散海床液化对近海结构物存在巨大的威胁。研究中，采用FSSI-CAS 2D数值计算模型，选用Pastor-Zienkiewicz-Mark Ⅲ弹塑性本构模型研究了波浪荷载诱发的松散海床液化问题。求解了典型海况下波浪荷载作用下海床的动力响应，计算结果进一步证实了波浪荷载作用下伴随着残余孔压的累积松散海床的液化从海床表面开始，逐步向海床深部方向发展直至达到液化深度。结果表明可以根据震荡孔压的竖向分布、应力角的时程作为判断海床发生液化或预测液化深度的间接依据。最后，通过应力状态分析，讨论了海床渐进式液化的发展过程和趋势。研究结果可为近岸海洋结构物的设计提供技术参考。

[1] Biot M A. Theory of propagation of elastic waves in a fluid-saturated porous solid. Part 1: Low frequency range[J]. Acoustical Society of America, 1956, 28(2): 168-179.

[2] Lee T C, Tsai C P, Jeng D S. Ocean wave propagating over a porous seabed of finite thickness[J]. Ocean Engineering, 2002, 29(12): 1577-1601.

[3] Jeng D S. Wave-induced sea floor dynamics[J]. Applied Mechanics Reviews, 2003, 56(4): 407-429.

[4] Seed H B, Martin P O, Lysmer J. Pore-water pressure changes during soil liquefaction[J]. Journal of Geotechnical Engineering ASCE, 1976, 102(4): 323-346.

[5] 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.

[6] Sassa S, Sekiguchi H, Miyamoto J. Analysis of progressive liquefaction as a moving boundary problem[J]. Géotechnique, 2001, 51(10): 847-857.

[7] Dunn S L, Vun P L, Chan A H C, et al. Numerical modeling of wave-induced liquefaction around pipelines[J]. Journal of Waterway, Port, Coastaland Ocean Engineering, 2006, 132(4): 276-288.

[8] Pastor M, Zienkiewicz O C, Chan A H C. Generalized plasticity and the modelling of soil behaviour[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1990, 14(3): 151-190.

[9] Zienkiewicz O C, Chang C T, Bettess P. Drained, undrained, consolidating and dynamic behaviour assumptions in soils[J]. Géotechnique, 1980, 30(4): 385-395.

[10] Jeng D S, Ou J. 3D models for wave-induced pore pressures near breakwater heads[J]. Acta Mechanica, 2010, 215(1/4): 85-104.

[11] Ou J. Three-dimensional numerical modelling of interaction between soil and pore fluid[D]. Birmingham, UK: Universtity of Birmingham, 2009.

[12] 王良民, 叶剑红, 朱长歧. 近海欠密实砂质海床内波致渐进液化特征研究[J]. 岩土力学, 2015, 36(12): 3583-3588.

Wang Liangmin, Ye Jianhong, Zhu Changqi. Investigation on the wave-induced progressive liquefaction of offshore loosely deposited sandy seabed[J]. Rock and Soil Mechanics,2015, 36(12): 3583-3588.

[13] Chan A H C. A unified finite element solution to static and dynamic problems of geomechanics[D]. Swansea Wales: University of Wales, 1988.

[14] Ye J H. Numerical analysis of wave-seabed-breakwater interactions[D]. Dundee, UK: Universtity of Dundee, 2012.

[15] Hsu T J, Sakakiyama T, Liu P L F. A numerical model for wave motions and turbulence flows in front of a composite breakwater[J]. Coastal Engineering, 2002, 46(1):25-50.

[16] Ye J H, Jeng D S, Wang R, et al. Validation of a 2D semi-coupled numerical model for Fluid-Structures-Seabed Interaction[J]. Journal of Fluids and Structures, 2013, 42(4): 333-357.

[17] Zienkiewicz O C, Chan A H C, Pastor M, et al. Computational Geomechanics with Special Reference to Earthquake Engineering[M]. Chichester: John Wiley & Sons Ltd, 1999.

[18] Miyamoto J, Sassa S, Sekiguchi H. Progressive solidification of a liquefied sand layer during continued wave loading[J]. Géotechnique, 2004, 54(10): 617-629.

[19] Sassa S, Sekiguchi H. Wave-induced liquefaction of beds of sand in a centrifuge[J]. Géotechnique, 1999, 49(5): 621-638.

[20] Hsu J R, 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.

[21] Lu H B. The research on pore water pressure response to waves in sandy seabed[D]. Changsha: Changsha University of Science & Technology, 2005.

[22] Tsai C P, Lee T L. Standing wave induced pore pressure in a porous seabed[J]. Ocean Engineering, 1995, 22(6): 505-517.

[23] Mizutani N, Mostarfa A, Iwata K. Nonliear regular wave, submerged breakwater and seabed dynamic interaction[J]. Coastal Engineering, 1998, 33(2/3): 177-202.

[24] Mostafa A, Mizutani N, Iwata K. Nonlinear wave, composite breakwater, and seabed dynamicinteraction[J]. Journal of Waterway, Port, Coastal, and Ocean Engineering, ASCE, 1999, 25(2): 88-97.

[25] Ye J H, Jeng D S, Wang R, et al. Numerical simulation of wave-induced dynamic response of poro-elasto-plastic seabed foundation and composite breakwater[J]. Applied Mathematical Modeling, 2015, 39(1): 322-347.

[26] Ye J H, Wang G. Seismic dynamics of offshore breakwater on liquefiable seabed foundation[J]. Soil Dynamics and Earthquake Engineering, 2015, 76: 86-99.

[27] Teh T C, Palmer A C, Damgaard J S. Experimental study of marine pipelines on unstable and liquefied seabed[J]. Coastal Engineering, 2003, 50(1/2): 1-17.

[28] Ye J H, Jeng D S. Response of porous seabed to nature loadings: Waves and currents[J]. Journal of Engineering Mechanics, ASCE, 2012, 138(6): 601-613.

[29] Ye J H. Numerical modelling of consolidation of 2-D porous unsaturated seabed under a composite breakwater[J]. Mechanika, 2012, 18(4): 373-379.

[30] Sassa S, Takayama T, Mizutani M, et al. Field observations of the build-up and dissipation of residual porewater pressures in seabed sands under the passage of stormwaves[J]. Journal of Coastal Research, 2006, 39(12): 410-414.

[31] Ishihara K. Liquefaction and flow failure during earthquakes[J]. Géotechnique, 1993, 43(3): 351-451.

[32] Wu J, Kammaerer A M, Riemer M F, et al. Laboratory study of liquefaction triggering criteria[C]// Proceedings of 13th World Conference on Earthquake Engineering. Vancouver, British Columbia, Canada. 2004.

Wave-induced progressive liquefaction in loosely deposited seabed

Luan Yixiao1

(1.WuhanUniversityofTechnology,Wuhan430070,China)

Quaternary newly deposited loose seabed soil widely distributes in offshore area in the world. Wave-induced residual liquefaction in loose seabed floor brings great risk to the stability of offshore structures in extreme climate. In this study, wave & current-induced residual liquefaction in loose seabed floor has been investigated comprehensively adopting FSSI-CAS 2D incorporating Pastor-Zienkiewicz-Mark Ⅲ(PZⅢ) soil model which is a validated integrated numerical model. The time history of wave & current-induced pore pressure, effective stress, stress angle are discussed. The variation process of progressive liquefactionis illustrated in detail. The computational results confirm that the wave & current-induced liquefaction in loose seabed soil is progressively downward, initiating at seabed surface. Besides, it is found that vertical distribution of oscillatory pore pressure, and time history of stress angle could be taken as indirect indicator to judge the occurrence of wave-induced residual liquefaction. The developing process of the progressive liquefaction is analyzed by stress statement data of seabed soil.

progressive liquefaction; loose seabed floor; wave loading; FSSI-CAS 2D; PZⅢ soil model

10.3969/j.issn.0253-4193.2017.09.010

2017-06-12；

2017-08-05。

栾一晓(1996—)，女，山东省招远市人，主要从事士力学数值模拟。E-mail：lyx1996yx@126.com

TU431

A

0253-4193(2017)09-0101-09

栾一晓. 波浪诱发松散海床渐进式液化的数值分析[J]. 海洋学报, 2017, 39(9): 101-109,

Luan Yixiao. Wave-induced progressive liquefaction in loosely deposited seabed[J]. Haiyang Xuebao, 2017, 39(9): 101-109, doi:10.3969/j.issn.0253-4193.2017.09.010