地震长期以来被认为是诱发滑坡、土体液化、岩土结构破坏等地质灾害的主要原因之一^[1]。特别是松散饱和砂土地基在地震荷载作用下，由于土体内部超孔隙水压的累积而容易发生液化破坏现象。这一现象往往伴随着土体的大变形运动及地面沉降，从而成为工程结构破坏的一种主要触发性因素。例如，在1995年日本阪神地震中，地面由于砂土液化发生侧向大变形及沉降，致使建筑物倒塌或严重受损^[2]。在2018年印度尼西亚帕卢地震中，土体由于液化而产生了约600 m的漂移，致使其上建筑物全部被损毁^[3]。

动力有限元方法(finite element method，FEM)被广泛用于分析土体的地震液化破坏问题^[4-5]。基于有限元法的地震分析可以考虑荷载的性质以及超孔压的产生和消散，从而能够评估土体的液化破坏风险。但是现有的有限元方法还不能在拉格朗日框架下恰当地计算分析土体液化后的大变形发展过程。当大变形发展时，由于网格纠缠、畸变等问题，经典有限元数值计算往往变得极不稳定^[6]。虽然可以通过重新划分网格在一定程度上克服这一缺点，但计算成本较高。近年来，物质点法(material point method，MPM)作为克服这一局限性的新型数值方法逐步流行。它结合了欧拉和拉格朗日框架在计算大变形效应方面的优点，消除了拉格朗日有限元法的网格畸变问题，可完全适用于岩土材料大变形问题的模拟分析^[7-11]。

在物质点法中，材料被视为一簇分散的物质点粒子，材料的相关属性信息统一附着在物质点上。计算网格作为背景网格，独立于真实材料的几何区域。物质点与背景网格节点通过映射建立联系。图 1为物质点法的主要计算流程。物质点法最初由Sulsky等^{[7, 12]}基于对胞内粒子法和流体隐式方法的扩展而提出的。目前，物质点法已在裂纹扩展、爆炸、岩土体大变形等方面得到应用。其中，在岩土体大变形方面，文献[13-14]采用单套物质点及两套物质点，用于研究饱和土在接触/冲击作用下的动力响应。文献[15]将MPM与FEM两种方法进行结合，编写了一个用于分析饱和土体变形问题的MPM-FEM耦合程序。文献[16]为研究饱和及非饱和条件下岩土体的大变形问题，开发了半隐式的物质点算法。但迄今为止，MPM还很少用于研究动力荷载作用下的土体液化问题。

本文首先基于多孔介质理论，推导了u-p格式下固液两相耦合物质点的理论算法公式，采用Fortran 95程序语言编制了固液两相耦合物质点法的三维显式程序MapleLeaf。其次，通过静水压力柱及太沙基一维固结理论模型的物质点法数值解与理论解的对比分析，检验了程序在分析固液两相介质相互作用问题方面的准确性。最后，引入可用于液化分析的弹塑性本构模型及动力算法模块，根据现有文献及现场振动台饱和地基的动力液化试验，验证了程序在分析饱和地基动力液化问题方面的有效性。

1 固液两相耦合物质点算法 1.1 饱和多孔介质理论

如图 2所示，假定岩土体为多相介质(固体骨架、液体及气体)组合体，在饱和状态下，忽略气体不计，所有孔隙均由液相填充，且各孔隙之间相互连通。设定土体单元总体积为V，孔隙率为n，固体骨架密度为ρ^sR，孔隙液体密度为ρ^wR，基于总体积的固、液两相的表观密度分别为ρ^s、ρ^w，则在满足质量守恒的条件下，满足：

$ \rho=\rho^{\mathrm{s}}+\rho^{\mathrm{w}} $

(1)

$ \rho^{\mathrm{s}}=(1-n) \rho^{\mathrm{sR}} $

(2)

$ \rho^{\mathrm{w}}=n \rho^{\mathrm{wR}} $

(3)

$ m^{\mathrm{s}}=\rho^{\mathrm{s}} V $

(4)

$ m^{\mathrm{w}}=\rho^{\mathrm{w}} V $

(5)

