0.   引言

在海洋结构物设计中，波浪引起的海床地基失稳是一个需要引起足够重视的问题。波浪引起的动水压力会在海床中产生循环剪切作用，引起土体中超孔隙水压力的动态累积，使海床土体的抗剪强度大幅度降低甚至发生液化，从而引起显著的变形和失稳。因此，对于近海结构物有必要针对波浪长期作用引起的海床地基液化问题进行深入研究^[1-2]。

目前国内外学者对波浪作用导致的海床液化问题已进行了一些相关研究，普遍认为波浪引起海床动力反应的机制主要取决于海床中超孔隙水压力的产生方式：①超孔隙水压力的循环累积，主要是由土体的残余体积变形引起的，类似于地震作用引起的饱和土体残余超孔压增长模式^[3-4]；②可恢复的振荡超孔隙水压力，沿海床深度存在着幅值衰减和相位迟滞^[5-6]，这种超孔压的本质是由土体的可恢复的压缩体积变形引起的。Zhang等^[7]提出了一种三维的多孔弹塑性模型能同时模拟振荡超孔压及累积超孔压的特性，王小雯等^[8]基于砂土液化大变形本构理论提出了随机波浪作用下饱和海床土体超孔压瞬态变化的弹塑性动力分析方法，Liao等^[9]在现有二维模型的基础上建立了一个三维多孔弹塑性模型模拟砂质海床中超孔压的累积及消散过程。然而，现有大多数研究均未能完全考虑波浪作用下海床土体产生复杂的流固耦合效应，特别是累积超孔压导致的土骨架循环软化效应。

在波浪传播对海床表面的循环波压力的作用下，海床土体内的有效应力场和超孔压场都随着时间和空间的位置而变化（图 1），所产生的循环剪切作用同时含有正应力差和剪应力两者的影响^[10-11]。与地震作用下水平场地主要受到纯剪切作用不同^[12]，不考虑正应力差和剪应力连续变化下土体的非线性特性显然无法真实反映波浪作用下海床的动力特性^[13]。

图  1  波浪作用下海床的受力状态

Figure  1.  Stress states of seabed under wave action

下载: 全尺寸图片 幻灯片

基于此，本文将描述土体动力非线性特性的Davidenkov骨架曲线与采用“扩展Masing”法则构造的土体应力–应变滞回曲线拓展至三维应力空间，基于Biot动力固结理论建立了一个可描述流固耦合效应的土体非线性动力本构模型。基于FLAC^3D计算平台实现了该本构模型的二次开发，并验证了有效性。

1.   饱和土体动力流固耦合本构模型

1.1   饱和土体动力流固耦合效应模拟

一般情况下，饱和土体为两相介质，由土骨架与孔隙流体组成。Biot动力固结理论考虑了可变形的多孔隙土骨架和孔隙内流体之间的相互耦合机制。假定土骨架和孔隙中的流体是可变形的。孔隙流体服从各向同性的达西定律，土骨架各向均质同性，控制方程为

式中：u为土骨架的位移矢量；$\varepsilon$为力学效应下的体积应变；M为比奥模量；$p$为超孔隙水压力，由振荡超孔压$\widetilde p$与累积超孔压$\overline p$两部分组成；ρ为土的密度；${\gamma _{\text{w}}}$为孔隙水重度；G为土体剪切模量；$\nu$为泊松比。对完全饱和土体，当认为土颗粒不可压缩时$M = {k_{\text{w}}}/n$，${k_{\text{w}}}$为流体体积模量，n为孔隙率。

式（2）中饱和土体的体应变${\varepsilon _{jj}}$有两个分量组成：

式中：$\varepsilon _{jj}^{{\text{re}}}$为可恢复的压缩体应变分量；$\varepsilon _{jj}^{{\text{ir}}}$为循环剪切产生的残余体应变分量。现有Masing型非线性动力本构模型无法描述该残余体应变分量，而求解循环荷载作用下累积超孔压的关键在于残余体应变分量的计算。本文在体积相容方程（2）中引入一个新的源项${F^*}$用以描述累积超孔压，定义为残余体应变增量，参照文献[14]中的表达式为

