光滑粒子流体动力学方法在饱和边坡地震滑移大变形中的适用性研究

黄帅; 刘传正; GODAKatsuichiro

doi:10.11779/CJGE20211274

光滑粒子流体动力学方法在饱和边坡地震滑移大变形中的适用性研究 English Version

1.
应急管理部国家自然灾害防治研究院, 北京 100085
2.
韦仕敦大学, 伦敦 N6A3K7, 加拿大

基金项目:

国家自然科学基金项目 51708516

博士后基金项目 2021M691391

地震数值预测联合实验室开放基金项目 2021LNEF04

详细信息

作者简介:
作者简介：黄帅(1987—)，男，博士，副研究员，主要从事工程防震减灾研究。E-mail: shuaihuang@ninhm.ac.cn

中图分类号: TU435
计量
- 文章访问数: 291
- HTML全文浏览量: 75
- PDF下载量: 120
出版历程
- 收稿日期: 2021-12-06
- 网络出版日期: 2023-02-23
- 刊出日期: 2023-01-31

Applicability of smooth particle hydrodynamics method to large sliding deformation of saturated slopes under earthquake action

1.
National Institute of Natural Hazards, MEMC, Beijing 100085, China
2.
Western University, London N6A 3K7, Canada

摘要

摘要: 地震作用下饱和边坡破坏会产生滑移大变形行为，一般难以采用有限元等网格方法进行分析，而光滑粒子流体动力学方法可以轻松处理大变形问题。基于此，通过在光滑粒子流体动力学方法（SPH方法）中加入瑞利衰减阻尼，基于有效应力本构模型改进了SPH方法的运动方程，建立了用于饱和边坡地震滑移大变形的数值计算方法；之后与饱和土样的空心圆柱扭剪试验获取的有效应力路径和剪应力-应变滞回曲线进行了对比，取得了较为一致的结果，验证了本改进方法的可行性。最后采用建立的数值计算方法，分析了不同SPH粒子密度下饱和边坡的地震滑移行为并与Newmark方法进行了对比。研究发现，改进的SPH方法确定的饱和边坡地震滑移临界面形状与Newmark方法对数螺旋滑移面形状吻合较好；因此改进的SPH方法能较好的再现饱和边坡地震作用下的滑移大变形过程，但该方法中边坡地震滑移大变形分析结果受设置的粒子密度影响。
- 地震作用 /
- 滑移大变形 /
- 光滑粒子流体动力学方法 /
- Newmark方法 /
- 临界滑移面
Abstract: The failure of a saturated slope under an intense earthquake ground motion will produce large sliding deformation, which is difficult to analyze using the finite element method. However, the smooth particle hydrodynamics method (SPH method) can deal with large deformation easily. In this study, to develop a new numerical method for calculating the seismic sliding large deformation of saturated slopes, the motion equation for the SPH method is improved by (1) introducing the Rayleigh damping into the smooth particle hydrodynamics and (2) adopting the effective stress constitutive model. To evaluate the effectiveness of the new method at soil element levels, the effective stress path and shear stress-strain hysteresis curve are obtained through the hollow cylindrical torsional shear tests on saturated soil samples, and are compared with the simulated results. The consistent results are obtained to verify the feasibility of the improved method. Finally, the seismic sliding behavior of saturated slopes under different SPH particle densities is analyzed by using the established numerical method and compared with that by the Newmark's sliding block method. It is found that the key parameter of the improved method is the particle density and that the shape of seismic sliding critical surface of saturated slopes determined by the improved SPH method is in good agreement with that of logarithmic spiral slip surface determined by the Newmark's method. Therefore, the improved SPH method can also be used to simulate a large sliding deformation process of saturated slopes under seismic actions.
- earthquake action /
- large sliding deformation /
- smooth particle hydrodynamics method /
- Newmark method /
- critical slip surface

HTML全文

0. 引言

膨胀性非饱和土是一种具有特殊结构的非饱和土，通常表现出吸水膨胀、失水收缩和浸水强度降低等特性，工程中容易使建筑物产生不均匀沉降，甚至会导致结构开裂和破坏。因此，对膨胀性非饱和土本构模型的研究具有重要工程意义，能够为预测膨胀土的变形特性从而避免工程灾害的发生提供有效手段。

1990年，Alonso等^[1]在修正剑桥模型的基础上，考虑了吸力的影响并引入了加载-湿陷屈服面（LC屈服面）和吸力增加屈服面（SI屈服面），进而提出了巴塞罗那弹塑性本构模型（BBM），能够描述非饱和土屈服应力随吸力增加而变大、湿化变形等力学特性。在此之后，Wheeler等^[2]、Sun等^[3-4]、缪林昌^[5]、赵成刚等^[6]、姚仰平等^[7-8]、Luo等^[9]也提出了多种形式的非饱和土本构模型，这些模型虽然可以描述包括湿陷变形在内的非饱和土的基本力学特性，但不能体现膨胀性非饱和土吸水膨胀的力学特性。

为更好地描述膨胀性非饱和土的力学特性，Gens等^[10]于1992年建立了膨胀性非饱和土弹塑性本构模型的理论框架模型（G-A模型）。在此基础上，Alonso等^[11]进一步建立了巴塞罗那膨胀性非饱和土弹塑性本构模型（BExM）。该模型认为膨胀土的膨胀性是由于其具有双重结构，即微观结构与宏观结构，当非饱和膨胀土湿化时，微观结构会发生膨胀，从而导致宏观结构发生塑性膨胀，最终表现出土体的整体膨胀。在此之后，许多学者在这一模型的理论基础上进行了改进和扩展，Sánchez等^[12]在G-A模型的基础上建立了考虑微观-宏观双重结构影响的广义塑性模型，并将其成功运用于有限元数值分析中；Mašín^[13]分别建立了能够描述微观结构和宏观结构的计算公式，通过引入双重结构的耦合方程建立了考虑水力耦合的非饱和膨胀土模型；卢再华等^[14]在改进的G-A模型中引入剪切屈服面（S_y），建立了非饱和原状膨胀土的弹塑性损伤本构模型；陈正汉等^[15-16]进一步针对原状膨胀性土、重塑膨胀性土和重塑膨润土等不同类型的膨胀土进行了全方位的研究，综合考虑了温度、持水特性等因素对膨胀土力学性质的影响，建立了考虑热-水-力等多场耦合特性的膨胀土模型；孙德安等^[17]考虑了饱和度对应力应变关系的影响及变形对土水特性的影响，并引入了等孔隙比线，建立了能够描述膨胀性非饱和土水力-力学耦合特性的弹塑性本构模型；曹雪山^[18]针对膨胀土干湿循环特性进一步完善模型，引入了微-宏观结构耦合参数，提出了考虑双重结构耦合的弹塑性本构模型；李舰等^[19-20]，Li等^[21]建立了适用于膨胀性非饱和土的边界面模型，并围绕其提出了有效的数值计算方法。

Yao等^[22]以修正剑桥模型为基础，提出了可以描述饱和超固结土力学特性的统一硬化（UH）本构模型，之后进一步考虑了吸力的影响，引入LC屈服面，扩展了UH模型使其可以描述非饱和土的特性^[8]。本文在非饱和土UH模型基础上，引入考虑超固结度和初始干密度影响的膨胀式，使其可以反映膨胀土湿化膨胀过程中产生的塑性变形，并进一步构建了考虑膨胀土膨胀性的非饱和土UH模型，由于本模型仅考虑了膨胀土的非饱和、胀缩性和超固结性等特性，没有考虑裂隙性^[23]，是对理想状态的膨胀土进行描述，故称理想膨胀性非饱和土UH模型。