式中：ρ为混合体整体表观密度，为各相表观密度之和；m^s、m^w分别为固、液两相质量。

对饱和土体内各相进行分析，并考虑求解域为单一封闭系统，则基于多孔介质理论所得到的质量守恒方程与动量守恒方程为

$ \nabla \cdot \boldsymbol{\sigma}^{\mathtt{α}}+\rho^{\mathtt{α}} \boldsymbol{b}^{\mathtt{α}}-\rho^{\mathtt{α}} \boldsymbol{a}^{\mathtt{α}}=0 $

(6)

$ \frac{\mathrm{d} \rho^{\mathtt{α}}}{\mathrm{d} t}+\nabla \cdot\left(\rho^{\mathtt{α}} \boldsymbol{v}^{\mathtt{α}}\right)=0 $

(7)

式中：ρ^α、v^α、a^α、σ ^α和b^α分别为α相(固相s、液相w)的相对密度、速度、加速度、内力以及体力。

在控制方程建立过程中，忽略液相相对于固相的加速度值，根据式(6)、(7)，建立固、液两相的u-p格式方程：

$ \nabla \cdot\left(\boldsymbol{\sigma}^{\prime}-p \boldsymbol{I}\right)+\rho \boldsymbol{b}-\rho \ddot{\boldsymbol{u}}=0 $

(8)

$ \dot{p}=\frac{\mathrm{d} p}{\mathrm{~d} t}=\frac{K^{\mathrm{w}}}{n}\left(\frac{k}{\gamma^{\mathrm{wR}}} \nabla \cdot \boldsymbol{R}^{\mathrm{e}}+\nabla \cdot \boldsymbol{v}^{\mathrm{s}}\right) $

(9)

$ \boldsymbol{R}^{\mathrm{e}}=-\nabla p+\rho^{\mathrm{wR}} \boldsymbol{b}-\rho^{\mathrm{wR}} \boldsymbol{a}^{\mathrm{s}} $

(10)

式中：σ ′为有效应力，p为孔隙水压力，$\ddot{\boldsymbol{u}} $为固相加速度(即a^s)，γ^wR为液相真实重度(γ^wR=ρ^wRg)，R^e为渗流过程中拖拽力。

1.2 u-p格式方程积分弱形式与离散

考虑采用一套物质点同时表示固液两相介质，并对u-p格式方程进行积分弱形式与离散转化。对应饱和多孔介质的混合体密度可表示为

$ \rho(\boldsymbol{X}, t)=\sum\limits_{P=1}^{N_P} M_P \delta\left(\boldsymbol{X}-\boldsymbol{X}_P\right) $

(11)

式中：δ(·)为delta函数，X_P为物质点的位置矢量，M_P为物质点所在子域的质量(M_P=ρV_P)，V_P为对应子域体积，大写指标P指代物质点。

基于式(8)、(10)、(11)，将饱和多孔介质混合体动量方程转化为积分弱形式：

$ \begin{gathered} \int_{\varOmega} \delta u_i \rho \boldsymbol{a}_i^{\mathrm{s}} \mathrm{d} \varOmega+\int_{\varOmega} \delta u_{i, j} \boldsymbol{\sigma}_{i j} \mathrm{~d} \varOmega-\int_{\varOmega} \delta u_i \rho \boldsymbol{b}_i \mathrm{d} \varOmega- \\ \int_{\partial \varOmega} \delta u_i \boldsymbol{\sigma}_{i j} \boldsymbol{n}_j \mathrm{d} A=0 \end{gathered} $

(12)

$ \begin{gathered} \int_{\varOmega} \delta u_i \boldsymbol{R}_i^{\mathrm{e}} \mathrm{d} V-\int_{\varOmega} \delta u_{i, j} p \mathrm{d} V-\int_{\varOmega} \delta u_i \rho^{\mathrm{wR}} \boldsymbol{b}_i \mathrm{d} V+ \\ \int_{\varOmega} \delta u_i \rho^{\mathrm{wR}} \boldsymbol{a}_i^{\mathrm{s}} \mathrm{d} V+\int_{\partial \varOmega} \delta u_i p \boldsymbol{n}_i \mathrm{d} A=0 \end{gathered} $