式中：C₁，C₂为常量，与相对密度相关的经验计算公式为${C_1} = 7600{({D_{\text{r}}})^{ - 2.5}}$，${C_2} = {C_1}/0.4$；$\gamma$为剪应变；${\gamma _{{\text{th}}}}$为弹性应变阈值，可通过循环排水剪切试验来确定。只有当循环荷载作用下$\gamma$大于给定的阈值${\gamma _{{\text{th}}}}$，残余体应变增量的计算才会被触发，否则为零。

1.2   土体循环软化模拟

采用Davidenkov本构模型描述骨架曲线的非线性，结合扩展Masing法则构造动应力-应变滞回关系曲线^[15-16]。初始加载时，骨架曲线表示为

式中：

t时刻的时变剪切模量可表示为

依据扩展Masing法则，当应变反转后，后续加载过程的剪应力-应变滞回曲线表达式为

t时刻的时变剪切模量可表示为

式中：τ为剪应力；G_max为最大动剪切模量；A，B，${\gamma _0}$为与土性有关的试验参数，通过拟合共振柱或循环三轴试验获得的G/G_{max, 0} -$\gamma$曲线得到^[14]，${\tau _{{\text{rev}}}}$，${\gamma _{{\text{rev}}}}$分别为加卸载转折点处的剪应力与剪应变。

对于饱和土体，骨架曲线的刚度衰退可认为是由于累积超孔压的发展引起的，这是控制土骨架变形与孔隙水流动耦合的关键。因此在应变反转点执行以下算法，即通过更新Davidenkov骨干曲线$f(\gamma )$中的G_max来描述土骨架的循环软化：

式中：G_{max, 0}为初始最大动剪切模量；${\sigma '_{{\text{p0}}}}$为土的初始平均有效应力；${r_{\text{u}}}$=$\bar p$/${\sigma '_{{\text{p0}}}}$为累积超孔压比。其中$G_{\max }^{}$计算过程中涉及的$\bar p$由每循环半周下的增量叠加而得，故$G_{\max }^{}$在每次应变反转时进行更新，设定土体液化后的残余剪切模量G_max=0.01G_{max, 0}。图 2为本文DavByrne模型描述的应力-应变关系曲线示意图。

图  2  DavByrne模型描述的应力-应变关系曲线示意图

Figure  2.  Shear stress-strain curves described by DavByrne model

下载: 全尺寸图片 幻灯片

Hashash等^[17]指出海相沉积土层的动剪切模量比$G/{G}_{\mathrm{max}，0}^{}\text{-}\gamma$曲线依赖于土体原位埋深，须对有效围压有关的${G_{\max , 0}}$和${\gamma _0}$进行修正：

式中：$G_{_{\max , 0}}^{{\text{ref}}} = \rho {V_{\text{s}}}^2$，${V_{\text{s}}}$为剪切波速，${\sigma '_{{\text{ref}}}}$为参考上覆土压力，取100 kPa；a₁，a₂为土性相关的试验参数，可通过室内循环剪切试验拟合获得^[18]。在相同的剪应变幅值下，初始平均有效应力${\sigma '_{{\text{p0}}}}$越大，动剪切模量比曲线的衰减越慢。

1.3   三维拓展

波浪作用下海床土体的正应力和剪应力连续变化，因此需要在三维应力空间中进行分析。由前述内容可知，累积超孔压的源项F^*与时变剪切模量${G^{\text{t}}}$的计算是本文提出的考虑流固耦合效应的非线性动力本构模型的核心，而两者均是关于剪应变$\gamma$的函数。显然，选择合适的等效剪应变是构造本文三维本构模型的关键。参照文献[16]，定义等效剪应变${\gamma _{{\text{eq}}}}$代替上述表达式中的剪应变$\gamma$，可简单有效地拓展至三维空间：

