考虑场地效应的多点地震动作用下边坡永久位移分析

宋健; 潘驭航; 陆朱汐; 姬建; 张飞; 高玉峰

doi:10.11779/CJGE20231282

考虑场地效应的多点地震动作用下边坡永久位移分析 English Version

宋健^{1, 2,},
潘驭航²,
陆朱汐²,
姬建^{1, 2},
张飞^{1, 2},
高玉峰^{1, 2}

1.
河海大学岩土力学与堤坝工程教育部重点实验室, 江苏南京 210024
2.
河海大学土木与交通学院, 江苏南京 210024

基金项目:

国家自然科学基金项目 52378335

国家自然科学基金项目 52178325

国家自然科学基金项目 U22A20594

详细信息

作者简介:
宋健（1988—），男，博士，副教授，主要从事岩土地震工程、边坡工程等方面的研究工作。E-mail：jiansonghh@163.com

中图分类号: TU435
计量
- 文章访问数: 352
- HTML全文浏览量: 64
- PDF下载量: 112
出版历程
- 收稿日期: 2023-12-27
- 网络出版日期: 2024-04-17
- 刊出日期: 2024-12-31

Permanent displacement of slopes under multi-point earthquake ground motions considering site effects

SONG Jian^{1, 2,},
PAN Yuhang²,
LU Zhuxi²,
JI Jian^{1, 2},
ZHANG Fei^{1, 2},
GAO Yufeng^{1, 2}

1.
Key Laboratory of Ministry of Education for Geomechanics and Embankment Engineering, Hohai University, Nanjing 210024, China
2.
College of Civil and Transportation Engineering, Hohai University, Nanjing 210024, China

摘要

摘要: 地震动引起的边坡土层非线性动力响应可能导致坡体不同位置的地震动加速度时程不一致，从而对整个滑动体的极限平衡状态及之后的地震累积永久位移产生影响。为了研究场地效应对边坡地震永久位移的影响，基于极限平衡条分法，通过对不同土条赋予不同的水平和竖向地震动时程，推导出一种考虑场地效应的圆弧和任意形状滑面的多点地震作用下边坡永久位移分析方法。通过与有限差分软件FLAC的数值模拟结果进行对比，验证了该方法能够合理考虑地震动的场地效应。通过研究不同分布形式的多点地震动及竖向地震动对边坡永久位移的影响，结果表明：场地效应引起的水平地震动放大以及多点地震动不同分布模式会导致整个滑体平均地震惯性力作用位置发生改变，从而影响地震作用下滑动体沿滑面的永久变形量，竖向地震动对边坡地震位移影响较小。将方法应用于Lexington土石坝震后变形案例分析，计算得到的坝坡地震永久位移与实际震后观测值相吻合，证明了考虑场地效应的重要性及本文方法的合理性。
- 地震边坡稳定性 /
- 场地效应 /
- 永久位移 /
- 极限平衡条分法 /
- 多点地震动
Abstract: The nonlinear dynamic response of slope soils caused by seismic motion may lead to inconsistent time histories of seismic acceleration at different positions of a slope. This will affect the limit equilibrium state and the cumulative seismic permanent displacement of the slope. In order to investigate the influences of site effects on the permanent displacement of slopes during earthquakes, a method for considering the multi-point earthquake ground motions is presented. The proposed method is derived based on the limit equilibrium slice method by considering different time histories of horizontal and vertical seismic motions for slices and can be used for circular and arbitrary-shaped slip surfaces. The method is capable to reasonably consider the site effects through a comparison with numerical results obtained by the FLAC. The effects of different distributions of multi-point ground motions and vertical ground motions on the permanent displacement of slopes are investigated. The results indicate that the amplification of horizontal ground motions and the different distributions of multi-point ground motions induced by the site effects have a significant impact on seismic slope displacement, while the influence of vertical ground motions are small. The method is then applied to a case study of post-earthquake deformation of the Lexington Dam. The calculated seismic permanent displacement agrees well with the observed post-earthquake values. This confirms the importance of considering the site effects and the rationality of the proposed method in this study.
- seismic slope stability /
- site effect /
- permanent displacement /
- limit equilibrium method /
- multi-point earthquake ground motion

HTML全文

0. 引言

高温冻土是指由土颗粒基质、冰晶、未冻水和气体组成的温度介于-1.5℃~0℃的土体，通常表现出温度敏感性突出、力学性能低劣等特性，工程中极易引发建构筑物冻土地基产生非均匀沉降变形，甚至会导致上部结构开裂和破坏，严重威胁冻土工程安全设计与施工^[1]。因此，研究高温冻土本构模型具有重要工程意义，可为预测高温冻土地层变形特性从而避免工程灾害的发生提供重要手段。

朱元林等^[2]利用大量单轴压缩试验得出冻土应力-应变关系，结合曲线拟合与塑弹性理论建立了单轴压缩本构模型。苗天德等^[3]从微观视角出发，定义了压缩过程中的损伤因子，获得了冻土蠕变损伤本构关系。随后，诸多学者基于损伤力学理论，从细微观角度建立了冻土单轴唯像本构模型^[4]、细观损伤本构模型^[5]、宏细观本构模型^[6]、宏微观黏弹塑性本构模型^[7]。冻土体常在低围压条件下呈现出应变软化的变形特征，通过CT技术识别冻土内部变形过程中的损伤变量，可以描述冻土体蠕变过程中的损伤特征^[8-9]。对不同温度不同应变率下的冻土试件进行冲击压缩试验，可构建动荷载条件下的冻土本构关系^[10-11]。将加载-破坏屈服面及统一硬化参数整合，可建立冻土多相多面本构模型^[12]。空心圆柱仪试验可以科学量化不同应力水平条件下的冻土变形特征，为复杂应力条件下的本构关系构建提供了重要参考^[13]。

对于常规低温冻结，上述研究成果得到的本构模型基本能表征冻土不同应力状态下的力学状态^[14-15]。然而对于近相变区的高温冻土，冻结过程中冻结锋面、冰分凝锋面、冻结缘厚度随时间动态变化^[16-17]。当高温冻土未冻水含量相对较高时，内部的孔隙水压力增加较大使得作用在土骨架上的有效应力降低，表现出剪胀、应变软化特征；当高温冻土未冻水含量相对较低时，内部的孔隙水压力增加值小于总应力增加值，表现出剪缩、应变硬化特征，因此高温冻土同时存在剪缩、剪胀、硬化和软化特性^[18-19]。针对高温冻土，Liao等^[20]从强度随机分布出发，建立了可以反映冻土力学性质较强离散性的随机损伤本构模型，但将其作为弹性体的损伤本构模型研究存在一定的缺陷。Lai等^[21]从概率论与数理统计的视角出发，定义了材料截面内部的损失变量，提出了高温冻土损伤统计本构模型，但该模型难以对高低围压三轴应力路径下的高温冻土应力应变行为进行统一描述。

本文考虑黏聚力及内摩擦角的影响，以整体变形 $\varepsilon _{\text{v}}^{}$ -lnp曲线对冻土试样变形特征进行表征，建立了能够有效反映高温冻土剪缩、剪胀、硬化和软化特性的本构模型。该模型采用应力路径相关因子对等向固结路径条件下屈服面硬化参数进行修正，获得与应力路径无关的当前屈服面硬化参数，结合临界状态参数推导得到了能够体现高温冻土固有临界特性的参考屈服面硬化参数。自编程序计算了不同围压条件下高温冻土应力应变关系，利用试验数据对本文所构建的双屈服面统一本构模型进行了验证。

1. 本构模型建立

1.1 屈服函数