1. 膨胀性非饱和土UH模型

1.1 吸力应力与屈服面方程

Alonso等^[1]在BBM中认为非饱和土的应变由吸力和净应力共同控制，并提出了LC屈服面和SI屈服面的概念来描述非饱和土的变形特性。Luo等 ^[9]在此基础上，根据对试验现象的分析，提出了新的吸力应力表达式，

$\begin{array}{*{20}{c}} {{p_{\text{s}}} = a\left( {1 - {e^{ - \frac{s}{a}}}} \right)} \end{array}\ \ \ \ \ \ \ 。$

(1)

式中： ${p_{\text{s}}}$ 为吸力应力，表示吸力对土体黏聚力的贡献值^[16]，可以通过恒定吸力下非饱和土的三轴试验或直剪试验测得； $a$ 为最大吸力应力，可以通过两个或两个以上不同吸力下非饱和土的三轴试验或直剪试验求得^[3]； $s$ 为吸力，是负孔隙水压力和表面张力在土骨架上的综合作用应力。区别于BBM，这里吸力应力与吸力并非正比关系，而是呈现非线性变化并存在一个最大值。

对于正常固结非饱和土，根据Sun等^[3]、Yao等^[8]和Luo等^[9]的研究，假设其等向压缩线在 $e{\text{ - ln}}(p + {p_{\text{s}}})$ 平面呈直线，如图 1所示，其等向压缩线公式可表示为

图 1 饱和土和非饱和土的等向压缩线

Figure 1. Normal compression lines for saturated and unsaturated soils

下载: 全尺寸图片幻灯片

$\begin{array}{*{20}{c}} {e = e(s) - \lambda \ln \frac{{p + {p_{\text{s}}}}}{{{p^{\text{c}}} + {p_{\text{s}}}}}} \end{array}\ \ \ \ \ \ \ 。$

(2)

式中： ${p^{\text{c}}}$ 为参考吸力应力，假设当净平均应力为 ${p^{\text{c}}}$ 时，吸力变化仅引起弹性应变， $\lambda$ 为正常固结饱和土在 $e{\text{ - ln}}p$ 平面上等向压缩线的斜率。当吸力恒定时，平均净应力变化所引起的变形可写为

$\begin{array}{*{20}{c}} {{\text{d}}e_{\text{p}}^{\text{e}} = - \kappa \frac{{{\text{d}}p}}{p}} \end{array}\ \ \ \ \ \ \ 。$

(3)

对于吸力变化引起的变形，假设历史最大吸力为 ${s_0}$ ，当吸力 $s$ 小于历史最大吸力时，吸力变化仅引起弹性变形，当吸力 $s$ 大于历史最大吸力 ${s_0}$ 时，吸力增加会引发弹塑性变形，定量表述如下：