式中：${\boldsymbol{J'}_2}$为偏应变张量的第二不变量，定义新的应变张量$\boldsymbol{e}_{ij}^{{\text{New}}}$为以反转点处的应变张量$\boldsymbol{e}_{ij, {\text{rev}}}^{}$为起始点，以反转后当前时刻的应变张量$\boldsymbol{e}_{ij}^{}$为终点的张量。

t时刻应力-应变滞回曲线段的时变体积模量${K^{\text{t}}}$可由下式确定：

由此t时刻的三维空间增量形式的应力张量${\text{d}}\boldsymbol{\sigma }_{ij}^{\text{t}}$按下式更新：

式中：$C_{ijkl}^{\text{t}}{\text{ = }} \left[ {\begin{array}{*{20}{c}} {({K^{\text{t}}}{\text{ + }}\frac{4}{3}{G^{\text{t}}})}&{({K^{\text{t}}} - \frac{2}{3}{G^{\text{t}}})}&{({K^{\text{t}}} - \frac{2}{3}{G^{\text{t}}})}&0&0&0 \\ {({K^{\text{t}}} - \frac{2}{3}{G^{\text{t}}})}&{({K^{\text{t}}}{\text{ + }}\frac{4}{3}{G^{\text{t}}})}&{({K^{\text{t}}} - \frac{2}{3}{G^{\text{t}}})}&0&0&0 \\ {({K^{\text{t}}} - \frac{2}{3}{G^{\text{t}}})}&{({K^{\text{t}}} - \frac{2}{3}{G^{\text{t}}})}&{({K^{\text{t}}}{\text{ + }}\frac{4}{3}{G^{\text{t}}})}&0&0&0 \\ 0&0&0&{{G^{\text{t}}}}&0&0 \\ 0&0&0&0&{{G^{\text{t}}}}&0 \\ 0&0&0&0&0&{{G^{\text{t}}}} \end{array}} \right] 。$

2.   程序流程与编程概要

FLAC^3D提供了内置FISH语言或C⁺⁺语言进行本构模型的二次开发程序编译。由于C⁺⁺语言面向对象，具有高兼容性、与硬件密切交互、且计算效率是FISH语言的（10~100）倍等优点^[19]，因此，在Visual Studio 2010环境中采用C⁺⁺语言进行本构模型的二次开发编译。本文的开发工作是编写头文件DavByrne.h和源文件DavByrne.cpp，编译生成DavByrne.dll动态链接库文件，主程序自动加载进行自定义本构模型的调用。

FLAC^3D采用显式有限差分法，为了避免数值计算中沙漏问题，利用混合离散技术将一个六面体单元（zone）划分成10个楔形子单元（subzone）。自定义本构模型要对每个单元进行材料参数、状态变量以及单元应变增量进行赋值。单元材料参数的赋予主要利用成员函数setProperty()与FLAC^3D程序命令Zone property实现，利用成员函数Initialize()执行CYCLE命令时进行一次状态变量初始化，单元应变增量由程序自行计算^[19]。自定义本构的主体是成员函数Run()，其功能是根据读取程序提供的应变增量，在每个子单元中调用自定义本构模型计算应力增量。

上述第一节理论与算法的重点是动剪切模量的更新与累积超孔压的计算，通过执行DavByrne: : Run()函数实现：

（1）根据初始应力状态修正G_{max, 0}与${\gamma _0}$。

（2）根据主程序提供的应变张量增量计算等效剪应变，依据加卸载准则判断是否发生应变反转。若发生应变反转，则根据式（4）计算$\varepsilon _{jj}^{{\text{ir}}}$，更新当前r_u与G_max；若不发生应变反转，则直接进入步骤（3）。

（3）更新${G^{\text{t}}}$与${K^{\text{t}}}$。

（4）更新$C_{ijkl}^{\text{t}}$，计算并返回主程序当前应力值。

重复上述步骤直至计算结束。具体的DavByrne模型计算流程如图 3所示。

图  3  考虑流固耦合效应的非线性动力模型程序流程图

Figure  3.  Program flowchart of nonlinear dynamic constitutive model considering pore-water-soil-skeleton coupling effects

下载: 全尺寸图片 幻灯片

在DavByrne模型开发流程中有以下3点须予以注意：

（1）应变反转问题的处理遵循原程序的方法，此处不再赘述，具体参考文献[19]。

（2）考虑到初始骨架曲线段的时变剪切模量与后续加卸载滞回曲线段的时变剪切模量计算不同，需要在程序中设置判断条件。为此，设置一个可以记录应变反转次数的状态变量，当值为零时，按初始骨架曲线计算，大于零时，按后续加卸载滞回曲线计算。

（3）当波浪较小时，可近似地认为波浪仅引起海床发生瞬时弹性应变和振荡超孔压，不会导致残余变形与累积超孔压的发展^[20]。为此，在程序中设置弹性应变阈值作为判断标准，即动荷载引起的应变小于该阈值时不执行超孔压累积的计算程序。

3.   模型验证

3.1   与模型理论值对比

在进行程序验证时，可仅采用一个单元进行计算。为了验证程序编制的可靠性，检验程序能否实现非线性滞回特性及超孔压增长规律的模拟，用一个1 m×1 m×1 m的立方体土单元模拟循环剪切试验，土体重度γ_sat=19.8 kN/m³，泊松比ν=0.3。设置100 kPa静孔压、200 kPa围压的初始条件，施加幅值为0.5 %的循环剪应变，模型及加载方式如图 4所示，模型参数见表 1。

图  4  数值模型及加载路径形式

Figure  4.  Numerical model and curve of loading path

下载: 全尺寸图片 幻灯片

表  1  模拟循环剪切试验模型参数

Table  1.  Model parameters of simulated cyclic shear tests

剪切模量G/MPa	Davidenkov模型					Byrne模型
	A	B	${\gamma _0}$/10-4	${\sigma '_{{\text{ref}}}}$/kPa		C1	C2	${\gamma _{{\text{th}}}}$/%
50	1.02	0.35	4	100		0.43	0.93	0.02

下载: 导出CSV 

| 显示表格

计算得到的应力-应变关系发展曲线如图 5所示。由图可知，土单元的骨架曲线随循环剪应变不断衰退，表现出了明显的非线性特征，结合图 6可知土单元在第80000分析步时发生液化，此时剪切模量并未衰退至0，而是保留了一定的残余剪切模量。此外，图 6中的超孔压增长仅表现出了循环累积特性，主要由循环剪切导致的残余体应变引起的，土骨架在纯剪切状态下未产生可恢复的压缩体应变。

图  5  模型计算的应力-应变关系曲线

Figure  5.  The calculated stress-strain curve

下载: 全尺寸图片 幻灯片

图  6  模型计算的超孔压变化规律

Figure  6.  Variation law of calculated excess pore pressure

下载: 全尺寸图片 幻灯片

为了验证程序的计算精度，同样以该单元模型为例，围压分别为200，400，600 kPa。施加剪应变的幅值为5%，计算得到的应力-应变关系发展曲线如图 7所示。由图 7可以看出，随着围压的增加，单元的初始剪切模量随之增大。在给定加载路径下，本文编制的本构模型程序计算值与理论值吻合较好，最大误差控制在1%以内。

图  7  计算与理论的应力-应变关系曲线对比

Figure  7.  Comparison of stress-strain relationship curvess between calculation and theory

下载: 全尺寸图片 幻灯片

3.2   与土单元循环扭剪试验结果对比

为说明DavByrne模型能有效模拟波浪作用下土体的动力反应，通过空心圆柱扭剪仪对试样施加波浪作用下土体的应力路径。试验粉土取自南通通州湾，土粒相对质量密度G_s=2.70，最大孔隙比e_max=1.29，最小孔隙比e_min=0.63，颗粒分析试验曲线如图 8所示，其中粉粒含量为50.7%，黏粒含量为1.3 %，细砂含量40.7 %，塑性指数I_p=3。

图  8  粉土级配曲线

Figure  8.  Grain-size distribution curve of silt

下载: 全尺寸图片 幻灯片

试样为空心圆柱试样，外半径r₀=50 mm，内半径r_i=30 mm，高度H=200 mm。试验设计相对密度D_r=50%，通过空中砂雨法制样。制样完成后，联合利用通CO₂、通无气水与施加反压等方法进行试样的饱和，确保试样固结前孔压系数B≥0.97。试样在有效围压为100 kPa的作用下进行均等固结。固结完成后，选取东海海域50 a一遇重现期的波浪要素，在不排水条件下，对试样施加频率为0.1 Hz，幅值为26 kPa的轴向动应力以及频率为0.1 Hz，幅值为13 kPa的扭剪应力来模拟波浪荷载，试样中一点的应力路径可由下式表示：

式中：${\sigma _x}$，${\sigma _z}$，${\tau _{xz}}$分别为水平正应力，竖向正应力及剪应力，${p_0}$为应力路径中应力圆的半径。

式（16）说明波浪荷载在土单元引起的动正应力和动剪应力呈周期变化，但动偏差应力不变，在坐标系$({\sigma _x} - {\sigma _z})/2 - {\tau _{xz}}$中为圆形。

数值计算可通过模拟单元的应力状态来与试验数据对比，单元初始阶段施加围压200 kPa，静孔压100 kPa，设置排水边界求解得到初始应力状态，加载时边界不排水。根据式（16），同时施加轴向${\sigma _{\text{d}}}$=26cos(ωt) kPa和切向${\tau _{\text{d}}}$=13sin(ωt) kPa的动应力，模型参数见表 1。

图 9给出了试验得到的应力路径与DavByrne模型计算的应力路径二者的对比，从图 9中可以看出DavByrne模型能准确地反映波浪作用下土体的应力路径。

图  9  试验与模拟单元的应力路径对比

Figure  9.  Comparison of stress paths between tests and calculation

下载: 全尺寸图片 幻灯片

图 10为本文DavByrne本构模型和FLAC^3D内置的Finn本构模型预测的超孔压与试验测得的超孔压对比结果。图中可知两种动力本构模型均能反映波浪荷载作用下土体超孔压的振荡及累积特性。试验中土体在120 s后液化，超孔压达到最大值。显然，本文提出的DavByrne动力本构模型能更好地反应波浪荷载作用下土体的超孔压增长特性，与试验实测数据较为吻合。然而，内置的Finn本构模型由于未考虑土骨架的循环软化特性，在超孔压上升至40 kPa时便达到稳定，并未出现液化现象，无法准确反映出波浪荷载作用下土体的动力反应。

图  10  波浪荷载作用下试验与DavByrne模型和Finn模型的超孔压结果对比

Figure  10.  Comparison of results among tests, DavByrne model and Finn model under wave loading

下载: 全尺寸图片 幻灯片

3.3   与土单元循环三轴试验结果对比

为验证DavByrne模型能有效模拟地震作用下土体的动力反应，对南京饱和细砂（相对密度D_r=45%）的不排水循环三轴试验进行模拟。建立一个底部固定竖直方向自由度的土单元，在单元的顶部节点施加一个频率为1 Hz、幅值为ε_a=0.15%的正弦轴向应变时程D_z(t)=-ε_a sin(2πt)，且边界设为不排水，初始围压${\sigma '_{x{\text{, 0}}}} = {\sigma '_{z{\text{, 0}}}}$=100 kPa，图 11给出了模型示意图，土体重度γ_sat=18.3 kN/m³，泊松比ν=0.25，具体模型参数如表 2所示。

图  11  地震作用下试验与DavByrne模型的超孔压比结果对比

Figure  11.  Comparison of excess pore pressure ratios between tests and DavByrne model under seismic loading

下载: 全尺寸图片 幻灯片

表  2  模拟循环三轴试验模型参数

Table  2.  Model parameters of simulated cyclic triaxial tests

剪切模量G/MPa	Davidenkov模型					Byrne模型
	A	B	${\gamma _0}$/10-4	${\sigma '_{{\text{ref}}}}$/kPa		C1	C2	${\gamma _{{\text{th}}}}$/%
53	1.02	0.43	4.1	100		0.55	1.38	0.02

下载: 导出CSV 

| 显示表格

本文DavByrne模型超孔压比计算结果与试验值对比如图 11所示。可以看出超孔压比的发展呈现出振荡与累积特性，与试验结果具有良好的一致性。在前20周循环荷载下超孔压快速增长，随后累积速度放缓且振荡幅值也随之减小，在第40周时，土体发生液化，超孔压无法累积，仅呈现出振荡特性。

4.   结论

本文建立了一个可描述流固耦合效应的土体非线性动力本构模型，基于FLAC^3D计算平台实现了该本构模型的二次开发，并验证了模型的有效性。

（1）DavByrne模型在给定加载路径下表现出了良好的非线性滞回特性和循环软化特性，与模型理论解吻合较好，说明了该本构模型的合理性与较高的计算精度。

（2）与南通通州湾的重塑粉土的室内循环扭剪试验结果对比表明DavByrne模型能有效模拟饱和海洋土在波浪作用下的应力路径和超孔压增长规律。

（3）与FLAC^3D内置的Finn模型结果对比表明忽略由超孔压累积引起的土骨架循环软化效应将导致低估饱和土体的动力响应。