英国剑桥大学Roscoe等^[22]从能量方程推导出应力比与应变增量比的关系，结合正交法则、相关联流动法则及等向固结试验确定塑性应力应变关系，于1963年提出了原始剑桥模型。鉴于原始剑桥模型获得的屈服函数在 $p{\text{ - }}q$ 坐标系中与 $p$ 轴不正交，Roscoe等^[23]于1968年提出了塑性势函数(屈服函数)与 $p$ 轴正交的修正剑桥模型，定义了应力空间内 $g(p,q)$ 和 ${\text{(d}}\varepsilon _{\text{v}}^{\text{p}}{\text{,d}}\varepsilon _{\text{d}}^{\text{p}}{\text{)}}$ 的正交条件，表示为

。

${\text{d}}p{\text{d}}\varepsilon _{\text{v}}^{\text{p}} + {\text{d}}q{\text{d}}\varepsilon _{\text{d}}^{\text{p}} = 0 。$

(1)

荷载作用下塑性变形所消耗的功为

。

${\text{d}}{W^{\text{p}}} = p{\text{d}}\varepsilon _{\text{v}}^{\text{p}} + q{\text{d}}\varepsilon _{\text{d}}^{\text{p}} 。$

(2)

考虑到材料破坏时，也即临界状态下有

。

${q_{\text{f}}} = M(p + {p_r})\;\;\;({\text{d}}\varepsilon _{\text{v}}^{\text{p}} = 0) 。$

(3)

式中：M = $6\sin \varphi /(3 - \sin \varphi )$ ； ${p_{\text{r}}} = c\cot \varphi$ ，c， $\varphi$ 分别为土的黏聚力和内摩擦角。

将式（3）代入式（2）可得

$\begin{aligned} \mathrm{d} W^{\mathrm{p}} & =\left(p+p_{\mathrm{r}}\right) \mathrm{d} \varepsilon_{\mathrm{v}}^{\mathrm{p}}+q \mathrm{~d} \varepsilon_{\mathrm{d}}^{\mathrm{p}}=M\left(p+p_{\mathrm{r}}\right) \mathrm{d} \varepsilon_{\mathrm{d}}^{\mathrm{p}} \\ & =\left(p+p_{\mathrm{r}}\right) \sqrt{\left(\mathrm{d} \varepsilon_{\mathrm{v}}^{\mathrm{p}}\right)^2+\left(M \mathrm{~d} \varepsilon_{\mathrm{d}}^{\mathrm{p}}\right)^2} 。\end{aligned}$

(4)

整理式（4）可得

$\frac{{{\text{d}}\varepsilon _{\text{v}}^{\text{p}}}}{{{\text{d}}\varepsilon _{\text{d}}^{\text{p}}}} = \frac{{{M^2}{{(p + {p_{\text{r}}})}^2} - {q^2}}}{{2(p + {p_{\text{r}}})q}} 。$

(5)

求解方程（5）可得塑性势函数为

$g = {q^2} + {M^2}{(p + {p_{\text{r}}})^2} - C(p + {p_{\text{r}}}) = 0 。$

(6)

式中： $C$ 为积分常数， $C$ 取不同的值可得到不同的塑性势函数。

本文得到考虑高温冻土黏聚力及内摩擦角的修正剑桥模型塑性势面如图 1所示。

图 1 塑性势函数

Figure 1. Plastic potential function

下载: 全尺寸图片幻灯片

采用相关流动法则，当广义剪应力 $q = 0$ 时，令平均正应力 $p = {p_x}$ ，积分常数 $C = {M^2}({p_x} + {p_r})$ ，高温冻土屈服函数可进一步表示为

$f = {q^2} + {M^2}p(p + {p_r}) - {M^2}{p_x}(p + {p_r}) = 0 。$

(7)

变换式（7），并对两边取自然对数可得

$\ln \left[ {1 + \frac{{{q^2}}}{{{M^2}p(p + {p_r})}}} \right] = \ln {p_x} - \ln p 。$

(8)

1.2 屈服面参量及其修正

将高温冻土的等向固结试验（ ${\sigma _1} = {\sigma _2} = {\sigma _3}$ ）结果绘制在 ${\varepsilon _{\text{v}}}{\text{ - }}\ln p$ 坐标系中，如所示，回弹再加载曲线的斜率为 $\kappa$ ，记等向固结曲线NCL为

${\varepsilon _{\text{v}}} = \varphi \ln p 。$

(9)

图 2 ε_v-lnp平面中等向固结线NCL

Figure 2. Normal consolidation line in ε_v-lnp plane

下载: 全尺寸图片幻灯片

对式（9）进行微分可得

${\text{d}}{\varepsilon _{\text{v}}} = \varphi '(\ln p){\text{d}}\ln p 。$

(10)

弹性体积应变变化量为

${\text{d}}\varepsilon _{\text{v}}^{\text{e}} = \kappa d\ln p 。$

(11)

综合式（10），（11），得到塑性体积应变变化量为