(13)

式中：δu为试函数，Ω代表当前连续体的构型，∂Ω表示构形边界。

在固液两相耦合物质点算法中，形函数N_IP=N_I(X_P)为节点I的形函数在物质点P的值，大写指标I、J和K表示背景网格节点。物质点坐标可以通过形函数插值得到

$ \boldsymbol{X}_{i P}=\sum\limits_{I=1}^{N_n} \boldsymbol{X}_{i I} N_I\left(\boldsymbol{X}_P\right)=\sum\limits_{I=1}^{N_n} \boldsymbol{X}_{i I} N_{I P} $

(14)

式中：X_iI为背景网格节点坐标，I表示定义在网格节点I上的变量。

同理，物质点的位移、速度、加速度及试函数亦可通过背景网格节点插值获得

$ \boldsymbol{u}_{i P}=\boldsymbol{u}_i\left(\boldsymbol{X}_P, t\right)=\sum\limits_{I=1}^{N_n} \boldsymbol{u}_{i I} N_I\left(\boldsymbol{X}_P\right)=\sum\limits_{I=1}^{N_n} \boldsymbol{u}_{i I} N_{I P} $

(15)

$ \boldsymbol{v}_{i P}=\dot{\boldsymbol{u}}_{i P}=\sum\limits_{I=1}^{N_n} \dot{\boldsymbol{u}}_{i I} N_I\left(\boldsymbol{X}_P\right)=\sum\limits_{I=1}^{N_n} \dot{\boldsymbol{u}}_{i I} N_{I P} $

(16)

$ \boldsymbol{a}_{i P}^{\mathrm{s}}=\ddot{\boldsymbol{u}}_{i P}=\sum\limits_{I=1}^{N_n} \ddot{\boldsymbol{u}}_{i I} N_I\left(\boldsymbol{X}_P\right)=\sum\limits_{I=1}^{N_n} \ddot{\boldsymbol{u}}_{i I} N_{I P} $

(17)

$ \delta u_{i P}=\sum\limits_{I=1}^{N_n} \delta u_{i I} N_I\left(\boldsymbol{X}_P\right)=\sum\limits_{I=1}^{N_n} \delta u_{i I} N_{I P} $

(18)

将式(11)、(14)~(18)代入到式(12)、(13)中，最终得固液两相耦合物质点算法的关键公式：

$ \sum\limits_{I=1}^{N_n} M_{I J} \ddot{\boldsymbol{u}}_{i J}=\boldsymbol{F}_{i I, \mathrm{int}}+\boldsymbol{F}_{i I, \mathrm{ext}}=\boldsymbol{F}_{i I} $

(19)

$ \sum\limits_{J=1}^{N_n} V_{I J} \boldsymbol{R}_{i J}^{\mathrm{e}}=\boldsymbol{R}_{i I, \mathrm{int}}+\boldsymbol{R}_{i I, \mathrm{ext}} $

(20)

$ M_{I J}=\sum\limits_{P=1}^{N_P} M_P N_{I P} N_{J P} $

(21)

$ \boldsymbol{F}_{i I, \text { int }}=-\sum\limits_{P=1}^{N_P} M_P \frac{\boldsymbol{\sigma}_{i j P}}{\rho} N_{I P, j} $

(22)

$ \boldsymbol{F}_{i I, \mathrm{ext}}=\sum\limits_{P=1}^{N_P} M_P \boldsymbol{b}_{i P} N_{I P}+\int_{\partial \varOmega} \boldsymbol{\sigma}_{i j} n_j N_I \mathrm{d} A $

(23)

$ \boldsymbol{R}_{i I, \mathrm{int}}=\sum\limits_{P=1}^{N_P} V_{P}p N_{i P}-\sum\limits_{P=1}^{N_P} V_P \boldsymbol{\rho}^{\mathrm{wR}} \boldsymbol{a}_{i P}^{\mathrm{s}} N_{i P} $

(24)