$\begin{array}{c}de=\{\begin{array}{c}-{\kappa }_{\text{s}}\frac{\text{d}s}{s+{p}_{\text{at}}}\text{(}s\le {s}_{0}\text{)}\\ -{\lambda }_{\text{s}}\frac{\text{d}s}{s+{p}_{\text{at}}}\text{(}s＞{s}_{0}\text{)}\end{array}\end{array} \ \ \ \ \ \ \ 。$

(4)

式中： ${\lambda _{\text{s}}}$ 为吸力超过历史最大吸力 ${s_0}$ 之后非饱和土在 $e{\text{ - ln}}s$ 平面内吸力加载线的斜率， ${\kappa _{\text{s}}}$ 为非饱和土在 $e{\text{ - ln}}s$ 平面内吸力回弹线的斜率， ${p_{{\text{at}}}}$ 为标准大气压。

基于上述假设，新的LC屈服面方程可以推导为

$\begin{array}{*{20}{c}} {{{\left( {\frac{{p_{\text{x}}^*}}{{{p^{\text{c}}}}}} \right)}^{\lambda - \kappa }} = {{\left( {\frac{{{p_{\text{x}}} + {p_{\text{s}}}}}{{{p^{\text{c}}} + {p_{\text{s}}}}}} \right)}^\lambda }{{\left( {\frac{{{p_{\text{x}}}}}{{{p^{\text{c}}}}}} \right)}^{ - \kappa }}} \end{array}\ \ \ \ \ \ \ 。$

(5)

式中： ${p_{\text{x}}}$ 为不同吸力下 $p{\text{ - }}q$ 平面屈服面与 $p$ 轴的交点值， $p_{\text{x}}^{\text{*}}$ 为饱和土 $p{\text{ - }}q$ 平面屈服面与 $p$ 轴的交点值，如图 2。

图 2 LC屈服面和SI屈服面

Figure 2. Yield curves of LC and SI

下载: 全尺寸图片幻灯片

SI屈服面采用BBM中与 $p$ 轴平行的直线，如图 2，方程可写为

$\begin{array}{*{20}{c}} {s = {s_0}} \end{array}\ \ \ \ \ \ \ 。$

(6)

式中： ${s_0}$ 为历史最大吸力。

$p{\text{ - }}q$ 平面采用椭圆屈服面，

$f = {\text{ln}}\hat p + \ln \left( {1 + \frac{{{{\hat q}^2}}}{{{M^2}{{\hat p}^2}}}} \right) - \ln ({p_{\text{x}}} + {p_{\text{s}}}) = 0\ \ \ \ \ \ \ 。$

(7)

式中： $\hat p = p + {p_{\text{s}}}$ ， $\hat q = q$ 。本文提出的模型采用相关联流动法则，假设塑性势面方程与屈服面方程相同。

1.2 硬化参量

本模型屈服面的硬化参数采用 $H$ ，硬化规律可表示为

$\begin{array}{*{20}{c}} {\frac{{{\text{d}}p_{\text{x}}^*}}{{p_{\text{x}}^*}} = \frac{{1 + {e_0}}}{{\lambda - \kappa }}dH} \end{array}\ \ \ \ \ \ \ 。$

(8)

Alonso等^[11]在提出的BExM中认为膨胀土吸水膨胀与其微观结构有关，在吸力变化时土体微观结构的改变会引起宏观结构的塑性变形。对于膨胀土，湿化会使得微观结构膨胀，进而引发宏观结构的塑性膨胀。膨胀土在湿化时塑性体应变由两部分组成，一部分是净应力变化和吸力变化综合引发的塑性体应变 ${\text{d}}\varepsilon _{{\text{vp}}}^{\text{p}}$ ，与非膨胀性非饱和土中的塑性体应变相同，另一部分为微观结构变化引起的宏观结构塑性体应变 ${\text{d}}\varepsilon _{{\text{vs}}}^{\text{p}}$ 。卢再华等^[24]提出对膨胀性非饱和土湿化时的膨胀量，可采用其与土体湿化产生的弹性体变的关系来求解，参照Alonso等^[11]和卢再华等^[24]对膨胀变形的求解思路，本文假定微观结构的应变与吸力引起的弹性应变相同，将宏观结构发生的塑性体应变 ${\text{d}}\varepsilon _{{\text{vs}}}^{\text{p}}$ 表示为微观结构体应变的函数：

$\begin{array}{*{20}{c}} {{\text{d}}\varepsilon _{{\text{vs}}}^{\text{p}} = {f_{\text{d}}}{\text{d}}\varepsilon _{{\text{vs}}}^{\text{e}}} \end{array}\ \ \ \ \ \ \ 。$

(9)

式中：

$\begin{array}{*{20}{c}} {{\text{d}}\varepsilon _{{\text{vs}}}^{\text{e}} = \frac{{{\kappa _{\text{s}}} \cdot {\text{d}}s}}{{\left[ {(1 + {e_0})(s + {p_{{\text{at}}}})} \right]}}} \end{array}\ \ \ \ \ \ \ \text{；}$

(10)

${f_{\text{d}}}$ 为放大系数，

$\begin{array}{*{20}{c}} {{f_{\text{d}}} = {J^{\frac{{{\rho _{\text{d}}}}}{{{\rho _{\text{w}}}}}}}{{(1 - R)}^2}} \end{array}\ \ \ \ \ \ \ 。$

(11)

式中： $J$ 为反映膨胀土膨胀能力的参数； ${\rho _{\text{d}}}$ 为土体初始干密度，根据高登辉等^[25]的研究，初始干密度对非饱和土的变形参数、强度参数和渗透参数等有重要影响，此处通过引入初始干密度来考虑其对非饱和土膨胀量的影响， ${\rho _{\text{w}}}$ 为水的密度， $R$ 为超固结系数，可通过孔隙状态参数 $\xi$ 求出，

$\begin{array}{*{20}{c}} {R = \exp \left( { - \frac{\xi }{{\lambda - \kappa }}} \right)} \end{array}\ \ \ \ \ \ \ 。$

(12)

为考虑膨胀性的影响，在屈服面中引入微观结构引起的塑性变形，即将塑性应变增量表示为两部分：

$\begin{array}{*{20}{c}} {{\text{d}}\varepsilon _{\text{v}}^{\text{p}} = {\text{d}}\varepsilon _{{\text{vp}}}^{\text{p}} + {\text{d}}\varepsilon _{{\text{vs}}}^{\text{p}}} \end{array}\ \ \ \ \ \ \ 。$

(13)

屈服面采用统一硬化参量 $H$ ，并认为统一硬化参量 $H$ 只与净应力引起的塑性体应变有关，具体增量表达式为

$\begin{array}{*{20}{c}} {{\text{d}}H = \frac{{M_{\text{f}}^4 - {{\hat \eta }^4}}}{{{M^4} - {{\hat \eta }^4}}}{\text{d}}\varepsilon _{{\text{vp}}}^{\text{p}}} \end{array}\ \ \ \ \ \ \ 。$

(14)

式中： $\hat \eta$ 为应力比， $\hat \eta = \hat q/\hat p$ ； ${M_{\text{f}}}$ 为潜在强度^[26]，

$\begin{array}{*{20}{c}} {{M_{\text{f}}} = 6\left[ {\sqrt {\frac{{12(3 - M)}}{{{M^2}}}\exp \left( { - \frac{\xi }{{\lambda - \kappa }}} \right) + 1} + 1} \right]} \end{array}\ \ \ \ \ \ \ 。$

(15)

式中： $\xi$ 为孔隙状态参数，反映土体当前孔隙比下的超固结程度，

$\begin{array}{*{20}{c}} {\xi = {e_{\hat \eta }} - e} \end{array}\ \ \ \ \ \ \ 。$

(16)

式中： ${e_{\hat \eta }}$ 为当前应力状态下土体正常固结时对应的孔隙比，

${e_{\hat \eta }} = e(s) - \lambda \ln \left[ {\frac{{\hat p}}{{{p^{\text{c}}} + {p_{\text{s}}}}}\left( {1 + \frac{{{{\hat \eta }^2}}}{{{M^2}}}} \right)} \right] + \kappa \ln \left( {1 + \frac{{\hat p}}{p}\frac{{{{\hat \eta }^2}}}{{{M^2}}}} \right) \ \ \ \ \ \ \ 。$

(17)

1.3 弹性应变增量

在净应力和吸力发生变化时，土体会发生变形，其中弹性应变可以表示为弹性体应变和弹性剪应变，本模型中弹性体应变可以表示为

$\begin{array}{*{20}{c}} {{\text{d}}\varepsilon _{\text{v}}^{\text{e}} = {\text{d}}\varepsilon _{{\text{vp}}}^{\text{e}} + {\text{d}}\varepsilon _{{\text{vs}}}^{\text{e}}} \end{array}\ \ \ \ \ \ \ 。$

(18)

其包含两部分，第一部分为净应力变化引起的弹性体应变，

$\begin{array}{*{20}{c}} {{\text{d}}\varepsilon _{{\text{vp}}}^{\text{e}} = \frac{k}{{1 + {e_0}}}\frac{{{\text{d}}\hat p}}{p} = \frac{{{\text{d}}\hat p}}{K}} \end{array}\ \ \ \ \ \ \ 。$

(19)

式中：

$\begin{array}{*{20}{c}} {K = \frac{E}{{3{\text{(}}1 - 2v{\text{)}}}}} \end{array} \ \ \ \ \ \ \ \text{；}$

(20)

$\begin{array}{*{20}{c}} {E = \frac{{3(1 - 2v)(1 + {e_0})}}{\kappa }\hat p} \end{array}\ \ \ \ \ \ \ 。$

(21)

第二部分为吸力变化引起的弹性体应变，方程表达式见式（10）。

本模型认为吸力变化不会引起弹性剪应变，故弹性剪应变仅由净应力中的剪应力引起，可表示为

$\begin{array}{*{20}{c}} {{\text{d}}\varepsilon _{\text{d}}^{\text{e}} = \frac{{{\text{d}}\hat q}}{{3G}}} \end{array}\ \ \ \ \ \ \ 。$

(22)

式中：

$\begin{array}{*{20}{c}} {G = \frac{E}{{2(1 + v)}}} \end{array}\ \ \ \ \ \ \ 。$

(23)

1.4 塑性应变增量

由屈服面公式可知，净应力变化引起的应变增量可表示为

$\begin{array}{*{20}{c}} {{\text{d}}\varepsilon _{{\text{vp}}}^{\text{p}} = \mathit{\Lambda } \frac{{\partial f}}{{\partial \hat p}}} \end{array}\ \ \ \ \ \ \ \text{，}$

(24)

$\begin{array}{*{20}{c}} {{d} \varepsilon _{\text{d}}^{\text{p}} = \mathit{\Lambda } \frac{{\partial f}}{{\partial \hat q}}} \end{array}\ \ \ \ \ \ \ 。$

(25)

式中： $\mathit{\Lambda }$ 为塑性标量因子。

对当前屈服面函数和LC屈服面函数进行全微分，联立求解后得到

$\begin{array}{*{20}{c}} {\mathit{\Lambda } = {A_0}\frac{{\frac{{\partial f}}{{\partial \hat p}}{\text{d}}\hat p + \frac{{\partial f}}{{\partial \hat q}}{\text{d}}\hat q + {f_{\text{s}}}{\text{d}}s}}{{\frac{{\partial f}}{{\partial \hat p}}}}} \end{array}\ \ \ \ \ \ \ 。$

(26)

式中：

$\begin{array}{*{20}{c}} {{A_0} = \frac{{{c_{\text{p}}}\hat {\mathit{\Omega }} }}{{p_{\text{x}}^{\text{*}}}}\frac{{\lambda {p_{\text{x}}} - \kappa \left( {{p_{\text{x}}} + {p_{\text{s}}}} \right)}}{{\lambda - \kappa }}{{\left( {\frac{{{p_{\text{x}}}}}{{{p^{\text{c}}}}}} \right)}^{\frac{{ - \lambda }}{{\lambda - \kappa }}}}{{\left( {\frac{{{p_{\text{x}}} + {p_{\text{s}}}}}{{{p^{\text{c}}} + {p_{\text{s}}}}}} \right)}^{\frac{\lambda }{{\lambda - \kappa }}}}} \end{array}\ \ \ \ \ \ \ \text{；}$

(27)

$\begin{array}{*{20}{c}} {{f_{\text{s}}} = \frac{{ - \lambda {p_{\text{x}}} + \kappa ({p^{\text{c}}} + {p_{\text{s}}})}}{{({p^{\text{c}}} + {p_{\text{s}}})\left[ {\lambda {p_{\text{x}}} - \kappa ({p_{\text{x}}} + {p_{\text{s}}})} \right]}}{{\text{e}}^{ - \frac{s}{a}}}} \end{array}\ \ \ \ \ \ \ \text{；}$

(28)

$\begin{array}{*{20}{c}} {p_{\text{x}}^{\text{*}} = {p^{\text{c}}}{{\left( {\frac{{{p_{\text{x}}} + {p_{\text{s}}}}}{{{p^{\text{c}}} + {p_{\text{s}}}}}} \right)}^{\frac{\lambda }{{\lambda - \kappa }}}}{{\left( {\frac{{{p_{\text{x}}}}}{{{p^{\text{c}}}}}} \right)}^{\frac{{ - \kappa }}{{\lambda - \kappa }}}}} \end{array}\ \ \ \ \ \ \ \text{；}$

(29)

$\begin{array}{*{20}{c}} {{p_{\text{x}}} = p + \frac{{{{\hat q}^2}}}{{{M^2}{{\hat p}^2}}}} \end{array}\ \ \ \ \ \ \ \text{；}$

(30)

$\begin{array}{*{20}{c}} {\hat{\mathit{\Omega }} = \frac{{{M^4} - {{\hat \eta }^4}}}{{M_{\text{f}}^4 - {{\hat \eta }^4}}}} \end{array}\ \ \ \ \ \ \ 。$

(31)

综上可得：

$\begin{array}{*{20}{c}} {{\text{d}}\varepsilon _{{\text{vp}}}^{\text{p}} = \mathit{\Lambda } \frac{{\partial f}}{{\partial \hat p}} = {A_0}\left( {\frac{{\partial f}}{{\partial \hat p}}{\text{d}}\hat p + \frac{{\partial f}}{{\partial \hat q}}{\text{d}}\hat q + {f_{\text{s}}}{\text{d}}s} \right)} \end{array}\ \ \ \ \ \ \ \text{，}$

(32)

$\begin{array}{*{20}{c}} {{\text{d}}\varepsilon _{\text{d}}^{\text{p}} = \mathit{\Lambda } \frac{{\partial f}}{{\partial \hat q}} = {A_0}\frac{{\frac{{\partial f}}{{\partial \hat p}}{\text{d}}\hat p + \frac{{\partial f}}{{\partial \hat q}}{\text{d}}\hat q + {f_{\text{s}}}{\text{d}}s}}{{\frac{{\partial f}}{{\partial \hat p}}}}\frac{{\partial f}}{{\partial \hat q}}} \end{array}\ \ \ \ \ \ \ 。$

(33)

对于膨胀性非饱和土，除了上述体应变外，还有微观结构变化引起的体应变，当吸力减小时（即湿化），相应塑性体应变表达式见式（9）。

当吸力增加时，考虑LC屈服面和SI屈服面的耦合，将吸力增加引起的塑性体应变并入 ${\text{d}}\varepsilon _{{\text{vp}}}^{\text{p}}$ 考虑，此时 ${f_{\text{d}}} = 0$ ，即得到

$\begin{array}{*{20}{c}} {{f_{\text{d}}} = \left\{ \begin{gathered} {J^{\frac{{{\rho _{\text{d}}}}}{{{\rho _{\text{w}}}}}}}{\left( {1 - {\text{exp}}\left( { - \frac{\xi }{{\lambda - \kappa }}} \right)} \right)^2}{\text{ (d}}s < 0) \hfill \\ 0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {\text{(d}}s \geqslant 0)\ \ \hfill \\ \end{gathered} \right.} \end{array} \ \ \ \ \ \ \ 。$

(34)

1.5 弹塑性应力应变关系

膨胀土塑性体应变由 ${\text{d}}\varepsilon _{{\text{vp}}}^{\text{p}}$ 和 ${\text{d}}\varepsilon _{{\text{vs}}}^{\text{p}}$ 两部分组成，弹性体应变由 ${\text{d}}\varepsilon _{{\text{vp}}}^{\text{e}}$ 和 ${\text{d}}\varepsilon _{{\text{vs}}}^{\text{e}}$ 两部分组成，将应力应变关系表示为矩阵得到

$\begin{array}{*{20}{c}} {\left\{ {\begin{array}{*{20}{c}} {{\text{d}}\hat p + {B_1}{\text{d}}s} \\ {{\text{d}}\hat q + {B_2}{\text{d}}s} \end{array}} \right\} = \left( {\begin{array}{*{20}{c}} {K{A_1}}&{3KG{A_2}} \\ {3KG{A_2}}&{3G{A_3}} \end{array}} \right)\left\{ {\begin{array}{*{20}{c}} {{\text{d}}{\varepsilon _{\text{v}}}} \\ {{\text{d}}{\varepsilon _{\text{d}}}} \end{array}} \right\}} \end{array}\ \ \ \ \ \ \ 。$

(35)

式中：

$\begin{array}{*{20}{c}} {{A_1} = \frac{{12G{A_0}{{\hat \eta }^2} + \hat p({M^4} - {{\hat \eta }^4})}}{{{T_0}}}} \end{array}\ \ \ \ \ \ \ \text{；}$

(36)

$\begin{array}{*{20}{c}} {{A_2} = \frac{{ - 2{A_0}\hat \eta ({M^2} - {{\hat \eta }^2})}}{{{T_0}}}} \end{array}\ \ \ \ \ \ \ \text{；}$

(37)

$\begin{array}{*{20}{c}} {{A_3} = \frac{{{A_0}K{{({M^2} - {{\hat \eta }^2})}^2} + \hat p({M^4} - {{\hat \eta }^4})}}{{{T_0}}}} \end{array}\ \ \ \ \ \ \ \text{；}$

(38)

$\begin{array}{*{20}{c}} {{B_1} = K\frac{{12G{A_0}{{\hat \eta }^2}({T_1} - {A_0}{f_{\text{s}}}) + \hat p{T_1}({M^4} - {{\hat \eta }^4})}}{{{T_0}}}} \end{array}\ \ \ \ \ \ \ \text{；}$

(39)

${B_2} = 6G{A_0}\hat \eta \frac{{K({M^2} - {{\hat \eta }^2})({A_0}{f_{\text{s}}} - {T_1}) + \hat p{f_{\text{s}}}({M^2} + {{\hat \eta }^2})}}{{{T_0}}}\ \ \ \ \ \ \ \text{；}$

(40)

$\begin{array}{*{20}{c}} {{T_0} = {A_0}K{{({M^2} - {{\hat \eta }^2})}^2} + 12G{A_0}{{\hat \eta }^2} + \hat p({M^4} - {{\hat \eta }^4})} \end{array}\ \ \ \ \ \ \ \text{；}$

(41)

$\begin{array}{*{20}{c}} {{T_1} = {A_0}{f_{\text{s}}} + (1 + {f_{\text{d}}})\frac{{{k_{\text{s}}}}}{{(1 + {e_0})(s + {p_{{\text{at}}}})}} - \frac{{{{\text{e}}^{ - \frac{s}{a}}}}}{K}} \end{array}\ \ \ \ \ \ \ 。$

(42)

2. 模型特性分析

本节通过对几种湿化应力路径的模拟和分析，来体现本文模型在描述膨胀性非饱和土湿化过程中体积膨胀和强度降低的特性。所取模型参数： $\lambda$ =0.2， $\kappa$ =0.02， ${\lambda _{\text{s}}}$ =0.08， $\nu$ =0.3，M=1.0，N=1.96， ${\kappa _{\text{s}}}$ =0.01， ${p^{\text{c}}}$ =100 kPa，a=150 kPa， ${p_{{\text{at}}}}$ =100 kPa，初始孔隙比取0.95。

本模型中膨胀性非饱和土膨胀性的体现主要与3个参数有关，即反映土体膨胀能力的系数J、反映土体密实程度的初始干密度 ${\rho _{\text{d}}}$ 以及反映土体超固结程度的系数 $\xi$ ，其中系数J和初始干密度 ${\rho _{\text{d}}}$ 的影响可以通过膨胀系数 $C = {J^{\frac{{{\rho _{\text{d}}}}}{{{\rho _{\text{w}}}}}}}$ 综合反映，系数 $\xi$ 的影响可以通过前期固结压力体现，相同应力路径和吸力状态下，前期固结压力越大， $\xi$ 越大，即土的超固结程度越高。

2.1 等向回弹再湿化

为了分析膨胀系数对膨胀量的影响，先在吸力不变的条件下将土样等向压缩后回弹，然后湿化至饱和，应力路径如中A—B( $B'$ )—C—D所示，其中A—B为等向压缩， $B'$ —C为回弹，C—D为湿化。

图 3 不同C值等向回弹湿化应力路径

Figure 3. Stress paths of isotropic unloading-wetting with different values of C

下载: 全尺寸图片幻灯片

分别取膨胀系数C为0，2，5，在吸力s=600 kPa不变的条件下，将土样从初始状态（围压p=200 kPa）等向压缩至围压p=1000 kPa，然后回弹至p=200 kPa，然后再湿化至饱和（s=0 kPa），由本文所提模型预测湿化过程中土样的体应变，如图 4所示。

图 4 等向回弹湿化中不同C值下体应变

Figure 4. Variation of volumetric strain under stress path of isotropic unloading-wetting with different values of C

下载: 全尺寸图片幻灯片

可以看出，随着土体中吸力降低，土体会发生膨胀变形（负体应变代表膨胀变形），随着膨胀系数C的提高，膨胀量会随之增大；同时在吸力较低的时候土体会发生湿陷，这是因为在吸力接近0 kPa时，土样应力状态接近屈服面，吸力降低使得土体产生的塑性湿陷体应变增量 ${\text{d}}\varepsilon _{{\text{vp}}}^{\text{p}}$ 在数值上超过微观结构引起的塑性膨胀体应变增量 ${\text{d}}\varepsilon _{{\text{vs}}}^{\text{p}}$ ，塑性湿陷变形在总变形中占据主导地位，土体整体表现出体积减小。

为分析等向回弹再湿化路径下不同超固结度对膨胀量的影响，取初始围压p=200 kPa，膨胀系数C=2，吸力s保持为600 kPa，将土样分别等向压缩至500，1000，1500 kPa，再回弹至200 kPa，应力路径如图 5所示A—B₁（B₂，B₃）为等向压缩过程，B₁（B₂，B₃）—C为回弹过程，C—D为湿化过程。

图 5 不同超固结度等向回弹湿化路径

Figure 5. Stress paths of isotropic unloading-wetting with different overconsolidation ratios

下载: 全尺寸图片幻灯片

3个土样湿化时的超固结程度不同，由本文所提模型预测体应变在湿化过程中的变化，如图 6所示。

图 6 等向回弹湿化不同超固结度下体应变

Figure 6. Variation of volumetric strain under stress path of isotropic unloading-wetting with different overconsolidation ratios

下载: 全尺寸图片幻灯片

可以看出，超固结程度越高，膨胀量越大，当超固结度过低时（例如前期固结压力500 kPa），甚至会发生土体整体湿陷的现象，即湿化至饱和后土体体积相比湿化前减小，体现超固结程度对膨胀性非饱和土膨胀量的重要影响。

2.2 回弹—等p剪切—湿化—加载

此路径中，吸力保持600 kPa不变，等向加载回弹后，进行等p剪切，剪切到一定程度时湿化，之后再剪切至破坏，分析剪切过程中体应变的变化。

首先预测不同膨胀系数C在三轴排水剪切过程中的影响，初始围压p取200 kPa，膨胀系数C取0，2，5，等向压缩至p=1000 kPa，然后再回弹至200 kPa，之后进行等p剪切，剪切到广义剪应力q=200 kPa时进行湿化，湿化至饱和后再继续剪切，直到破坏，如图 7所示为剪切过程中体应变的变化。

图 7 等p剪切中湿化不同C值下体应变

Figure 7. Variation of volumetric strain under stress path of shearing- wetting-shearing with different values of C and constant mean effective stress

下载: 全尺寸图片幻灯片

可以看出膨胀系数C的取值不影响湿化前土体的应力应变关系，但在湿化过程中膨胀系数大的土体表现出更大的膨胀变形，甚至会使土样体积大于初始体积，这一变形导致湿化后剪切过程中土样剪缩量降低，但在不断剪切过程中，不同参数的土样体积应变趋于一致。同时通过分析湿化过程剪应力的变化，发现随着膨胀系数C值的增大，湿化后的剪切强度会降低。

其次分析此路径下不同超固结度对膨胀特性的影响。设置初始围压为200 kPa，膨胀系数C取2，分别将土样等向压缩至500，1000，1500 kPa，然后等向回弹至200 kPa，之后等p剪切，当q=200 kPa时湿化至饱和，再剪切至破坏，分析此剪切过程中体应变与轴向应变的关系，如图 8所示。

图 8 等p剪切中湿化不同超固结度下体应变

Figure 8. Variation of volumetric strain under stress path of shearing- wetting-shearing with different overconsolidation ratios and constant mean effective stress

下载: 全尺寸图片幻灯片

可以看出不同超固结度下湿化时产生的体应变不同，当超固结度较大时，土体剪缩量会减小，使得土样提前进入剪胀阶段，同时增大剪切过程中产生的总膨胀量。

2.3 等体积湿化中的膨胀力

膨胀力是评价膨胀性非饱和土膨胀特性的基本参数，试验中可以通过三向胀缩仪测定^[27]。在土体各向同性条件下，模型中的膨胀力与 ${\text{d}}{\varepsilon _{\text{v}}} = 0$ 时的平均正应力相等。本节研究土体膨胀过程中膨胀力的变化特性。

首先预测不同膨胀系数下土体膨胀力变化，初始吸力600 kPa，初始围压200 kPa，膨胀系数C分别取0，2，5，将土样等向压缩至1000 kPa，后回弹至200 kPa，改变围压维持体积不变，将土样湿化至饱和，预测湿化过程中膨胀力与吸力的关系，如图 9。

图 9 等体积湿化不同C值膨胀力变化

Figure 9. Variation of swelling pressure during wetting under constant volume with different values of C

下载: 全尺寸图片幻灯片

可以看出随着膨胀系数的提高，土体最大膨胀力也相应提高，同时在湿化至低吸力阶段时，土体膨胀力会降低，反映此阶段土体整体表现出湿陷趋势，相应的，膨胀系数C越大，低吸力阶段土体膨胀力降低越多，体现土体膨胀系数的增大会加剧产生土体塑性湿陷体应变增量 ${\text{d}}\varepsilon _{{\text{vp}}}^{\text{p}}$ 的趋势，使得土体整体表现出的体缩趋势更大，在本模拟中由于限制了土体变形，体缩趋势表现为膨胀力的减小。

其次预测不同超固结程度下膨胀力的变化，初始吸力仍为600 kPa，膨胀系数C取2，初始围压为200 kPa，将土样分别等向压缩至500，1000，1500 kPa，然后回弹至200 kPa，保持体积不变湿化至饱和，其膨胀力与吸力的关系如图 10。

图 10 等体积湿化不同超固结度膨胀力变化

Figure 10. Variation of swelling pressure during wetting under constant volume with different overconsolidation ratios

下载: 全尺寸图片幻灯片

可见超固结度越高土体最大膨胀力越大，因土体在低吸力时有体缩趋势，当超固结度较低时，最终膨胀力甚至会低于初始围压，但相比膨胀系数C，超固结度的提高对低吸力阶段土体膨胀力降低的影响较小，不同超固结度下土体低吸力阶段膨胀力降低量相近。

2.4 回弹—三轴压缩—湿化—三轴压缩

此路径中先将土样等向压缩后回弹，之后进行三轴压缩，到达某一状态后湿化至饱和，之后再三轴压缩至破坏。首先预测不同膨胀系数对压缩过程的影响，取初始围压p=200 kPa，膨胀系数C分别取为0，2，5，初始吸力为600 kPa，先将土样等向压缩至1000 kPa，然后再回弹至200 kPa，之后进行三轴排水剪切，围压保持200 kPa，当q=200 kPa时，将土样湿化至饱和，再三轴压缩至破坏，分析此压缩过程中体应变与轴向应变的关系，如图 11。

图 11 三轴压缩中湿化不同C值下体应变

Figure 11. Variation of volumetric strain under stress path of triaxial compression-wetting-triaxial compression with different values of C

下载: 全尺寸图片幻灯片

可以看出膨胀系数并不影响湿化前过程，当土体湿化时会发生膨胀，C值越大膨胀量越大，但湿化至饱和后随着三轴压缩的不断进行，体应变会逐渐趋于一致。同时通过分析剪切过程中广义剪应力q的变化，发现随着C值的提高，湿化后土体峰值强度会降低。

其次预测分析不同超固结度的影响，初始吸力仍为600 kPa，膨胀系数C取2，初始围压为200 kPa，将土样分别等向压缩至500，1000，1500 kPa，然后再回弹至200 kPa，之后进行三轴压缩，当q=200 kPa时湿化至饱和，再三轴压缩至破坏，分析压缩过程中体应变与轴向应变之间的关系，如图 12。

图 12 三轴压缩中湿化不同超固结度体应变

Figure 12. Variation of volumetric strain under stress path of triaxial compression-wetting-triaxial compression with different overconsolidation ratios

下载: 全尺寸图片幻灯片

可以看出，超固结度越高，湿化过程中土体膨胀量越大，湿化至饱和后土体三轴剪切阶段体应变的变化趋势与一般非饱和土基本一致。

3. 模型的试验验证

为验证模型的合理性，分别对Zhan等^[28]和Alonso等^[29]的湿化膨胀试验数据进行预测。

3.1 模型的试验验证Ⅰ

Zhan等^[28]对中国湖北枣阳某边坡的膨胀黏土进行了室内试验，土体干密度为1.52~1.60 g/cm³，经测定，土样中砂砾含量3%、粉土颗粒含量58%，黏土颗粒含量39%，液限为50.5%，塑性指数为31，Zhan等^[28]利用非饱和土三轴试验仪对土样进行了多组试验研究，本节取其不同围压下的湿化路径作为验证试验数据。

本文建立的膨胀性非饱和模型需要10个参数，其中力学模型参数 $\lambda ，\kappa ，M，\nu ，N$ 与修正剑桥模型相同，参数 $\lambda ，\kappa ，N$ 可通过饱和土等向压缩回弹试验求得，参数 $M$ 可通过饱和土三轴压缩试验求得，参数 ${\lambda }_{s}，{\kappa }_{s}，{p}^{c}$ 与巴塞罗那模型相同，参数 ${\lambda }_{s}，{\kappa }_{s}$ 可由等围压的干化和湿化试验得到，参数 ${p^{\text{c}}}$ 可通过不同吸力下等向压缩回弹再湿化试验得到，参数 $a$ 可通过不同吸力下的三轴压缩试验得到，反映土体膨胀能力的参数J可根据土体不同围压下的湿化试验求得，根据文献[，]得到参数： $\lambda$ =0.1， $\kappa$ =0.01， ${\lambda _{\text{s}}}$ =0.1， $\nu$ =0.3，M=1.0，N=1.18， ${\kappa _{\text{s}}}$ =0.01， ${p^{\text{c}}}$ =125 kPa，a=150 kPa，J=1.55。大气压 ${p_{{\text{at}}}}$ 取为100 kPa。

此试验取了3组土样，经过吸力平衡，控制土样吸力保持在200 kPa，然后维持吸力不变，分别在40, 70, 100 kPa围压下进行等向固结，在固结完成后，进行湿化试验。湿化过程中保持各自围压不变，吸力从200 kPa降至0 kPa，湿化前土体干密度取1.56 g/cm³，不同围压下湿化过程的试验结果与预测结果对比如图 13所示。

图 13 Zhan等^[28]试验预测对比

Figure 13. Comparison between predicted values and experimental results by Zhan et al. ^[28]

下载: 全尺寸图片幻灯片

可以看出3组土样在试验过程中都发生体积膨胀，但围压不同体积膨胀的幅度不同，围压较大的土样膨胀量较小，预测曲线与试验数据符合较好，可以反映出土体湿化过程中体积膨胀的过程，并体现不同围压下体积膨胀量的不同。

3.2 模型的试验验证Ⅱ

Alonso等^[29]利用Boom膨胀黏土进行了室内试验，其液限为55.9%，塑性指数为29.2，土体的高岭土、伊利石和蒙脱石的含量分别为20%，30%, 10%。试验土样初始干密度为1.40 g/cm³，Alonso等^[29]利用土样进行了不同等向围压下的干湿循环试验，本次验证选取其试验中第一次湿化的试验数据进行验证。

同上一节参数确定方法，本次试验根据文献[，]及参数拟合确定参数： $\lambda$ =0.15， $\kappa$ =0.014， ${\lambda _{\text{s}}}$ =0.1， $\nu$ =0.3，M=1.0，N=1.6， ${\kappa _{\text{s}}}$ =0.008， ${p^{\text{c}}}$ =45 kPa，a=120 kPa，J=1.65。

Chen^[30]将土样分为3组，控制吸力分别为500，500，700 kPa，保持吸力不变，分别在等向压力20，100，400 kPa下等向固结，固结完成后保持围压不变，进行湿化试验，将土样吸力从500 kPa或700 kPa降低到0 kPa，预测土样的体积变化，试验结果与预测结果如图 14。

图 14 Alonso等^[29]试验预测对比

Figure 14. Comparison between predicted values and experimental results by Alonso et al. ^[29]

下载: 全尺寸图片幻灯片

可以看出，由于围压不同，土体在湿化过程中体积变化也有区别，低围压下土体湿化会发生体积膨胀，但当围压增大到足够大时，土体在湿化过程中会发生体积收缩，且随着围压的进一步提高体积收缩也会加剧，这是因为随着围压的提高，土体超固结程度降低，湿化过程中更容易发生屈服，产生更大的湿陷塑性体应变 $\varepsilon _{{\text{vp}}}^{\text{p}}$ ，引发更大的湿陷变形，当围压足够大时这种湿陷变形会占主导地位，使得土体整体表现出体积的收缩，预测曲线可以反映出膨胀土的这一特性，且与试验结果符合较好。

4. 结论

本文提出了考虑膨胀土超固结度和初始干密度的膨胀式，将其融入非饱和土UH模型得到可以考虑理想膨胀性非饱和土的UH模型，该模型具有以下3个特点。

（1）通过引入与土体超固结度和初始干密度有直接关系的膨胀式，可以考虑不同的超固结度和初始干密度对膨胀土应力应变关系的影响，同时可以反映湿化过程中土体超固结度的变化。

（2）以UH模型为基础，并采用硬化参量H，可以体现膨胀土在超固结状态下不同于正常固结状态的应力软化特性和剪胀特性。

（3）本模型在一般非饱和土模型基础上引入的新参数较少，可以通过基础土力学试验获取，模型应用简便。经过模型特性分析与试验结果的验证可以显示出其预测膨胀土应力应变关系的合理性。

另一方面，本模型未考虑水量变化等水力作用对膨胀性非饱和土应力应变关系的影响，需在之后的研究中进一步完善。

图 1 有效应力路径

Figure 1. Effective stress path

下载: 全尺寸图片幻灯片

图 2 SPH方法的计算流程图

Figure 2. Flow chart of SPH method

下载: 全尺寸图片幻灯片

图 3 边坡概化模型

Figure 3. Conceptual model for slope

下载: 全尺寸图片幻灯片

图 4 室内试验测试

Figure 4. Laboratory tests

下载: 全尺寸图片幻灯片

图 5 SPH计算模型

Figure 5. Computational model for SPH

下载: 全尺寸图片幻灯片

图 6 地震动记录

Figure 6. Ground motion recording

下载: 全尺寸图片幻灯片

图 7 有效应力路径对比

Figure 7. Comparison of effective stress paths

下载: 全尺寸图片幻灯片

图 8 剪切应力-应变关系对比

Figure 8. Comparison of shear stress-strain relationship

下载: 全尺寸图片幻灯片

图 9 边坡临界滑移面计算示意图

Figure 9. Calculation diagram of critical slip surface of slope

下载: 全尺寸图片幻灯片

图 10 由Newmark方法计算的滑移面

Figure 10. Slip surface calculated by Newmark method

下载: 全尺寸图片幻灯片

图 11 地震持续10 s时滑移面

Figure 11. Slip surface under earthquake at time 10 s

下载: 全尺寸图片幻灯片

图 12 地震结束后滑移面

Figure 12. Slip surface after earthquake

下载: 全尺寸图片幻灯片

表 1 土体物理力学参数

Table 1 Physical and mechanical parameters of soils

泊松比弹性量/
MPa 重度/
(kN·m^-3) 黏聚力/
kPa 内摩擦角/(°) 残余黏聚力/
kPa 残余内摩擦角/(°)

0.3 100.4 16.8 50.42 15 5.95 15

下载: 导出CSV

参考文献(30)

[1]	彭孔曙, 胡敏云.有限元法在临水边坡设计中的应用[J].科技通报, 2012, 28(9): 142-146. doi: 10.3969/j.issn.1001-7119.2012.09.033 PENG Kongxu, , HU Minyun. Finite element method in the design of slope adjacent to water[J]. Bulletin of Science and Technology, 2012, 28(9): 142-146. (in Chinese) doi: 10.3969/j.issn.1001-7119.2012.09.033
[2]	ISHII Y, OTA K, KURAOKA S, et al. Evaluation of slope stability by finite element method using observed displacement of landslide[J]. Landslides, 2012, 9(3): 335-348. doi: 10.1007/s10346-011-0303-7
[3]	黄帅, 吕悦军.强震作用下动孔隙水压力对砂质边坡动力响应的影响[J].水运工程, 2015(10): 158-167. doi: 10.3969/j.issn.1002-4972.2015.10.028 HUANG Shuai, LYU Yuejun. Influence of dynamic pore water pressure on dynamic response of sandy slope under strong earthquake[J]. Port & Waterway Engineering, 2015(10): 158-167. (in Chinese) doi: 10.3969/j.issn.1002-4972.2015.10.028
[4]	董士杰, 魏红卫.地震作用下土工合成材料加筋土边坡动力分析[J].铁道科学与工程学报, 2015, 12(4): 778-783. doi: 10.3969/j.issn.1672-7029.2015.04.010 DONG Shijie, WEI Hongwei. Dynamic analysis of geosynthetic reinforced soil slope under seismic action[J]. Journal of Railway Science and Engineering, 2015, 12(4): 778-783. (in Chinese) doi: 10.3969/j.issn.1672-7029.2015.04.010
[5]	王飞, 吴红刚, 郭春香.碎石土路堑高边坡地震动力响应过程分析[J].中国地质灾害与防治学报, 2020, 31(1): 18-24. doi: 10.16031/j.cnki.issn.1003-8035.2020.01.03 WANG Fei, WU Honggang, GUO Chunxiang. Dynamic response of high cut based a numerical simulation slope to earthquake[J]. The Chinese Journal of Geological Hazard and Control, 2020, 31(1): 18-24. (in Chinese) doi: 10.16031/j.cnki.issn.1003-8035.2020.01.03
[6]	YIN Y P, LI B, WANG W P. Dynamic analysis of the stabilized Wangjiayan landslide in the Wenchuan Ms 8.0 earthquake and aftershocks[J]. Landslides, 2015, 12(3): 537-547. doi: 10.1007/s10346-014-0497-6
[7]	WU H, ATANGANA NJOCK P G, CHEN J J, et al. Numerical simulation of spudcan-soil interaction using an improved smoothed particle hydrodynamics (SPH) method[J]. Marine Structures, 2019, 66: 213-226. doi: 10.1016/j.marstruc.2019.04.007
[8]	BUI H H, FUKAGAWA R, SAKO K, et al. Slope stability analysis and discontinuous slope failure simulation by elasto-plastic smoothed particle hydrodynamics (SPH)[J]. Géotechnique, 2011, 61(7): 565-574. doi: 10.1680/geot.9.P.046
[9]	HUANG Y, ZHANG W J, XU Q, et al. Run-out analysis of flow-like landslides triggered by the Ms 8.0 2008 Wenchuan earthquake using smoothed particle hydrodynamics[J]. Landslides, 2012, 9(2): 275-283. doi: 10.1007/s10346-011-0285-5
[10]	唐宇峰, 施富强, 廖学燕.基于SPH的边坡稳定性计算中失稳判据研究[J].岩土工程学报, 2016, 38(5): 904-908. http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract16564.shtml TANG Yufeng, SHI Fuqiang, LIAO Xueyan. Failure criteria based on SPH slope stability analysis[J]. Chinese Journal of Geotechnical Engineering, 2016, 38(5): 904-908. (in Chinese) http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract16564.shtml
[11]	NONOYAMA H, MORIGUCHI S, SAWADA K, et al. Slope stability analysis using smoothed particle hydrodynamics (SPH) method[J]. Soils and Foundations, 2015, 55(2): 458-470. doi: 10.1016/j.sandf.2015.02.019
[12]	ZHANG Z Y, JIN X G, BI J. Development of an SPH-based method to simulate the progressive failure of cohesive soil slope[J]. Environmental Earth Sciences, 2019, 78(17): 537. doi: 10.1007/s12665-019-8507-6
[13]	张卫杰, 高玉峰, 黄雨, 等.水土耦合SPH数值模型的正则化修正及其应用[J].岩土工程学报, 2018, 40(2): 262-269. http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract17280.shtml ZHANG Weijie, GAO Yufeng, HUANG Yu, et al. Normalized correction of soil-water-coupled SPH model and its application[J]. Chinese Journal of Geotechnical Engineering, 2018, 40(2): 262-269. (in Chinese) http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract17280.shtml
[14]	张卫杰, 郑虎, 王占彬, 等.基于三维并行SPH模型的土体流滑特性研究[J].工程地质学报, 2018, 26(5): 1279-1284. https://www.cnki.com.cn/Article/CJFDTOTAL-GCDZ201805021.htm ZHANG Weijie, ZHENG Hu, WANG Zhanbin, et al. Study on flowing behavior of soil based on three dimen-sional and parallelized sph model[J]. Journal of Engineering Geology, 2018, 26(5): 1279-1284. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-GCDZ201805021.htm
[15]	DAI Z L, HUANG Y, CHENG H L, et al. 3D numerical modeling using smoothed particle hydrodynamics of flow-like landslide propagation triggered by the 2008 Wenchuan earthquake[J]. Engineering Geology, 2014, 180: 21-33. doi: 10.1016/j.enggeo.2014.03.018
[16]	HUANG Y, DAI Z L. Large deformation and failure simulations for geo-disasters using smoothed particle hydrodynamics method[J]. Engineering Geology, 2014, 168: 86-97. doi: 10.1016/j.enggeo.2013.10.022
[17]	FULK D A. A numerical analysis of smoothed particle hydrodynamics[D]. Wright-Patterson AFB: Thesis Air Force Inst Tech, 1994.
[18]	LIU G R, LIU M B. Smoothed Particle Hydrodynamics-A Meshfree Particle Method[M]. Singapore: World Scientific Publishing Co Pte Ltd, 2003.
[19]	MAO Z, LIU G R, DONG X. A comprehensive study on the parameters setting in smoothed particle hydrodynamics (SPH) method applied to hydrodynamics problems. Computers and Geotechnics, 2017, 92: 77-95. doi: 10.1016/j.compgeo.2017.07.024
[20]	岩本哲也, 小野祐輔.弾性波伝播問題に対する粒子法の適用性[J].応用力学論文集, 2009, 12: 611-622. TETSUYA I, YUSUKE O. Applicability of meshfree particle method to elastic wave propagation analysis[J]. Journal of Applied Mechanics, 2009, 12: 611-622. (in Japanese)
[21]	CHEN J K, BERAUN J E, JIH C J. Completeness of corrective smoothed particle method for linear elastodynamics[J]. Computational Mechanics, 1999, 24(4): 273-285.
[22]	马红权, 张学莹. SPH的核近似和粒子近似[J].信息技术, 2012, 36(7): 170-171,175. https://www.cnki.com.cn/Article/CJFDTOTAL-HDZJ201207049.htm MA Hongquan, ZHANG Xueying. Kernel approximation and particle approximation about SPH[J]. Information Technology, 2012, 36(7): 170-171,175. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-HDZJ201207049.htm
[23]	社本康広, 時松孝次, 有泉浩蔵.一次元有効応力解析の実地盤に対する適用性[J].日本建築学会構造系論文報告集, 1992, 433: 113-119. YDSUHIRO S, KOHJI T, KOUXO A. Applicability of a one-dimensional effective stress analysis to an existing soil deposit[J]. Journal of Structure Construction Engineering, 1992, 433: 113-119. (in Japanese)
[24]	JENNINGS P C. Periodic response of a general yielding structure[J]. Journal of the Engineering Mechanics Division, 1964, 90(2): 131-166.
[25]	GRAY J P, MONAGHAN J J, SWIFT R P. SPH elastic dynamics[J]. Computer Methods in Applied Mechanics and Engineering, 2001, 190(49): 6641-6662.
[26]	HA H, BUI K, SAKO R, et al. Numerical simulation of soil–water interaction using smoothed particle hydrodynamics (SPH) method[J]. Journal of Terramechanics, 2007, 44(5): 339-346.
[27]	冷艺, 栾茂田, 许成顺, 等.应力历史对饱和砂土力学性状影响的试验研究[J].岩土力学, 2009, 30(5): 1257-1263. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX200905012.htm LENG Yi, LUAN Maotian, XU Chengshun, et al. Experimental study of effect of stress history on mechanical properties of saturated sand under complex stress conditions[J]. Rock and Soil Mechanics, 2009, 30(5): 1257-1263. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX200905012.htm
[28]	黄帅, 宋波, 牛立超, 等.地震作用下动孔隙水压力对边坡永久位移影响的简便计算方法[J].建筑结构学报, 2014, 35(3): 215-221. https://www.cnki.com.cn/Article/CJFDTOTAL-JZJB201403028.htm HUANG Shuai, SONG Bo, NIU Lichao, et al. Simple calculation method of permanent displacement of slope influenced by dynamic pore water pressure under earthquake[J]. Journal of Building Structures, 2014, 35(3): 215-221. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-JZJB201403028.htm
[29]	HUANG S, LYU Y J, SHA H J, et al. Seismic performance assessment of unsaturated soil slope in different groundwater levels[J]. Landslides, 2021, 18: 2813-2833.
[30]	霍沿东.基于极限分析上限方法的海底黏性土边坡地震稳定性评[D].大连: 大连理工大学, 2018. HUO Yandong. Evaluation of Seismic Stability of Submarine Clay Slopes Based on Upper Bound Approach of Limit Analysis[D]. Dalian: Dalian University of Technology, 2018. (in Chinese)

施引文献(3)

期刊类型引用(2)

1.	王涛，凡红，王康任，周国庆，王亮亮. 高温冻土双屈服面统一本构模型及其验证. 岩土工程学报. 2025(01): 135-143 . 本站查看
2.	柳伟，徐长节，杜昊东，朱怀龙，王长虹. 基于修正统一硬化模型的超固结非饱和红黏土边坡稳定分析. 岩土力学. 2024(04): 1233-1241 . 百度学术