${\text{d}}\varepsilon _{\text{v}}^{\text{p}} = {\text{d}}{\varepsilon _{\text{v}}} - {\text{d}}\varepsilon _{\text{v}}^{\text{e}} = (\varphi '\ln p - \kappa )d\ln p 。$

(12)

对式（12）进行变换并积分得到

$\int {\frac{{{\text{d}}\varepsilon _{\text{v}}^{\text{p}}}}{{\varphi '(\ln p) - \kappa }}} = \ln {p_x} - \ln {p_0} 。$

(13)

式中： ${p_0}$ 为与初始体积应变 ${\varepsilon _{{\text{v0}}}}$ 相对应的当前屈服面上的球应力，其决定了初始当前屈服面与 $p$ 轴的交点，其数值的确定在模型参数确定方法中详细介绍。

将式（13）代入式（8）得到屈服面函数为

$f = \ln \frac{p}{{{p_0}}} + \ln \left( {1 + \frac{{{q^2}}}{{{M^2}p(p + {p_r})}}} \right) - \int {\frac{{{\text{d}}\varepsilon _{\text{v}}^{\text{p}}}}{{\varphi '(\ln p) - \kappa }}} = 0 。$

(14)

根据式（14），高温冻土屈服面硬化参量可定义为

$H = \int {\frac{{{\text{d}}\varepsilon _{\text{v}}^{\text{p}}}}{{\varphi '(\ln p) - \kappa }}} 。$

(15)

为了使所要构建的高温冻土本构模型能够反映常压下的软化剪胀特性及塑性体积应变由剪缩阶段发展至剪胀阶段的特征状态分界线，模型中的硬化参量须包含能够有效反映特征状态的参数。根据姚仰平等^[24]提出的应力路径相关因子的方法，将应力路径相关因子设为

$R(\eta ) = \frac{{M_{\text{f}}^4}}{{{M^4}}}\frac{{{M^4} - {\eta ^4}}}{{M_{\text{f}}^4 - {\eta ^4}}} 。$

(16)

依据式（16）表征的应力路径相关因子，屈服面硬化参量修正为

$H = \int {\frac{1}{{R(\eta )}}\frac{{{\text{d}}\varepsilon _{\text{v}}^{\text{p}}}}{{\varphi '(\ln p) - \kappa }}} = \int {\frac{{{M^4}(M_{\text{f}}^4 - {\eta ^4})}}{{M_{\text{f}}^4({M^4} - {\eta ^4})}}\frac{{{\text{d}}\varepsilon _{\text{v}}^{\text{p}}}}{{\varphi '(\ln p) - \kappa }}} 。$

(17)

由此得到修正后的高温冻土当前屈服面函数：

$f=\mathrm{ln}\frac{p}{{p}_{0}}+\mathrm{ln}\left(1+\frac{{q}^{2}}{{M}^{2}p(p+{p}_{\text{r}})}\right)-{\displaystyle \int \frac{{M}^{4}({M}_{\text{f}}^{4}-{\eta }^{4})}{{M}_{\text{f}}^{4}({M}^{4}-{\eta }^{4})}\frac{\text{d}{\epsilon }_{\text{v}}^{\text{p}}}{{\phi }^{\prime }(\mathrm{ln}p)-\kappa }}。$

(18)

临界状态线在 $q{\text{ - }}p$ 空间中表示为一条不可逾越的应力比直线状态，在 ${\varepsilon _{\text{v}}}{\text{ - }}\ln p$ 空间中表示为一个稳定的非线性状态。通过三轴压缩试验可获得不同围压条件下由剪切开始至临界状态整个应力路径过程中冻土材料体积应变改变量与球应力p的关系，对达到临界状态时的球应力p取自然对数并统计与其相对应的体积应变可获得 ${\varepsilon _{\text{v}}}{\text{ - }}\ln p$ 空间中的CSL曲线，如图 3所示，记临界状态曲线CSL为

${\varepsilon _{\text{v}}} = \psi (\ln \bar p) 。$

(19)

图 3 ε_v-lnp平面中临界状态线CSL

Figure 3. Critical state line in ε_v-lnp plane

下载: 全尺寸图片幻灯片

可以看出，图 3中的高温冻土临界状态线与图 2中的等向压缩曲线形式相同，因此对其硬化参量的推导同式（13），可以得到

$\int {\frac{{{\text{d}}\varepsilon _{\text{v}}^{\text{p}}}}{{\psi '(\ln \bar p) - \kappa }}} = \ln {\bar p_x} - \ln {\bar p_0} 。$

(20)

式中： ${\bar p_0}$ 为与初始体积应变 ${\varepsilon _{{\text{v0}}}}$ 相对应的参考屈服面上的球应力，其决定了初始参考屈服面与 $p$ 轴的交点，其数值的确定在模型参数确定方法中详细介绍。

根据式（20）以及屈服函数式（8）便可推导得出参考屈服面函数为

$f = \ln \frac{{\bar p}}{{{{\bar p}_0}}} + \ln \left[ {1 + \frac{{{{\bar q}^2}}}{{{M^2}\bar p(\bar p + {p_{\text{r}}})}}} \right] - \int {\frac{{{\text{d}}\varepsilon _{\text{v}}^{\text{p}}}}{{\psi '(\ln \bar p) - \kappa }}} = 0 。$

(21)

参考屈服面硬化参量为

$\bar H = \int {\frac{{{\text{d}}\varepsilon _{\text{v}}^{\text{p}}}}{{\psi '\ln \bar p - \kappa }}} 。$

(22)

1.3 弹性应变增量

对于高温冻土本构模型的弹性部分，根据各向同性的假设，弹性应变增量可用张量表示为

${\text{d}}\varepsilon _{ij}^{\text{e}} = \frac{{1 + \nu }}{E}{\text{d}}{\sigma _{ij}} - \frac{\nu }{E}{\text{d}}{\sigma _{kk}}{\delta _{ij}} 。$

(23)

式中：E为弹性模量； $\nu$ 为泊松比。

根据标记法规则， ${\text{d}}{\sigma _{kk}} = {\text{d}}{\sigma _{11}} + {\text{d}}{\sigma _{22}} + {\text{d}}{\sigma _{33}}$ ，因此高温冻土体积应变与剪切应变增量形式为

${\text{d}}\varepsilon _{\text{v}}^{\text{e}} = \frac{{3(1 - 2\nu )}}{E}{\text{d}}p{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\text{d}}\varepsilon _{\text{d}}^{\text{e}} = \frac{{2(1 + \nu )}}{{3E}}{\text{d}}q 。$

(24)

将体积变形模型与剪切变形模量作为其计算参数，通过高温冻土等向固结试验（ ${\sigma _1} = {\sigma _2} = {\sigma _3}$ ），将结果绘制在 ${\varepsilon _{\text{v}}}{\text{ - }}\ln p$ 坐标系中，得到加载后回弹再加载曲线斜率 $\kappa$ ，则弹性体积应变增量公式为

${\text{d}}\varepsilon _v^e = \kappa d\ln p = \frac{\kappa }{p}{\text{d}}p 。$

(25)

对比式（24），（25）的数学表达计算式，模量K与模量G可进一步表示为

$K = \frac{{{\text{d}}p}}{{{\text{d}}\varepsilon _{\text{v}}^{\text{e}}}} = \frac{E}{{3(1 - 2\nu )}} = \frac{p}{\kappa } \text{，}$

(26)

$G = \frac{{{\text{d}}q}}{{3{\text{d}}\varepsilon _{\text{d}}^{\text{e}}}} = \frac{E}{{2(1 + \nu )}} = \frac{{3(1 - 2\nu )}}{{2(1 + \nu )}}K 。$

(27)

1.4 塑性应变增量

高温冻土在复杂应力状态下发生变形时，总应变可表示为

${\text{d}}{\varepsilon _{ij}} = {\text{d}}\varepsilon _{ij}^{\text{e}} + {\text{d}}\varepsilon _{ij}^{\text{p}} 。$

(28)

应力应变关系为

${\text{d}}{\sigma _{ij}} = C_{ijkl}^{\text{e}}{\text{d}}\varepsilon _{kl}^{\text{e}} = C_{ijkl}^{\text{e}}({\text{d}}{\varepsilon _{kl}} - {\text{d}}\varepsilon _{kl}^{\text{p}}) = C_{ijkl}^{{\text{ep}}}{\text{d}}{\varepsilon _{kl}} 。$

(29)

式中： $C_{ijkl}^{\text{e}}$ 为弹性张量； $C_{ijkl}^{{\text{ep}}}$ 为弹塑性张量。

由广义胡克定律可知弹性张量 $C_{ijkl}^{{\text{ep}}}$ 的具体表达式为

${C}_{ijkl}^{\text{e}}=\left(K-\frac{2}{3}G\right){\delta }_{ij}{\delta }_{kl}+G({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk})\text{ }。$

(30)

塑性应变增量由塑性位势理论确定, 因此存在如下数学关系：

${\text{d}}\varepsilon _{ij}^{\text{p}} = {\text{d}}\lambda \frac{{\partial g}}{{\partial {\sigma _{ij}}}} 。$

(31)

式中： ${\text{d}}\lambda$ 为塑性标量因子。

将屈服函数f简记为

$f({\sigma _{ij}},H(\varepsilon _{ij}^{\text{p}})) = 0 。$

(32)

对式（32）微分可得

${\text{d}}f = \frac{{\partial f}}{{\partial {\sigma _{ij}}}}{\text{d}}{\sigma _{ij}} + \frac{{\partial f}}{{\partial H}}\frac{{\partial H}}{{\partial \varepsilon _{ij}^{\text{p}}}}{\text{d}}\varepsilon _{ij}^{\text{p}} = 0 。$

(33)

将式（29），（31）代入式（33），并整理可得塑性标量因子 ${\text{d}}\lambda$ 的具体表达式为

${\text{d}}\lambda = \frac{{\frac{{\partial f}}{{\partial {\sigma _{ij}}}}C_{ijkl}^{\text{e}}{\text{d}}{\varepsilon _{kl}}}}{{ - \frac{{\partial f}}{{\partial H}}\frac{{\partial H}}{{\partial \varepsilon _{ij}^{\text{p}}}}\frac{{\partial g}}{{\partial {\sigma _{ij}}}} + \frac{{\partial f}}{{\partial {\sigma _{ij}}}}C_{ijkl}^{\text{e}}\frac{{\partial g}}{{\partial {\sigma _{kl}}}}}} 。$

(34)

将式（31），（34）代入式（29），并整理可得弹塑性张量 $C_{ijkl}^{{\text{ep}}}$ 的具体表达式为

$C_{ijkl}^{{\text{ep}}} = C_{ijkl}^{\text{e}} - \frac{{C_{ijmn}^{\text{e}}\frac{{\partial g}}{{\partial {\sigma _{mn}}}}\frac{{\partial f}}{{\partial {\sigma _{st}}}}C_{stkl}^{\text{e}}}}{{A + \frac{{\partial f}}{{\partial {\sigma _{ij}}}}C_{ijkl}^{\text{e}}\frac{{\partial g}}{{\partial {\sigma _{kl}}}}}} 。$

(35)

式中： $A = - \frac{{\partial f}}{{\partial H}}\frac{{\partial H}}{{\partial \varepsilon _{ij}^{\text{p}}}}\frac{{\partial g}}{{\partial {\sigma _{ij}}}}$ 。

1.5 应力应变关系

采用相关流动法则，只需对屈服函数f进行微分求导，根据复合函数求导法则有

$\frac{{\partial f}}{{\partial \sigma }} = \frac{{\partial f}}{{\partial {I_1}}}\frac{{\partial {I_1}}}{{\partial \sigma }} + \frac{{\partial f}}{{\partial \sqrt {{J_2}} }}\frac{{\partial \sqrt {{J_2}} }}{{\partial \sigma }} + \frac{{\partial f}}{{\partial \theta }}\frac{{\partial \theta }}{{\partial \sigma }} 。$

(36)

式中： ${I_1}$ 为第一主应力不变量； ${J_2}$ 为第二偏应力不变量； $\theta$ 为Lode角。

$\frac{{\partial \theta }}{{\partial \sigma }} = \frac{{\sqrt 3 }}{{2sin3\theta }}\left[ {\frac{1}{{{{({J_2})}^{{3 \mathord{\left/ {\vphantom {3 2}} \right. } 2}}}}}\frac{{\partial {J_3}}}{{\partial \sigma }} - \frac{{3{J_3}}}{{{{({J_2})}^2}}}\frac{{\partial \sqrt {{J_2}} }}{{\partial \sigma }}} \right] 。$

(37)

将式（37）代入式（36），整理可得

$\begin{gathered} \frac{\partial f}{\partial \sigma}=\frac{\partial f}{\partial I_1} \frac{\partial I_1}{\partial \sigma}+\left[\frac{\partial f}{\partial \sqrt{J_2}}+\frac{\partial f}{\partial \theta} \frac{\cot 3 \theta}{\sqrt{J_2}}\right] \frac{\partial \sqrt{J_2}}{\partial \sigma}+\frac{\partial f}{\partial \theta} \frac{\sqrt{3}}{2 \sin 3 \theta} . \\ \frac{1}{\left(J_2\right)^{3 / 2}} \frac{\partial J_3}{\partial \sigma} \end{gathered} 。$

(38)

考虑高温冻土当前屈服面函数形式的计算公式（14），（15），根据复合函数求导法则有

$\begin{gathered} \frac{\partial f}{\partial I_1}=\frac{\sqrt{3}}{2 \sin 3 \theta} \frac{1}{\left(J_2\right)^{3 / 2}} \frac{-2 q^2}{M(\theta)^3 p\left(p+p_r\right)+M(\theta) q^2}. \\ \frac{\partial M(\theta)}{\partial \theta} \end{gathered} 。$

(39)

由屈服函数中临界应力比 $M(\theta )$ 定义有

$\left.\begin{array}{l} M(\theta)=M\left(\theta=0^{\circ}\right)\left[\frac{2 \alpha^4}{1+\alpha^4-\left(1-\alpha^4\right) \cos 3 \theta}\right]^{\frac{1}{4}}, \\ \alpha=\frac{M\left(\theta=60^{\circ}\right)}{M\left(\theta=0^{\circ}\right)}=\frac{3-\sin \varphi}{3+\sin \varphi} 。 \end{array}\right\}$

(40)

对式（40）微分得

$\begin{aligned} \frac{\partial M(\theta)}{\partial \theta}= & {\left[\frac{2 \alpha^4}{1+\alpha^4-\left(1-\alpha^4\right) \cos 3 \theta}\right]^{\frac{5}{4}} } \\ & \frac{3\left(\alpha^4-1\right)}{8 \alpha^4} \sin 3 \theta \cdot M\left(\theta=0^{\circ}\right) 。\end{aligned}$

(41)

将式（40），（41）代入式（38）即可求得屈服函数f任一应力分量的一阶导数值。

根据式（16），（17）可得

$\begin{aligned} \frac{\partial f}{\partial H} \frac{\partial H}{\partial \varepsilon_{i j}^{\mathrm{p}}} \frac{\partial g}{\partial \sigma_{i j}} & =-\frac{M(\theta)^4\left(M(\theta)_{\mathrm{f}}^4-\eta^4\right)}{M(\theta)_{\mathrm{f}}^4\left(M(\theta)^4-\eta^4\right)} \frac{1}{\varphi^{\prime}(\ln p)-\kappa} .\\ & \frac{\partial \varepsilon_{\mathrm{v}}^{\mathrm{p}}}{\partial \varepsilon_{i j}^{\mathrm{p}}} \frac{\partial g}{\partial \sigma_{i j}} 。\end{aligned}$

(42)

因为 $\varepsilon _{\text{v}}^{\text{p}} = \varepsilon _{11}^{\text{p}} + \varepsilon _{22}^{\text{p}} + \varepsilon _{33}^{\text{p}}$ ，所以存在

$\left[ {\frac{{\partial \varepsilon _{\text{v}}^{\text{p}}}}{{\partial \varepsilon _{ij}^{p} }}} \right] = \left[ {1\;\;1\;\;1\;\;0\;\;0\;\;0} \right] 。$

(43)

将式（43）代入式（42），整理可得

$A=\frac{M{(\theta )}^{4}(M{(\theta )}_{\text{f}}^{4}-{\eta }^{4})}{M{(\theta )}_{\text{f}}^{4}(M{(\theta )}^{4}-{\eta }^{4})}\frac{1}{{\phi }^{\prime }(\mathrm{ln}p)-\kappa }\left(\frac{\partial g}{\partial {\sigma }_{11}}+\frac{\partial g}{\partial {\sigma }_{22}}+\frac{\partial g}{\partial {\sigma }_{33}}\right)。$

(44)

将式（30），（44）代入式（35）即可求得弹塑性张量 $C_{ijkl}^{{\text{ep}}}$ 的具体表达式。

1.6 模型参数及确定方法

根据高温冻土强度包线（直线）可确定内摩擦角 $\varphi$ 和黏聚力 $c$ ；泊松比 $\nu$ 可通过 $- \Delta {\varepsilon _r}/\Delta {\varepsilon _z}$ 获得；对于等向固结试验 ${\varepsilon _{\text{v}}}$ -lnp空间中回弹线斜率 $\kappa$ 可通过加载后回弹再加载试验确定； ${\varepsilon _{{\text{v0}}}}$ 为压缩固结完成开始剪切所对应的体积应变；伏斯列夫面(Hvorslev)斜率 ${M_{\text{h}}}$ 是 $q{\text{ - }}p$ 空间中不排水强度包线的斜率，通过冻土等向固结试验 $({\sigma _1} = {\sigma _2} = {\sigma _3})$ ，将结果绘制在 ${\varepsilon _{\text{v}}}{\text{ - }}p$ 坐标系中，便可获得等向固结曲线NCL，与初始体积应变 ${\varepsilon _{{\text{v0}}}}$ 对应的即为 ${p_0}$ 。通过多组不同压力等级下的三轴压缩试验，获得不同压力条件下的冻土临界体积应变与临界状态时球应力p的关系，将结果绘制于 ${\varepsilon _{\text{v}}}{\text{ - }}p$ 空间中，获得其临界状态曲线CSL，与初始体积应变 ${\varepsilon _{{\text{v0}}}}$ 对应的即为 ${\bar p_0}$ 。

2. 模型特性分析

2.1 固结参数与潜在强度

由式（21）推导可得参考应力点的球应力 $\bar p$ 为

$\bar p = {\bar p_0}\left( {\frac{{{M^2}}}{{{M^2} + \eta \eta '}}} \right)\exp \left( {\int {\frac{{{\text{d}}\varepsilon _{\text{v}}^{\text{p}}}}{{\psi '(\ln \bar p) - \kappa }}} } \right) 。$

(45)

式中： $\eta = q/(p + {p_r}) = \bar q/(\bar p + {p_r})$ ， $\eta ' = \bar q/\bar p$ 。

依据平面解析几何中直线与椭圆的基本关系，结合方程式（21）及 $\eta$ ，可求得 $\eta '$ 的具体表达式为

$\eta ' = \frac{{({M^2} + {\eta ^2})\eta {p_r}}}{{{M^2}{{\bar p}_0}\exp \left( {\int {\frac{{{\text{d}}\varepsilon _{\text{v}}^{\text{p}}}}{{\psi '(\ln \bar p) - \kappa }}} } \right) - {\eta ^2}{p_r}}} + \eta 。$

(46)

参照超固结土概念，定义高温冻土固结参数 $R$ 为参考应力点的球应力 $\bar p + {p_r}$ 与当前应力点的球应力 $p + {p_r}$ 的比值，即

$R = \frac{{\bar p + {p_r}}}{{p + {p_r}}} 。$

(47)

综合式（45），（46）可得高温冻土固结参数 $R$ 的具体表达式为

$R = \frac{1}{{p + {p_r}}}\left\{ {{{\bar p}_0}\left( {\frac{{{M^2}}}{{{M^2} + \eta \eta '}}} \right)\exp \left[ {\int {\frac{{{\text{d}}\varepsilon _{\text{v}}^{\text{p}}}}{{\psi '(\ln \bar p) - \kappa }}} } \right] + {p_r}} \right\} 。$

(48)

将高温冻土当前屈服面、参考屈服面、伏斯列夫面(Hvorslev)、临界状态线CSL绘制于同一坐标系（p-q）中，如所示。对于特定的应力状态点 $(p,q)$ ，潜在强度M_f与其相对应，并由参考应力点 $(\bar p,\bar q)$ 与冻土材料伏斯列夫面(Hvorslev)斜率共同决定，剪应力 ${q_{\text{f}}}$ 可表示为

${q_{\text{f}}} = M(\bar p + {p_r}) - {M_{\text{h}}}(\bar p - p) \text{，}$

(49)

${M_{\text{f}}} = \frac{{M(\bar p + {p_r}) - {M_{\text{h}}}(\bar p - p)}}{{p + {p_r}}} = M + (R - 1)(M - {M_{\text{h}}}) 。$

(50)

图 4 基于Hvorslev面的潜在强度M_f

Figure 4. Potential strength M_f based on Hvorslev surface

下载: 全尺寸图片幻灯片

图 5表明了高温冻土固结参数R、潜在强度M_f与硬化参数H之间的相互关系，由于不同的固结程度的高温冻土具有不同的潜在强度，所以固结参数R决定潜在强度M_f的大小；式（17）的硬化参数H中包含了潜在强度M_f，故M_f的变化直接影响硬化参数H；高温冻土在复杂应力作用下的硬化参数H决定着当前屈服面与参考屈服面状态及表达式，进而间接影响固结参数R。

图 5 高温冻土模型参数R，M_f，H相互关系

Figure 5. Relationship among model parameters R, M_f and H of warm frozen soil

下载: 全尺寸图片幻灯片

2.2 模型演化分析

微分当前屈服面硬化参量式（15）及参考屈服面硬化参量式（22）可得

${\text{d}}H = \frac{{{M^4}(M_{\text{f}}^4 - {\eta ^4})}}{{M_{\text{f}}^4({M^4} - {\eta ^4})}}\frac{{{\text{d}}\varepsilon _{\text{v}}^{\text{p}}}}{{\varphi '(\ln p) - \kappa }} \text{，}$

(51)

${\text{d}}\bar H = \frac{{{\text{d}}\varepsilon _{\text{v}}^{\text{p}}}}{{\psi '(\ln p) - \kappa }} 。$

(52)

对式（51），（52）进行整个应力路径中的分析，可得表 1所示的变化规律。

表 1 模型状态演化过程

Table 1. Evolution process of model states

状态阶段	应力比η、特征状态应力比M及潜在强度 ${M_{\text{f}}}$ 间关系	塑性体变增量 ${\text{d}}\varepsilon _{\text{v}}^{\text{p}}$ 、屈服面硬化参量增量增量 ${\text{d}}H$ 、 ${\text{d}}\overline H$
初始状态	$0 = \eta < M < {M_{\text{f}}}$	${\text{d}}\varepsilon _{\text{v}}^{\text{p}} > 0\;,\;{\text{d}}H > 0\;,\;{\text{d}}\bar H > 0$
剪缩、硬化阶段	$0 < \eta < M < {M_{\text{f}}}$	${\text{d}}\varepsilon _{\text{v}}^{\text{p}} > 0\;,\;{\text{d}}H > 0\;,\;{\text{d}}\bar H > 0$
特征状态	$0 < \eta = M < {M_{\text{f}}}$	${\text{d}}\varepsilon _{\text{v}}^{\text{p}} = 0\;,\;{\text{d}}H > 0\;,\;{\text{d}}\bar H = 0$
剪胀、硬化阶段	$0 < M < \eta < {M_{\text{f}}}$	${\text{d}}\varepsilon _{\text{v}}^{\text{p}} < 0\;,\;{\text{d}}H > 0\;,\;{\text{d}}\bar H < 0$
峰值应力状态	$0 < M < \eta = {M_{\text{f}}}$	${\text{d}}\varepsilon _{\text{v}}^{\text{p}} < 0\;,\;{\text{d}}H = 0\;,\;{\text{d}}\bar H < 0$
第一剪胀软化阶段	$0 < M < {M_{\text{f}}} < \eta$	${\text{d}}\varepsilon _{\text{v}}^{\text{p}} < 0\;,\;{\text{d}}H < 0\;,\;{\text{d}}\bar H < 0$
第二剪胀软化阶段	$0 < M = {M_{\text{f}}} < \eta$	${\text{d}}\varepsilon _{\text{v}}^{\text{p}} < 0\;,\;{\text{d}}H < 0\;,{\text{d}}\bar H < 0$
临界状态	$0 < M = {M_{\text{f}}} = \eta$	${\text{d}}\varepsilon _{\text{v}}^{\text{p}} = 0\;,\;{\text{d}}H = 0\;,\;{\text{d}}\bar H = 0$

下载: 导出CSV

| 显示表格

由此可见，建立在潜在强度M_f、特征状态应力比M和应力比η有关的硬化参数H的基础上，本文模型能够有效反映高温冻土的剪缩、剪胀、硬化和软化特性。

3. 模型试验验证

3.1 模型的试验验证Ⅰ

Ma等^[17]对兰州某地饱和冻结粉土进行了一系列三轴剪切室内试验，粉土干密度为1.75 g/cm³，三轴试验MTS810围压设定分别为1.0，2.0，6.0 MPa，高温冻结温度为-1.5℃；试验过程分为等向加载阶段和轴向应力加载阶段，围压约束的加载在5 min内完成，轴向应力以1.108 mm/min的速度加载，以试样轴向应变达到15%时作为试样破坏条件；本文所建模型需要6个参数和2条试验曲线，根据文献[，]及参数拟合： $c$ =112.3 kPa， $\varphi$ =25.8℃， $\nu$ =0.32， $\kappa$ =-8.8×10^-7， ${\varepsilon _{{\text{v0}}}}$ =1.72×10^-2， ${M_{\text{h}}}$ =1.17，等向固结曲线曲线 ${\varepsilon _{\text{v}}} =$ ${\text{2}}{\text{.18}} \times {\text{1}}{{\text{0}}^{ - {\text{6}}}} \cdot {\text{exp(0}}{\text{.726ln}}p{\text{)}}$ ，临界状态曲线曲线 ${\varepsilon _{\text{v}}} =$ ${\text{1}}{\text{.83}} \times {\text{1}}{{\text{0}}^{ - {\text{5}}}} \cdot {\text{exp(0}}{\text{.687ln}}p{\text{)}}$ ，

图 6（a）为-1.5℃条件下高温冻结粉土偏应力与轴向应变的模型预测结果与试验数据对比，可以看出，本文模型对1.0，2.0，6.0 MPa围压条件下剪切行为的预测结果与试验数据吻合度较高。在1.0 MPa围压条件下，高温冻结粉土偏应力一直增加（硬化），在6.0 MPa围压条件下，高温冻结粉土偏应力先增加后减小（软化）；在相同围压下，应力应变曲线主要呈三阶段变化，在初始阶段曲线呈线性变化，随着轴向应变的增加，曲线逐渐过渡到弹塑性阶段，此时试样内部以不可恢复的塑性变形为主。对于应变硬化的情况，应力逐渐到达最大值，曲线斜率近视趋近于零；对于应变软化的情况，当偏应力超过一定值时，曲线斜率由正转为负值，这是由于高温冻土特有的压融现象使得未冻水含量增加，试样内部超静孔隙水压力增加，颗粒之间的有效应力降低，剪切破坏极易发生。图 6（b）为-1.5℃条件下体积应变与轴向应变的模型预测结果与试验数据对比，可以看出，本文模型对1.0，2.0，6.0 MPa围压条件下体变行为的预测结果与试验数据吻合度较好。在围压较低时，试样因侧向变形受到的约束力较小而表现为稳定的剪缩特征；在围压较高时，试样因内部冰晶结构微裂隙及孔隙损伤开裂破坏而表现为特定的剪胀特征。综合图 6（a）和图 6（b）的试验数据和预测曲线，可以得出，-1.5℃条件下的高温冻结粉土在低围压条件下总体表现出剪缩硬化特性，在高围压条件下总体表现出剪胀软化特性，本文所建模型可以有效反映-1.5℃条件下的高温冻土剪缩、剪胀、硬化和软化特性。

图 6 模型预测结果与试验数据对比（-1.5℃）

Figure 6. Comparison between model results with experimental data (-1.5℃)

下载: 全尺寸图片幻灯片

3.2 模型的试验验证Ⅱ

宋丙堂等^[18]对兰州市安宁区某地冻结粉土进行了一系列三轴压缩加卸载试验，粉土干密度为1.70 g/cm³，液限和塑限分别为27.6%，19.4%，初始含水率为1.36%；饱和冻结处理后，试验含水率为22.28%，重度为0.017 kN/m³，孔隙比为0.73；MTS810三轴试验围压设定为0.3~6.0 MPa（共8组），高温冻结温度为-1.0℃；试验过程分为等向加载阶段和轴向应力加载阶段，围压约束的加载在5 min内完成，轴向应力以1.108 mm/min的速度加载，以试样轴向应变达到20%时作为试样破坏条件；当高温冻土轴向应变 ${\varepsilon _{\text{a}}}$ =3%时，卸载偏差应力到零，然后将轴向应变加载到5.75%，再进行第二次卸载，应变的加卸载增量步取2.75%，本文所建模型需要6个参数和2条试验曲线，根据文献[，]及参数拟合确定参数： $c$ =89.3 kPa， $\varphi$ =23.4℃， $\nu$ =0.35， $\kappa$ =-8.6×10^-7， ${\varepsilon _{{\text{v0}}}}$ = 2.14× 10^-2， ${M_{\text{h}}}$ =1.09，等向固结曲线 ${\varepsilon _{\text{v}}} = {\text{2}}{\text{.23}}$ × ${\text{1}}{{\text{0}}^{ - {\text{6}}}} \cdot$ ${\text{exp(0}}{\text{.689ln}}p{\text{)}}$ ，临界状态曲线 ${\varepsilon _{\text{v}}} = {\text{1}}{\text{.76}}$ × ${\text{1}}{{\text{0}}^{ - {\text{5}}}} \cdot$ ${\text{exp}}$ ${\text{(0}}{\text{.734ln}}p{\text{)}}$ ，

从图 7（a）可以看出，本文模型对1.4，2.0，4.0，6.0 MPa围压条件下剪切行为的预测结果与试验数据吻合度较好。在1.4，2.0 MPa围压条件下，高温冻结粉土偏应力一直增加（硬化），在4.0，6.0 MPa围压条件下，高温冻结粉土偏应力先增加后减小（软化）；在同一围压条件下，初始阶段曲线呈线性变化，试样内部以可恢复的弹性变形为主；随着轴向应变的增加，曲线逐渐过渡到弹塑性阶段，此时试样内部以不可恢复的塑性变形为主。对于应变硬化的情况，应力逐渐增加到最大值，曲线斜率近视趋近于零，可认为该强度即为-1.0℃的高温冻结粉土在该围压条件下的极限强度；对于应变软化的情况，当偏应力超过一定值时，曲线斜率由正转为负值，这是由于高温冻土中形成的不稳定冰晶在一定压力条件下发生破坏，试样内孔隙水压力增加，颗粒之间的有效应力降低，进而发生剪切破坏，强度降低。从图 7（b）可以看出，本文模型对体变行为的预测结果与试验数据吻合度较高。在1.4，2.0 MPa较低围压时，试样因侧向变形受到的约束力较小而表现为稳定的剪缩特征；在4.0，6.0 MPa较高围压时，试样在初始阶段因内部冰晶结构微裂隙及孔隙损伤开裂破坏而表现为特定的剪胀特征，随着轴向应变的增加，试样逐渐表现为稳定的剪缩特征。综合图 7（a），（b）可以得出，-1.0℃条件下的高温冻结粉土在低围压条件下总体表现出剪缩硬化特性，在高围压条件下总体表现出剪胀软化特性，本文所建模型可以有效反映-1.0℃条件下的高温冻土剪缩、剪胀、硬化和软化特性。

图 7 模型预测结果与试验数据对比（-1.0℃）

Figure 7. Comparison between model results with experimental data (-1.0℃)

下载: 全尺寸图片幻灯片

4. 结论

（1）与传统模型中使用的 $e{\text{ - }}\ln p$ 变形特征曲线相比，以 ${\varepsilon _{\text{v}}}{\text{ - }}\ln p$ 曲线对高温冻土试样变形进行描述后，NCL和CSL曲线表达简便，有利于双屈服面函数形式的简化表征及模型参数的方便获取。

（2）由当前屈服面与参考屈服面间的关系重新定义了能够反映高温冻土应力状态的固结参数及潜在强度，得到的双屈服面统一本构模型能够有效反映高温冻土的剪缩、剪胀、硬化和软化特性。

（3）本模型在经典修正剑桥模型基础上引入的新参数较少，可以通过常规三轴试验、等向固结试验、回弹再加载试验等基本土工试验获取，模型应用简便，试验结果验证了所建模型的合理性。

图 1 边坡分析模型

Figure 1. Slope analysis model

下载: 全尺寸图片幻灯片

图 2 位移相容条件

Figure 2. Displacement compatibility conditions

下载: 全尺寸图片幻灯片

图 3 边坡地震位移分析流程图

Figure 3. Flow chart of seismic slope displacement analysis

下载: 全尺寸图片幻灯片

图 4 地震动的加速度时程（水平方向）

Figure 4. Time histories of acceleration of seismic motion (horizontal direction)

下载: 全尺寸图片幻灯片

图 5 水平方向地震动加速度的傅里叶谱曲线（半对数坐标系）

Figure 5. Curves of Fourier spectrum of seismic acceleration in horizontal direction (in semilog coordinate)

下载: 全尺寸图片幻灯片

图 6 不同土条的最大滑体等效加速度

Figure 6. Maximum equivalent accelerations of different slices in sliding mass

下载: 全尺寸图片幻灯片

图 7 不同地震作用下边坡水平位移云图

Figure 7. Contour map of horizontal displacement of slopes under different earthquake ground motions

下载: 全尺寸图片幻灯片

图 8 不同地震作用下边坡水平相对位移

Figure 8. Horizontal relative displacements of slopes under different earthquake ground motions

下载: 全尺寸图片幻灯片

图 9 地震动的加速度时程（竖直方向）

Figure 9. Time histories of acceleration of seismic motion (vertical direction)

下载: 全尺寸图片幻灯片

图 10 不同地震作用下边坡水平相对位移

Figure 10. Horizontal relative displacements of slope under different earthquake ground motions

下载: 全尺寸图片幻灯片

图 11 多点地震动的不同分布

Figure 11. Different distributions of multi-point seismic motions

下载: 全尺寸图片幻灯片

图 12 不同分布多点地震动下边坡永久位移

Figure 12. Permanent displacements of slopes under multi-point seismic motions with different distributions

下载: 全尺寸图片幻灯片

图 13 Lexington土石坝剖面图^[41]

Figure 13. Main section of Lexington Dam^[41]

下载: 全尺寸图片幻灯片

图 14 输入的地震动时程

Figure 14. Time histories of input motion

下载: 全尺寸图片幻灯片

图 15 不同土条的最大滑体等效加速度

Figure 15. Maximum equivalent accelerations of different slices in sliding mass

下载: 全尺寸图片幻灯片

图 16 不同地震动作用下坝坡永久位移

Figure 16. Permanent displacements of Lexington Dam under different seismic excitations

下载: 全尺寸图片幻灯片

表 1 地震动记录主要参数

Table 1 Main parameters of earthquake ground motions

地震动编号(RSN)	地震名称	场地名称	矩震级	场地-断层距离/km	峰值加速度 PGA/g		峰值速度 PGV/(cm·s^-1)		平均主导频率^[37] f_m/Hz
地震动编号(RSN)	地震名称	场地名称	矩震级	场地-断层距离/km	水平向	竖向	水平向	竖向	水平向	竖向
4097	Parkfield-02, CA	Slack Canyon	6.00	2.99	0.21	0.11	25.94	8.72	1.40	1.74
1012	Northridge-01	LA 00	6.69	19.07	0.26	0.18	28.85	6.05	2.29	5.95
3932	Tottori, Japan	OKYH14	6.61	26.51	0.45	0.17	23.22	3.13	6.28	10.99

下载: 导出CSV

表 2 Lexington土石坝主要土体参数^[41]

Table 2 Main soil parameters adopted for Lexington Dam^[41]

材料	重度γ/(kN·m^-3)	黏聚力c/kPa	内摩擦角φ/(°)	剪切波速v_s/(m·s^-1)
心墙	18	12.8	27.5	450
坝壳	18	0	36	549
坝基土层1	20	0	40	762
坝基土层2	20	0	40	914
基岩	22	—	—	1219

下载: 导出CSV

参考文献(42)

[1]	NEWMARK N M. Effects of earthquakes on dams and embankments[J]. Géotechnique, 1965, 15(2): 139-160. doi: 10.1680/geot.1965.15.2.139
[2]	MAKDISI F I, SEED H B. Simplified procedure for estimating dam and embankment earthquake-induced deformations[J]. Journal of the Geotechnical Engineering Division, 1978, 104(7): 849-867. doi: 10.1061/AJGEB6.0000668
[3]	RATHJE E M, BRAY J D. One- and two-dimensional seismic analysis of solid-waste landfills[J]. Canadian Geotechnical Journal, 2001, 38(4): 850-862. doi: 10.1139/t01-009
[4]	PAPADIMITRIOU A G, BOUCKOVALAS G D, ANDRIANOPOULOS K I. Methodology for estimating seismic coefficients for performance-based design of earthdams and tall embankments[J]. Soil Dynamics and Earthquake Engineering, 2014, 56: 57-73. doi: 10.1016/j.soildyn.2013.10.006
[5]	RATHJE E M, BRAY J D. Nonlinear coupled seismic sliding analysis of earth structures[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2000, 126(11): 1002-1014. doi: 10.1061/(ASCE)1090-0241(2000)126:11(1002)
[6]	SONG J, FAN Q Q, FENG T G, et al. A multi-block sliding approach to calculate the permanent seismic displacement of slopes[J]. Engineering Geology, 2019, 255: 48-58. doi: 10.1016/j.enggeo.2019.04.012
[7]	SONG J, WU K L, FENG T G, et al. Coupled analysis of earthquake-induced permanent deformations at shallow and deep failure planes of slopes[J]. Engineering Geology, 2020, 274: 105688. doi: 10.1016/j.enggeo.2020.105688
[8]	SONG J, LU Z X, JI J, et al. A fully nonlinear coupled seismic displacement model for earth slope with multiple slip surfaces[J]. Soil Dynamics and Earthquake Engineering, 2022, 159: 107353. doi: 10.1016/j.soildyn.2022.107353
[9]	SONG J, LU Z X, PAN Y H, et al. Investigation of seismic displacements in bedding rock slopes by an extended Newmark sliding block model[J]. Landslides, 2024, 21(3): 461-477. doi: 10.1007/s10346-023-02170-z
[10]	BANDINI V, BIONDI G, CASCONE E, et al. A GLE-based model for seismic displacement analysis of slopes including strength degradation and geometry rearrangement[J]. Soil Dynamics and Earthquake Engineering, 2015, 71: 128-142. doi: 10.1016/j.soildyn.2015.01.010
[11]	JI J, ZHANG W J, ZHANG F, et al. Reliability analysis on permanent displacement of earth slopes using the simplified Bishop method[J]. Computers and Geotechnics, 2020, 117: 103286. doi: 10.1016/j.compgeo.2019.103286
[12]	陈昌富, 李瑞芳. 基于Spencer法边坡地震永久位移计算斜条分法[J]. 公路工程, 2017, 42(3): 1-5, 25. doi: 10.3969/j.issn.1674-0610.2017.03.001 CHEN Changfu, LI Ruifang. Seismic displacements analysis of slopes using inclined slices and newmark method[J]. Highway Engineering, 2017, 42(3): 1-5, 25. (in Chinese) doi: 10.3969/j.issn.1674-0610.2017.03.001
[13]	王家鑫, 夏元友, 王智德. 考虑滑面应变软化效应的边坡震后位移计算方法[J]. 计算力学学报, 2023, https://link.cnki.net/urlid/21.1373.O3.20231208.1400.032. WANG Jia-xin, XIA Yuan-you, WANG Zhi-de. Seismic displacement calculation method of slope with considering strain-softening effect of slip surface[J]. Chinese Journal of Computational Mechanics, 2023. (in Chinese) https://link.cnki.net/urlid/21.1373.O3.20231208.1400.032
[14]	刘爱娟, 崔玉龙, 刘铁新. 考虑动态临界加速度的地震边坡永久位移预测模型研究[J]. 地震工程学报, 2021, 43(2): 445-452. doi: 10.3969/j.issn.1000-0844.2021.02.445 LIU Aijuan, CUI Yulong, LIU Tiexin. A prediction model for permanent displacement of seismic slopes considering dynamic critical acceleration[J]. China Earthquake Engineering Journal, 2021, 43(2): 445-452. (in Chinese) doi: 10.3969/j.issn.1000-0844.2021.02.445
[15]	宋健, 陆朱汐, 谢华威, 等. 基于极限平衡原理的地震边坡浅层和深层耦合滑移分析[J]. 岩土工程学报, 2023, 45(6): 1141-1150. doi: 10.11779/CJGE20220035 SONG Jian, LU Zhuxi, XIE Huawei, et al. Analysis of coupled shallow and deep sliding of slopes induced by earthquake based on limit equilibrium method[J]. Chinese Journal of Geotechnical Engineering, 2023, 45(6): 1141-1150. (in Chinese) doi: 10.11779/CJGE20220035
[16]	高玉峰, 范昭平. 多点、多向地震作用下非圆弧滑面边坡稳定分析通用条分法[J]. 岩土力学, 2010, 31(12): 3816-3822. doi: 10.3969/j.issn.1000-7598.2010.12.020 GAO Yufeng, FAN Zhaoping. Generalized slice method for stability analysis of slope with non-circular slide face under multipoint and multidirection seismic ground motions[J]. Rock and Soil Mechanics, 2010, 31(12): 3816-3822. (in Chinese) doi: 10.3969/j.issn.1000-7598.2010.12.020
[17]	高玉峰. 河谷场地地震波传播解析模型及放大效应[J]. 岩土工程学报, 2019, 41(1): 1-25. doi: 10.11779/CJGE201901001 GAO Yufeng. Analytical models and amplification effects of seismic wave propagation in canyon sites[J]. Chinese Journal of Geotechnical Engineering, 2019, 41(1): 1-25. (in Chinese) doi: 10.11779/CJGE201901001
[18]	蒋乐英, 廖意辉, 王志明, 等. 河谷差异地震作用下拱式倒虹吸地震响应分析[J]. 长江科学院院报, 2022, 39(12): 117-121, 140. JIANG Leying, LIAO Yihui, WANG Zhiming, et al. Seismic response of arch inverted siphon under the excitation of river valley differential earthquake[J]. Journal of Yangtze River Scientific Research Institute, 2022, 39(12): 117-121, 140. (in Chinese)
[19]	WANG G Q, ZHOU X Y, ZHANG P Z, et al. Characteristics of amplitude and duration for near fault strong ground motion from the 1999 Chi-Chi, Taiwan Earthquake[J]. Soil Dynamics and Earthquake Engineering, 2002, 22(1): 73-96. doi: 10.1016/S0267-7261(01)00047-1
[20]	WANG H Y, XIE L L, WANG S Y, et al. Site response in the Qionghai Basin in the Wenchuan earthquake[J]. Earthquake Engineering and Engineering Vibration, 2013, 12(2): 195-199. doi: 10.1007/s11803-013-0162-4
[21]	BERESNEV I A, WEN K L; Nonlinear soil response: a reality?[J]. Bulletin of the Seismological Society of America 1996, 86(6): 1964-1978.
[22]	袁吉. 基于日本KiK-net强震动数据的场地效应研究[D]. 哈尔滨: 中国地震局工程力学研究所, 2021. YUAN Ji. Study on site effect based on Japanese KiK-net strong earthquake data[D]. Harbin: Institute of Engineering Mechanics, China Earthquake Administration, 2021. (in Chinese)
[23]	DING Y, WANG G X, YANG F J. Parametric investigation on the effect of near-surface soil properties on the topographic amplification of ground motions[J]. Engineering Geology, 2020, 273: 105687. doi: 10.1016/j.enggeo.2020.105687
[24]	HE J X, FU H Y, ZHANG Y B, et al. The effect of surficial soil on the seismic response characteristics and failure pattern of step-like slopes[J]. Soil Dynamics and Earthquake Engineering, 2022, 161: 107441. doi: 10.1016/j.soildyn.2022.107441
[25]	王志颖, 郭明珠, 曾金艳, 等. 地震作用下含软弱夹层顺层岩质斜坡动力响应的试验研究[J]. 岩土力学, 2023, 44(9): 2566-2578, 2592. WANG Zhiying, GUO Mingzhu, ZENG Jinyan, et al. Experimental study on dynamic response of bedding rock slope with weak interlayer under earthquake[J]. Rock and Soil Mechanics, 2023, 44(9): 2566-2578, 2592. (in Chinese)
[26]	JIBSON R W. Methods for assessing the stability of slopes during earthquakes: a retrospective[J]. Engineering Geology, 2011, 122(1/2): 43-50.
[27]	栾茂田, 李湛, 范庆来. 土石坝拟静力抗震稳定性分析与坝坡地震滑移量估算[J]. 岩土力学, 2007, 28(2): 224-230, 236. doi: 10.3969/j.issn.1000-7598.2007.02.003 LUAN Maotian, LI Zhan, FAN Qinglai. Analysis and evaluation of pseudo-static aseismic stability and seism-induced sliding movement of earth-rock dams[J]. Rock and Soil Mechanics, 2007, 28(2): 224-230, 236. (in Chinese) doi: 10.3969/j.issn.1000-7598.2007.02.003
[28]	李红军, 迟世春, 钟红, 等. 考虑时程竖向加速度的Newmark滑块位移法[J]. 岩土力学, 2007, 28(11): 2385-2390. doi: 10.3969/j.issn.1000-7598.2007.11.026 LI Hongjun, CHI Shichun, ZHONG Hong, et al. Effects of vertical acceleration time-histories on Newmark sliding block analyses[J]. Rock and Soil Mechanics, 2007, 28(11): 2385-2390. (in Chinese) doi: 10.3969/j.issn.1000-7598.2007.11.026
[29]	STAMATOPOULOS C A, DI B F. Simplified multi-block constitutive model predicting earthquake-induced landslide triggering and displacement along slip surfaces of saturated sand[J]. Soil Dynamics and Earthquake Engineering, 2014, 67: 16-29. doi: 10.1016/j.soildyn.2014.08.008
[30]	SARMA S K, CHLIMINTZAS G O. Co-seismic & post-seismic displacements of slopes[C]//15th ICSMGE TC4 Satellite Conference on" Lessons Learned from Recent Strong Earthquakes. Adapazari, 2001.
[31]	SPENCER E. A method of analysis of the stability of embankments assuming parallel inter-slice forces[J]. Géotechnique, 1967, 17(1): 11-26. doi: 10.1680/geot.1967.17.1.11
[32]	BISHOP A W. The use of the slip circle in the stability analysis of slopes[J]. Géotechnique, 1955, 5(1): 7-17. doi: 10.1680/geot.1955.5.1.7
[33]	KIM J, SITAR N. Direct estimation of yield acceleration in slope stability analyses[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2004, 130(1): 111-115. doi: 10.1061/(ASCE)1090-0241(2004)130:1(111)
[34]	KOTTKE A R, RATHJE E M. Technical manual for Strata. PEER Report 2008/10, Pacific Earthquake Engineering Research Center, College of Engineering[R]. Berkeley: University of California, Berkeley, 2009.
[35]	SEED H B, MARTIN G R. The seismic coefficient in earth dam design[J]. Journal of the Soil Mechanics and Foundations Division, 1966, 92(3): 25-58. doi: 10.1061/JSFEAQ.0000871
[36]	VUCETIC M, DOBRY R. Effect of soil plasticity on cyclic response[J]. Journal of Geotechnical Engineering, 1991, 117(1): 89-107. doi: 10.1061/(ASCE)0733-9410(1991)117:1(89)
[37]	RATHJE E M, ABRAHAMSON N A, BRAY J D. Simplified frequency content estimates of earthquake ground motions[J]. Journal of Geotechnical and Geoenvironmental Engineering, 1998, 124(2): 150-159. doi: 10.1061/(ASCE)1090-0241(1998)124:2(150)
[38]	ITASCA CONSULTING GROUP. Fast Lagrangian Analysis of Continua in 3 Dimensions[M]. MN, USA: Itasca Consulting Group, Minneapolis, 2018.
[39]	DU W Q. Effects of directionality and vertical component of ground motions on seismic slope displacements in Newmark sliding-block analysis[J]. Engineering Geology, 2018, 239: 13-21. doi: 10.1016/j.enggeo.2018.03.012
[40]	HARDER JR L F, BRAY J D, VOLPE R L, et al. Performance of earth dams during the Loma Prieta earthquake proceedings[C]//Second International Conference on Recent Advances in Geotechnical Earthquake Engineering and Soil Dynamics. St Louis, 1991.
[41]	TROPEANO G, CHIARADONNA A, D'ONOFRIO A, et al. An innovative computer code for 1D seismic response analysis including shear strength of soils[J]. Géotechnique, 2016, 66(2): 95-105. doi: 10.1680/jgeot.SIP.15.P.017
[42]	SEED H B, WONG R T, IDRISS I M, et al. Moduli and damping factors for dynamic analyses of cohesionless soils[J]. Journal of Geotechnical Engineering, 1986, 112(11): 1016-1032. doi: 10.1061/(ASCE)0733-9410(1986)112:11(1016)