$ \boldsymbol{R}_{i I, \mathrm{ext}}=\sum\limits_{P=1}^{N_P} V_P \rho^{\mathrm{wR}} \boldsymbol{b}_{i P}-\int_{\partial \varOmega} p n_i N_{I P} \mathrm{d} A $

(25)

式中：F_iI、F_{iI, int}与F_{iI, ext}分别代表背景网格节点合力、内力与外力矢量，R_{iI, int}与R_{iI, ext}分别表示背景网格节点渗流过程中的拖拽内力与外力矢量。

1.3 边界条件与物质点自由边界追踪

针对问题域边界条件，在程序中采用位移边界(∂Ω_u)、速度边界(∂Ω_v)、外力边界(∂Ω_f)、压强边界条件(∂Ω_p)和流量边界条件(∂Ω_q)，由Dirichlet和Neumann两种不同类型的边界共同组成:

$ \left\{\begin{array}{l} \boldsymbol{u}_i=\overline{\boldsymbol{u}}_i, \partial \varOmega_u \\ \boldsymbol{v}_i=\overline{\boldsymbol{v}}_i, \partial \varOmega_v \\ \boldsymbol{\sigma}_{i j} n_{j, \tau}^b=\boldsymbol{\tau}_i, \varOmega_f \\ p=\bar{p}, \partial \varOmega_p \\ \boldsymbol{q}_n=n_{i, q}^b \boldsymbol{R}_i^{\text{e}} \cdot k / \gamma_f, \partial \varOmega_q \end{array}\right. $

(26)

式中：u_i为指定固体骨架位移，v为边界指定速度，τ_i为指定边界荷载，n_{j, τ}^b为外力边界外法线矢量值，p为边界指定压强，q_n为指定流体流量，n_{i, q}^b为流量边界外法线矢量值。

当存在自由渗流边界时，除上述边界设定外，还需对自由渗流边界位置处的背景网格节点或物质点孔压值施加Dirichlet边界条件。在不考虑表面张力影响下，施加0压强条件，即p^t+Δt|∂Ω_p=0。

由于在物质点法数值计算中，经常会涉及大变形计算，因此自由渗流边界位置会经常发生改变，因此在程序每循环步施加0压强条件前，需先对自由边界进行算法识别。本文采用体积分数法(图 3)，预先定义体积分数阈值，并通过式(27)计算网格节点单元体积分数，当单元计算所得体积分数小于阈值，且包含有与空单元共用节点时，则共用节点即为自由边界网格点，进而对其施加0压强条件。

$ \gamma=\frac{\sum\limits_{P=1}^{N_P} V_P}{V_{\mathrm{e}}} $

(27)

式中：γ为体积分数阈值，V_P为物质点体积，V_e为背景网格单元体积。

1.4 物质点法程序

图 4为所编制的物质点程序MapleLeaf的简要流程，图 5为程序内各子模块的结构关系。程序主要由3个部分组成：第1部分(Part 1)为数据读取及前处理阶段，主要用于模型数据文件的读入、模型生成及存储关键控制变量，生成物质点及背景网格必要信息，并进行数值模型自重应力初始化；第2部分(Part 2)为数值程序的核心，主要基于核心计算模块Core_subrou.f90，完成应力更新及孔隙水压力计算；第3部分(Part 3)为后处理模块，按要求生成可视化输出文件，对数值计算结果进行图形处理分析。

2 固液两相耦合物质点算法验证 2.1 静水压力柱模型

图 6(a)为长×宽×高=0.06 m×0.06 m×1.0 m的三维静水压力柱模型，模型内部土体视为饱和状态，初始有效应力及孔隙水压力场均为0。模型两侧和底面为非渗透边界，顶面为自由渗透边界。固相边界条件侧边为铰接，底面为固定位移。该模型将仅在自重作用下，从零初始应力场随时间最终获得稳定孔隙水压场。根据三维静水压力柱模型建立对应物质点模型，立方体网格单元尺寸0.02 m，共划分450个背景网格单元，每个单元格内设置单个物质点。数值计算时长为2.5 s，时间步长为1.0×10^-5 s，模型采用弹性本构，固液两相物理力学及本构参数见表 1。

在自重作用下，模型内不同高度位置处最终孔隙水压力理论值可采用式p=(h－z)γ^wR计算，其中h和z分别为柱体模型高度和物质点纵坐标。图 6(b)为物质点孔隙水压0 s及2.5 s时的数值结果分布。图 7(a)为不同高度位置处(z=0.01~ 0.99 m)物质点孔隙水压随时间的数值曲线。图 7(b)为不同数值时间时(t=0.001~2 s)模型孔隙水压沿高度方向上的分布曲线。由图 7分析可知，由于显式算法自身数值解的稳定性不足，不同高度处物质点孔隙水压力值将随时间振荡变化，最终振荡趋势减弱并达到稳定，模型最终形成稳定的孔隙水压场，不同高度水压数值解与理论值一致。

2.2 太沙基一维固结模型

如图 8所示，在静力水压力柱物质点模型的基础上进一步建立与太沙基一维固结理论相适应的数值模型。模型尺寸、物质点网格单元、物质点总数、渗透边界及位移边界条件设置均与上节静水压力柱模型一致。模型内部水压初始值为10 kPa，且模型顶面施加10 kPa均布超载，不考虑重力作用，在固结条件下，模型内部水压值将随时间不断耗散直至0 kPa，而模型内部竖向有效应力则不断增加，并最终承担所有外荷载。模型数值计算时长2.0 s，时间步长为1.0×10^-5 s，采用弹性本构，固液两相物理力学及本构参数详见表 1。

依据太沙基一维固结理论，模型内部不同高度水压理论值可根据式(28)~(30)计算。表 2为不同渗透系数(10^-2 ~ 10^-6 m/s)条件下土体底部在不同数值计算时间的物质点程序数值解与解析解。由表 2数据分析可知，当渗透系数较小时，土体排水条件较差，水压值消散速度较慢，完整模拟土体整个消散过程需设定较长的数值计算时间段及程序计算时间。当渗透系数较大时，土体排水条件较好，水压值消散速度较快，且对于物质点数值程序来说，渗透系数的增大意味着孔压增量值在单时间步的增大，为确保数值计算稳定性，需设定更小的时间步，因此，程序计算时间也会较长。图 9为渗透系数10^-3 m/s条件下不同高度处孔隙水压与初始水压比随时间的数值解与理论解对比曲线。由图 9分析可知，由于采用显式固液两相耦合物质点算法，在数值计算初期，物质点数值解随时间振荡明显，与理论解存在明显偏差，但随后随时间持续增加，数值解与理论值始终保持一致。

$ \frac{\partial p}{\partial t}=C_{\mathrm{v}} \frac{\partial^2 p}{\partial z^2} p=\sum\limits_{m=0}^{+\infty}\left(\frac{2 p_0}{M} \mathrm{e}^{-M^2 T_{\mathrm{v}}} \sin \frac{M z}{H}\right) $

(28)

$ T_{\mathrm{v}}=\frac{C_{\mathrm{v}} t}{H^2} $

(29)

$ M=\frac{{\mathtt{π}}}{2}(2 m+1), m=0, 1, 2, \cdots, +\infty $

(30)

式中：C_v为固结系数，T_v为无量纲时间，H为模型最大排水距离，在本文中为模型高度。

3 液化本构模型

由于物质点法可以描述大变形问题，而在饱和地基动力液化过程中，经常伴随大变形问题，且现有针对地基液化的MPM数值案例较少，因此，本程序在固液两相耦合物质点算法的基础上进一步提出了动力条件下的饱和地基液化数值计算方法。首先在程序中引入了能够描述土体动力液化行为的Cyclic Mobility本构模型(CM模型)^[20-22]。如图 10所示，CM模型在修正剑桥模型的基础上，引入了上、下负荷面和应力诱导各向异性的概念，并结合土体超固结性和结构性的发展准则，能够较好地描述土体在循环荷载作用下的力学行为，特别是砂土的液化力学特性。

在CM模型中，分别定义用于描述土体结构特性及土体超固结性的状态参量R^*和R：

$ R^*=\frac{\tilde{p}^{\prime}}{\bar{p}^{\prime}}=\frac{\tilde{q}}{\bar{q}}, 0<R^* \leqslant 1 $

(31)

$ R=\frac{p^{\prime}}{\bar{p}^{\prime}}=\frac{q}{\bar{q}}, 0<R \leqslant 1 $

(32)

式中：应力坐标对(p′, q)、($\tilde{p}^{\prime}, \tilde{q} $)和($\bar{p}^{\prime}, \bar{q} $)分别代表下负荷面(即当前土体)应力状态、正常固结应力状态和具有上负荷面(即具有结构性的自然土体)应力状态。

土体当前应力状态始终位于下负荷面上，对应屈服面本构方程(33)，该模型包含8个参数，其中5个参数与剑桥模型参数相同。针对CM模型参数的具体介绍，详见文献[20-22]。

$ f=c \ln \frac{p^{\prime}}{\tilde{p}_0^{\prime}}+c \ln \frac{M^2-\zeta^{-}+\eta^{* 2}}{M^2-\zeta^2}+c \ln R^*-c \ln R-\varepsilon_V^p=0 $

(33)

式中：$c=\frac{\lambda-\kappa}{1+e_0}, \eta=\sqrt{\frac{3}{2} \boldsymbol{\eta} \cdot \boldsymbol{\eta}}, \zeta=\sqrt{\frac{3}{2} \boldsymbol{\beta} \cdot \boldsymbol{\beta}} $，$\boldsymbol{S}=\boldsymbol{\sigma}-p^{\prime} \boldsymbol{I}, \boldsymbol{\eta}=\frac{\boldsymbol{S}}{p^{\prime}}, \boldsymbol{\eta}^*=\sqrt{\frac{3}{2} \hat{\boldsymbol{\eta}} \cdot \hat{\boldsymbol{\eta}}}, \hat{\boldsymbol{\eta}}=\boldsymbol{\eta}-\boldsymbol{\beta}, \tilde{p}_0 $为参考应力98 kPa，p′为平均有效应力，M为临界状态线斜率，ζ为描述应力诱导异性大小的参量，β 为各向异性张量，η 为应力比张量，η^*为应力比张量与各向异性张量之差，S为偏应力张量，I为Kronecker符号，ε_V^p为塑性体积应变，λ和κ分别为土体的压缩指数和膨胀指数，e₀为$\tilde{p}_0 $=98 kPa时对应的孔隙比。

4 饱和地基液化数值计算 4.1 动力算法

在固液两相物质点耦合算法基础上，进一步考虑动力条件下的模型数值计算结果。在动力条件下，基于达朗贝尔原理对动力加速度值进行静态等效，在模型整体物质点颗粒上统一施加一组附加力，即惯性力。设定动力加速度为a_out，则式(8)可进一步扩展为$\nabla \cdot\left(\boldsymbol{\sigma}^{\prime}-p \boldsymbol{I}\right)+\rho \boldsymbol{b}-\rho \ddot{\boldsymbol{u}}-\rho \boldsymbol{a}_{\text {out }}=0 $。进一步对该公式进行积分弱形式及离散化转化，将动力算法模块在程序中添加后，针对现有文献及现场制备的饱和地基模型进行MPM动力液化数值计算。

4.2 基于现有文献的饱和地基液化数值计算

文献[20]采用拍波形式正弦加速度输入波形，对模型箱内饱和砂土进行了振动台动力液化分析，给出了不同深度处超孔隙水压(excess pore water pressure, EPWP)试验结果及采用CM模型的FEM数值计算结果。本文采纳该文献中的相关模型材料参数及结果数据，建立MPM数值模型，并基于CM模型进行动力数值计算，并与文献中相关结果进行对比分析。如图 11(a)所示，振动台上模型箱长×宽×高=1.2 m×1.0 m ×0.8 m，模型箱内饱和地基采用日本Toyoura液化砂，层高0.6 m，并埋置5个孔隙水压力传感器，间隔0.1 m。图 11(b)为对应三维MPM数值模型，背景立方体网格单元尺寸为0.1 m，总共建立12×10×6=720个物质点，固液两相物理力学参数及CM模型参数见表 3。图 11(c)为本次数值计算所采用的15 s正弦拍波输入加速度波形，最大振幅值0.3g，频率为4 Hz。MPM数值计算时长100.0 s，时间步长为1.0×10^-5 s。

图 12(a)为饱和地基不同深度处EPWP值的现场振动台试验结果与MPM数值结果对比结果。图 12(b)为饱和地基不同深度处EPWP的FEM与MPM数值结果对比结果。由图 12(a)、12(b)分析可知，在拍波形式动力作用下，饱和地基内部不同深度处的EPWP值随振动时间呈先上升后下降的趋势，且埋置深度越深，EPWP值越高。同一埋置深度处，基于MPM与FEM的EPWP数值结果与现场振动台试验结果存在一定的差别，其中基于现场振动台试验的EPWP值累积及消散速度较快，这可能与数值计算中所设渗透系数^[23-24]、边界条件^[25]有关，此外，现场试验过程中砂土内部产生的导水通道^{[20, 26]}也将进一步影响EPWP的消散速度。结合图 12(c)由MPM算法计算得到的不同埋置深度处的竖向有效应力结果曲线可知，在整体振动过程中，EPWP值的逐渐累积伴随着相应位置处竖向有效应力值的逐渐减小，直至达到0点附近，产生液化状态。当振动停止后，随着EPWP值的消散，竖向有效应力值逐步恢复至初始状态。从整体上看，基于MPM的固液两相耦合数值计算方法能够反应动力条件下饱和地基液化前后EPWP的增长和消散的基本规律。

4.3 基于现场制备的饱和地基液化数值计算

利用图 13中所示的小型振动台系统现场制备饱和砂土地基，并对饱和地基进行振动台动力液化分析。模型箱长×宽×高=0.4 m×0.2 m×0.3 m，模型箱内底部为2 cm厚的非液化黏土层，上部为14 cm厚的液化层，液化层由石英砂经水中落砂法配置而成，不均匀系数为1.0，曲率系数为2.0，最大干密度为1.77 g/cm³，最小干密度为1.37 g/cm³，相对密度为32.8%，孔隙比为0.78。在液化层内深度7 cm位置放置单个孔隙水压力传感器。图 14为对应建立的饱和砂层MPM数值模型，背景立方体网格单元尺寸0.1 m，共建立20×10×7=1 400个物质点，数值计算所需的模型材料物理力学参数及本构参数与表 3一致。

图 15(a)为振动过程中的加速度波形，振动时间为20 s，频率为4 Hz，最大振幅值为0.2g左右。图 15(b)为饱和地基埋置深度7 cm位置处EPWP的现场试验数据及MPM模型数值计算结果图。该位置处初始竖向有效应力约为0.63 kPa。

由图 15可知，饱和地基振动过程中，传感器位置EPWP随时间变化曲线与MPM模型数值计算得到的曲线值基本一致，且随超孔压的不断累积，EPWP值达到初始竖向有效应力线，传感器所在位置土体产生液化现象，随振动时间的继续持续增加，由于土体内孔隙水的消散，EPWP曲线呈下降趋势。

5 结论

1) 推导了固液两相耦合物质点算法的基本公式，采用Fortran 95程序语言编制了固液两相耦合物质点法三维显式程序MapleLeaf。分别建立了静水压力柱MPM数值模型与太沙基一维固结理论MPM数值模型，通过模型数值解与理论值的对比分析，验证了数值方法在分析固液两相相互作用问题的可靠性。

2) 在固液两相耦合物质点算法的基础上，引入弹塑性本构CM模型，并在程序中增添动力算法模块，将程序进一步用于分析包含大变形的动力条件下土体液化行为。分别依据现有文献中的饱和地基模型及现场制备的饱和地基模型的相关参数及结果数据，与对应建立的MPM模型数值计算结果进行对比分析，检验了程序在分析动力条件下饱和地基在液化分析方面的准确性。

3) 编制的基于固液两相耦合物质点算法程序能够较好地模拟固液两相介质在静力及动力条件下的相互作用，考虑到物质点法相比与传统有限元算法的独特优点，后续将针对水下边坡液化等大变形问题进一步开展研究工作。