Seismic response of free-field earthquakes in unsaturated soils with P-wave incidence under thermal effects
-
摘要: 基于单相热弹性介质和非饱和多孔热弹性介质中波的传播理论,建立了平面P波入射下非饱和土自由场地模型,采用亥姆霍兹矢量分解原理,对热效应作用下非饱和土自由场地中的波场进行分析,获得了热效应作用下平面P波入射非饱和土自由场地地震响应的解析解答。通过数值计算,分析了热传导系数、介质温度、热膨胀系数等热物性参数对非饱和土自由场地地震响应所产生的影响规律。研究结果表明:考虑热效应和不考虑热效应两种理论模型下所得到的地表位移放大系数有着明显差异;水平位移放大系数随着热传导系数的增大而逐渐增大,但竖向位移放大系数几乎不发生变化;随着热膨胀系数的增大,水平位移放大系数逐渐增大,竖向位移放大系数逐渐减小;随着介质温度的增大,水平位移放大系数逐渐减小,竖向位移放大系数逐渐增大;此外,随着饱和度的增大,水平位移放大系数逐渐减小,竖向位移放大系数逐渐增大,说明土体中气相的存在对地表位移放大系数有着显著的影响。
-
关键词:
- 非饱和多孔热弹性介质 /
- 自由场地 /
- 地震响应 /
- 热物性参数 /
- 平面P波
Abstract: Based on the theories of wave propagation in single-phase and unsaturated porous thermoelastic media, a free-field model for unsaturated soils under plane P-wave incidence is established, and the Helmholtz vector decomposition principle is employed to analyze the wave fields in the free field of unsaturated soils under the action of thermal effects, and analytical solutions for the seismic response of the free field of unsaturated soils under the action of thermal effects with plane P-wave incidence are obtained. The influences of thermal physical parameters such as thermal conductivity, medium temperature, and thermal expansion coefficient on the seismic response of a free field of unsaturated soils is analyzed through numerical calculations. The results show that the displacement magnification factors of the ground obtained under the two theoretical models considering and not considering thermal effects have apparent differences. The horizontal displacement magnification factor gradually increases with the increase of the thermal conductivity, but the vertical one is almost unchanged. The horizontal displacement magnification factor gradually increases with the increase of the coefficient of thermal expansion, while the vertical one decreases gradually. With the increase of the medium temperature, the horizontal displacement magnification factor decreases, and the vertical one gradually increases. In addition, with the increase of the saturation degree, the horizontal displacement magnification factor gradually decreases, and the vertical one gradually increases, which indicates that the existence of gas phase in the soils has a significant influence on the surface displacement magnification factor. -
0. 引言
边坡抗震问题是岩土工程界的经典问题。地震频发引起的边坡失稳通常伴随明显的三维特征[1]。目前,关于边坡稳定性的研究方法主要有极限平衡法、数值分析法、极限分析法。极限平衡法简单易懂,可以扩展二维方法分析三维边坡。然而,在这些方法中大都包含与二维方法相同的应力和破坏面的分布假设,不满足实际的三维问题,而且忽略了本构关系,得到的结果不是严格的上限解和下限解。数值模拟法不需要假设坍塌机制,通过地震加速度时程研究地震效应,但该方法十分复杂,并且受所用数值方法中规定的人工边界条件和网格尺寸影响[2]。极限分析采用流动法则考虑了岩土材料的本构关系,与极限平衡法相比,理论更严谨,分析过程更简单,已成为解决边坡稳定性问题的一种有效途径[3]。
边坡稳定性分析通常假设土完全干燥或饱和,而实际上边坡处于非饱和状态。饱和土与非饱和土的物理和力学性质存在显著差异,导致边坡稳定性条件不同。因此,需要建立一个更通用的非饱和土屈服准则,以便在非饱和土边坡稳定性分析中使用。Lu等[4]在估算稳定非饱和流条件下土的吸应力方面取得了成功。在该框架下,许多学者对稳定非饱和渗水作用下边坡稳定性进行了研究[5-7]。
此外,表征地震加速度的方法也至关重要。工程常用的方法是拟静力法[7-8],拟静力法(PSM)将地震加速度视为稳定均匀分布的惯性力。然而这种方法没有考虑地震加速度的时空效应,将边坡中任何位置的地震加速度视为常数,这与真实的地震波不符。为了克服以上缺点,Steedman等[9]提出了一种传统的拟动力法(CPDM),该方法用正弦函数研究土中地震加速度的时空变化,可以较好反映地震波的动力特性。虽然传统拟动力法在边坡稳定性分析中得到了一定的应用,但仍存在局限性,如:(a)不满足零应力边界条件;(b)没有考虑土体的阻尼特性。Bellezzal等[10]提出了一种考虑土体阻尼特性的修正拟动力法(MPDM),该方法克服了传统拟动力法的局限性。
综上所述,目前已有的非饱和土边坡抗震稳定性分析中,并未全面分析土的阻尼特性和共振现象。为此,将吸应力方程纳入极限分析运动学中。采用Michalowski等[1]提出的三维旋转破坏模型,提出一种新的三维水平切片法,用于计算非饱和土重度、地震加速度外部功率和内能耗散。通过优化程序得到边坡安全系数的最小上限解,验证了该方法的有效性。最后,通过参数研究,详细讨论了归一化频率、地震系数、三维效应和吸力对边坡三维稳定性的影响。
1. 理论模型
1.1 非饱和土吸应力与有效重度分布
在一维稳定入渗条件下,结合达西定律和Gardner渗透率模型,以及地下水位处零吸力的边界条件,可以得到吸力的解析解[5]:
(ua−uw)=−1αln[(1+q/ks)e−αγwz−q/ks]。 (1) 式中:ua为孔隙气压力;uw为孔隙水压力;α为进气值的倒数;q/ks为垂直比流量;γw为水的重度;z为边坡点到地下水位的距离。
Lu等[11]建立了基质吸力表示的吸应力的闭合形式方程,表示为
σs=−(ua−uw){1+[α(ua−uw)]n}(1−1/1nn) (ua−uw)>0。 (2) 式中:n为描述土壤孔径分布的参数,对于大多数土壤,n的取值范围为1.1~8.5[5]。将式(1)代入式(2)中,可得到在稳定非饱和渗流条件下,吸应力的闭合形式函数[12]:
σs=1αln[(1+q/ks)e−αγwz−q/ks]{1+{−ln[(1+q/ks)e−αγwz−q/ks]}n}(1−1/1nn)。 (3) 非饱和土边坡中的含水率随深度变化,不仅影响着吸力的大小,还会改变土的有效重度,从而影响边坡的稳定性[13]。根据非饱和土的三相比例关系,非饱和土的有效重度可表示为[12]
γ′=(γsat−γw)Gs+[Se−(1−Se)Sr](Gsγw−γsat)Gs−1。 (4) 式中:γsat为土的饱和重度;Gs为土的相对质量密度;Se为有效饱和度;Sr为残余饱和度。Gs的取值范围[7]为2.60~2.80。由于缺乏现场数据,本研究中假设γsat和Gs的值分别为20 kN/m3和2.70。
利用Van-Genuchten的土水特征曲线,有效饱和度与基质吸力之间的关系可以写成[12]
Se={11+[α(ua−uw)]n}1−1/1nn。 (5) 1.2 修正拟动力法
Bellezzal等[10]提出了修正的拟动力法。将土体视为Kelvin-Voigt黏弹性介质,根据Kramer模型,运动方程可表示为
σij=2Gεij+2η∂εij∂t。 (6) 式中:G为土的剪切模量;t为时间;σij和εij分别为应力和应变;η为黏度,η=2Gξ/ωs,ξ为土壤阻尼比,ωs为地震角速度。研究表明:当kv⩽时,可忽略竖向地震力的影响[14]。因此,本文仅考虑水平地震波,则剪切波的一维运动方程表示为
{\rho _{\text{s}}}\frac{{{\partial ^2}{u_{\text{h}}}}}{{\partial {t^2}}} = G\frac{{{\partial ^2}{u_{\text{h}}}}}{{\partial {{z'}^2}}} + \eta \frac{{{\partial ^3}{u_{\text{h}}}}}{{\partial {{z'}^2}\partial t}} 。 (7) 式中: {\rho _{\text{s}}} 为土体密度; {u_{\text{h}}} 为水平位移。将边界条件:(a)边坡坡顶处剪应力为零;(b)边坡坡脚处位移{u_{{\text{toe}}}} = {u_{{\text{h0}}}}\cos ({\omega _{\text{s}}}t)代入方程(7),可得到任意时刻t距离坡脚z'处的水平位移:
{u_{\text{h}}}(z',t) = \frac{{{u_{{\text{h0}}}}}}{{C_{\text{S}}^{\text{2}} + S_\text{s}^2}}[({C_{\text{s}}}{C_{{\text{sz}}}} + {S_{\text{s}}}{S_{{\text{sz}}}})\cos (2{\text{π }}t{\text{/}}T) + \\ ~~~~~~~~~~~~~~~~~({S_{\text{s}}}{C_{{\text{sz}}}} - {C_s}{S_{{\text{sz}}}})\sin (2{\text{π }}t{\text{/}}T)] 。 (8) 式中:T为水平地震周期;{u_{{\text{h0}}}}为初始水平位移;{C_{\text{s}}},{S_{\text{s}}},{C_{{\text{sz}}}},{S_{{\text{sz}}}}的表达式如下:
{C_{\text{s}}} = \cos ({y_{{\text{s1}}}})\cosh ({y_{{\text{s2}}}}) \text{,} (9) {S_{\text{s}}} = - \sin ({y_{{\text{s1}}}})\sinh ({y_{{\text{s2}}}}) \text{,} (10) {C_{{\text{sz}}}} = \cos \left[ {\frac{{{y_{{\text{s1}}}}(H - z')}}{H}} \right]\cosh \left[ {\frac{{{y_{{\text{s2}}}}(H - z')}}{H}} \right] \text{,} (11) {S_{{\text{sz}}}} = - \sin \left[ {\frac{{{y_{{\text{s1}}}}(H - z')}}{H}} \right]\sinh \left[ {\frac{{{y_{{\text{s2}}}}(H - z')}}{H}} \right] \text{,} (12) {y_{{\text{s1}}}} = \frac{{2{\text{π }}H}}{{T{v_{\text{s}}}}}{\left[ {\frac{{\sqrt {1 + 4{\xi ^2}} + 1}}{{2(1 + 4{\xi ^2})}}} \right]^{0.5}} \text{,} (13) {y_{{\text{s2}}}} = - \frac{{2{\text{π }}H}}{{T{v_{\text{s}}}}}{\left[ {\frac{{\sqrt {1 + 4{\xi ^2}} - 1}}{{2(1 + 4{\xi ^2})}}} \right]^{0.5}} 。 (14) 式中:{y_{s1}},{y_{{\text{s2}}}}为归一化频率H/(T \cdot {v_{\text{s}}})与阻尼比\xi 的函数;{v_{\text{s}}}为水平地震波波速。
通过对时间t进行两次微分,可得到滑动面上任一点距离边坡坡脚垂直距离z'处随时间变化的地震加速度表达式:
{a_{\text{h}}}(z',t) = \frac{{{k_{\text{h}}}g}}{{C_{\text{S}}^{\text{2}} + S_{\text{s}}^{\text{2}}}}[({C_{\text{s}}}{C_{{\text{sz}}}} + {S_{\text{s}}}{S_{{\text{sz}}}})\cos (2{\text{π }}t{\text{/}}T) + \\~~~~~~~~~~~~~~~~~{\text{ }}({S_{\text{s}}}{C_{{\text{sz}}}} - {C_{\text{s}}}{S_{{\text{sz}}}})\sin (2{\text{π }}t{\text{/}}T)] 。 (15) 式中: {k_{\text{h}}} 为水平地震系数;g为重力加速度。
2. 三维非饱和边坡地震分析
2.1 三维失效机制
对于摩擦性土,塑性屈服必然伴随着膨胀,这为构造运动许可速度场带来了困难。对于刚性体,构建运动许可的三维机构相对简单,满足相关联流动法则的局部表面必须与顶点为2\varphi '的角相切。Michalowski等[1]构建了一种牛角状模型,如图 1所示,该机构有一个对称面,上、下轮廓由对数螺旋线A'C'和AC定义:
\rho ' = {\rho '_0}{{\text{e}}^{ - \left( {\theta - {\theta _0}} \right)\tan \varphi '}} \text{,} (16) \rho = {\rho _0}{{\text{e}}^{\left( {\theta - {\theta _0}} \right)\tan \varphi '}} 。 (17) 式中: {\rho '_0} 和\rho _0分别表示极径{OA '},OA;{\theta _0}为对数螺旋线的初始夹角; \rho ' 和 \rho 分别表示为旋转中心O到曲线A'C'和AC上一般点的距离。
如图 1所示,破坏机制的截面由不断扩大的圆旋转而成,图中{r_{\text{c}}}为点O到圆截面圆心的距离,R为圆的半径,由以下方程定义:
{r_{\text{c}}} = (\rho + \rho ')/2 \text{,} (18) R = (\rho + \rho ')/2 。 (19) 当边坡的宽高比较小时,土坡会呈现明显的三维特征,如图 2所示。为了允许向平面应变机制过渡,加入了“平面插入块”b,当B \to \infty 时,该破坏机构可近似看作二维破坏。设无插入块时的边坡宽度为B',则平面插入块宽度b可表示为
b = B - B' 。 (20) 2.2 上限定理
由虚功原理可知,假想破坏机构的外力做功不会超过内能耗散功率,即
\int_V \boldsymbol{\sigma}_{i j} \dot{\boldsymbol{\varepsilon}}_{i j} \mathrm{~d} V=\int_A T_i \dot{u}_i \mathrm{~d} A+\int_V F_i \dot{u}_i \mathrm{~d} V。 (21) 式中: \dot{\boldsymbol{\varepsilon}}_{i j} 为运动许可速度场的应变率张量; {\boldsymbol{\sigma} _{ij}} 为应变率张量对应的应力张量;{T_i}和{F_i}分别为面力和体力; {\dot u_i} 为运动许可速度场;A和V分别为积分对象的面积和体积。
2.3 能量平衡方程
外部功率W包括土重做功{W_{\gamma '}}和地震功{W_{\text{e}}},内能耗散包括土的有效黏聚力内能耗散{D_{c'}}和毛细黏聚力引起的内能耗散{D_{c''}}。能量平衡方程表示为
{W_{\gamma '}} + {W_{\text{e}}} = {D_{c'}} + {D_{c''}} 。 (22) 本文提出新的三维水平切片积分法计算非饱和土的外部功率。如图 3(a)所示,将非饱和土体划分为m个离散层,这些层的厚度足够薄,以确保有效重度和地震加速度可以作为每一层内的常数。
简化假设有效考虑了土的有效重度和地震加速度的时空变化特性。通过累加所有土层单元土体重力做功率,可得到非饱和土块自重总功率:
\begin{array}{l} {W_{\gamma '}} = \int_V {\gamma 'v\cos \theta {\text{d}}V} \hfill \\ ~~~~~{\text{ }} = \sum\nolimits_{k = 1}^m {\omega \int_{{\theta _{k'}}}^{{\theta _k}} {\gamma '({z_1})} } \cdot \rho _{k'}^2 \cdot {\sin ^2}{\theta _{k'}} \cdot \hfill \\ \end{array}{} \\ ~~~~~~~\frac{{\cos \theta }}{{{{\sin }^3}\theta }} \cdot \left( {2\sqrt {{R^2} - {{({r_{{s_{k'}}}} - {r_{\text{c}}})}^2}} + b} \right) \cdot {\text{d}}h{\text{d}}\theta 。 (23) 式中: \omega 为角速度;{\text{d}}h为无穷小体积元素的高度;\gamma '为某一离散土层的有效重度,通过将{z_1}代入式(4)可以很容易得到;{\rho _{k'}},{\theta _{k'}}分别为微元体(图 3阴影部分)至旋转轴的极径和极角;{z_0}为边坡坡脚至地下水位地垂直距离;变量{r_{{s_{k'}}}},{\theta _k},{\theta _{k'}},{z_1}都可以很容易地从图 3(b)中的几何关系得到,表达式如下:
{z_1} = {z_0} + \frac{k}{m}\left[ {{\rho _{\text{h}}} \cdot \sin {\theta _{\text{h}}} - {\rho _0} \cdot \sin {\theta _0}} \right]\text{,} (24) {\rho _0}{{\text{e}}^{\left( {{\theta _{k'}} - {\theta _0}} \right)\tan \varphi '}} \cdot \sin {\theta _{k'}} - (1 - k/m) \cdot H = {\rho _0}\sin {\theta _0}\text{,} (25) {z_{k'}} = {\rho _{\text{h}}}\sin {\theta _{\text{h}}} - {\rho _{k'}}\sin {\theta _{k'}} \text{,} (26) {\theta _k} = \arctan \frac{{{\rho _{k'}} \cdot \sin {\theta _{k'}} \cdot \sin \beta }}{{{z_{k'}}\cos \beta + {\rho _{\text{h}}}\cos {\theta _{\text{h}}}\sin \beta }}\text{,} (27) {r_{{s_{k'}}}} = \frac{{{\rho _0}\sin {\theta _0} + (1 - k/m) \cdot H}}{{\sin \theta }},\theta \in \left[ {{\theta _{k'}},{\theta _k}} \right]。 (28) 同样地,地震功率表示为
\begin{array}{l} {W_{\text{e}}} = \int_V {\frac{{\gamma '}}{g}} \cdot {a_{\text{h}}} \cdot v\sin \theta {\text{d}}V \hfill \\ {\text{ }}~~~~~ = \sum\nolimits_{k = 1}^m {\frac{\omega }{g}} \int_{{\theta _{k'}}}^{{\theta _k}} {\gamma '({z_1}) \cdot {a_{\text{h}}}({z_2})} \cdot \rho _{k'}^2 \cdot {\sin ^2}{\theta _{k'}} \cdot \hfill \\ \end{array}{} \\~~~~~~~~~~~~~\frac{1}{{{{\sin }^2}\theta }} \cdot \left( {2\sqrt {{R^2} - {{({r_{{s_{k'}}}} - {r_{\text{c}}})}^2}} + b} \right) \cdot {\text{d}}h{\text{d}}\theta 。 (29) 式中:{a_{\text{h}}}为无穷小体积元素的地震加速度,将{z_2}代入方程(15)可得到,{z_2}表达式如下:
{z_2} = \frac{k}{m}\left[ {{\rho _{\text{h}}} \cdot \sin {\theta _{\text{h}}} - {\rho _0} \cdot \sin {\theta _0}} \right]。 (30) 总内能耗散包括有效黏聚力和毛细黏聚力部分:
D = D_{c'}^{{\text{3D}}} + D_{c''}^{{\text{3D}}} + D_{c'}^{{\text{2D}}} + D_{c''}^{{\text{2D}}} 。 (31) 式中:有效黏聚力耗散部分的表达式可以在文献[1]中找到,对于非饱和土块,每一层的毛细管内黏聚力可以写成
{c''_k} = - {\sigma _{\text{s}}}(z){|_{z = {z_1}}} \cdot \tan \varphi ' 。 (32) 式中: {\sigma _{\text{s}}}z{|_{z = z1}} 表示第k层土的吸应力,通过将z = z1带入方程(3)得到。
对于黏性摩擦土边坡(c′>0,φ′>0),每层土的毛细黏聚力产生的总耗散功率(体积内耗散{D_{\text{V}}}和速度不连续面上耗散{D_{\text{t}}})基于高斯散度定理可以推导为 D_{c''}^{{\text{3D}}} = \sum\nolimits_{k = 1}^m {{{c''}_k}} \cot \varphi ' \cdot
\left\{ {\int_{{s_{k'k}}} {{v_i}{n_i}{\text{d}}S + \int_{{s_{kk - 1}}} {{v_i}{n_i}{\text{d}}S + } \int_{{s_{k - 1k' - 1}}} {{v_i}{n_i}{\text{d}}S} } } \right\} 。 (33) 式中:{S_{k'k}}为土层k的上表面;{S_{k - 1k' - 1}}为土层k的下表面;{S_{kk - 1}}为土层k的坡面;非饱和土块ABC的毛细黏聚力引起的内能耗散可以分别写成
D_{c''}^{{\text{2D}}} = - b\omega \rho _0^2\tan \varphi '\int_{{\theta _0}}^{{\theta _{\text{h}}}} {{{\left. {{\sigma _s}(z)} \right|}_{z = z3}}{{\text{e}}^{2\left( {\theta - {\theta _0}} \right)\tan \varphi '}}{\text{d}}\theta } \text{,} (34) \begin{array}{l} D_{c''}^{{\text{3D}}} = 2\omega \rho _0^2\left\{ {{{\sin }^2}{\theta _0}\int_{{\theta _0}}^{{\theta _B}} {{\sigma _{\text{s}}}} (z){|_{z = z4}}\frac{{\cos \theta }}{{{{\sin }^3}\theta }}} \right. \cdot \hfill \\ ~~~~~{\text{ }}\sqrt {{R^2} - {{({r_{{s_1}}} - {r_{\text{c}}})}^2}} {\text{d}}\theta + {{\text{e}}^{2\left( {{\theta _{\text{h}}} - {\theta _0}} \right)\tan \varphi '}} \cdot \hfill \\ ~~~~~{\text{ }}{\sin ^2}({\theta _h} + \beta )\int_{{\theta _B}}^{{\theta _{\text{h}}}} {{\sigma _s}} (z){|_{z = {z5}}}\frac{{\cos (\theta + \beta )}}{{{{\sin }^3}(\theta + \beta )}} \cdot \hfill \\ ~~~~~{\text{ }}\left. {\sqrt {{R^2} - {{({r_{{{\text{s}}_{\text{2}}}}} - {r_{\text{c}}})}^2}} {\text{d}}\theta } \right\} - 2\omega \cot \varphi ' \cdot \hfill \\ ~~~~~{\text{ }}\sum\nolimits_{k = 1}^{m - 1} {\rho _{k'}^2} {\sin ^2}{\theta _{k'}}({{c''}_k} - {{c''}_{k + 1}}) \cdot \hfill \\ \end{array}{} \int_{{\theta _{k'}}}^{{\theta _k}} {\frac{{\cos \theta }}{{{{\sin }^3}\theta }}} \sqrt {{R^2} - {{({r_{{s_{k'}}}} - {r_{\text{c}}})}^2}} {\text{d}}\theta 。 (35) 式中: {\sigma _{\text{s}}}(z){|_{z = {z_3}}} , {\sigma _{\text{s}}}(z){|_{z = {z_4}}} , {\sigma _{\text{s}}}(z){|_{z = {z_5}}} 分别描述了插入件、坡顶、坡面特定点处的吸应力,通过将变量{z_3},{z_4},{z_5}分别代入方程(3)中得到。变量{r_{{\text{s1}}}},{r_{{\text{s2}}}},{\theta _B},{z_3},{z_4},{z_5}的具体表达式可表示为
{r_{{\text{s1}}}} = {\rho _0}\frac{{\sin {\theta _0}}}{{\sin \theta }}\text{,} (36) {r_{s2}} = {\rho _0}{{\text{e}}^{\left( {{\theta _{\text{h}}} - {\theta _0}} \right)\tan \varphi '}}\frac{{\sin ({\theta _{\text{h}}} + \beta )}}{{\sin (\theta + \beta )}}\text{,} (37) {\theta _B} = \arctan \frac{{\sin {\theta _0}}}{{\cos {\theta _0} - A}} \text{,} (38) \begin{aligned} A & =\frac{\sin \left(\theta_{\mathrm{h}}-\theta_0\right)}{\sin \theta_{\mathrm{h}}}- \\ & \frac{\mathrm{e}^{\left(\theta_{\mathrm{h}}-\theta_0\right) \tan \varphi^{\prime}} \sin \theta_{\mathrm{h}}-\sin \theta_0}{\sin \theta_{\mathrm{h}} \sin \beta} \sin \left(\theta_{\mathrm{h}}+\beta\right), \end{aligned} (39) {z_3} = {z_0} + {\rho _h}\sin {\theta _{\text{h}}} - {\rho _\theta }\sin \theta \text{,} (40) {z_4} = {z_0} + {\rho _{\text{h}}}\sin {\theta _{\text{h}}} - {\rho _0}\sin {\theta _0} \text{,} (41) {z_5} = {z_0} + {\rho _{\text{h}}}\sin {\theta _{\text{h}}} - {r_{{{\text{s}}_{\text{2}}}}}\sin \theta 。 (42) 2.4 安全系数和优化程序
安全系数计算方法包括强度折减法(SRM)和重度增加法(GIM)。Yang等[15]讨论了这两种方法的特点,发现两种安全系数{F_{\text{s}}}都是有效的。SRM仍是边坡安全评估中最受欢迎的方法,但它只能给出一个关于{F_{\text{s}}}的隐式表达式,且在三维条件下求解耗时。相反,GIM可以获得{F_{\text{s}}}的显式表达式,定义简单明确。GIM将安全系数{F_{\text{s}}}定义为内部能量耗散功率D与实际外部功率W的比值,其中内能耗散功率D包括有效黏聚力耗散{t \mathord{\left/ {\vphantom {t T}} \right. } T}和毛细黏聚力耗散{b \mathord{\left/ {\vphantom {b H}} \right. } H},外力做功部分包括重力所做功率{W_{\gamma '}}和地震力功率{W_{\text{e}}},即
{F_{\text{s}}} = \frac{{D_{c'}^{{\text{3D}}} + D_{c''}^{{\text{3D}}} + D_{c'}^{{\text{2D}}} + D_{c''}^{{\text{2D}}}}}{{{W_{\gamma '}} + {W_{\text{e}}}}} 。 (43) 边坡安全系数 {F_{\text{s}}} 是关于变量{\theta _0},{\theta _{\text{h}}}, {\rho '_0}/{\rho _0} ,{t \mathord{\left/ {\vphantom {t T}} \right. } T},{b \mathord{\left/ {\vphantom {b H}} \right. } H}的目标函数,利用Matlab软件编写代码进行计算。根据三维失效机制,这些变量应满足约束:
\left.\begin{array}{l} 0<\theta_0<\theta_B<\theta_{\mathrm{h}}<\pi, \\ 0<\rho_0^{\prime} / \rho_0<1, \\ 0<t / T \leqslant 1, \\ 0 \leqslant b / H<B / H 。\quad \end{array}\right\} (44) 求边坡安全系数的最小上限解可以转化为非线性约束条件下的最小值优化问题。为了防止陷入局部最优解,本文采用随机搜索法[16]。同时为了考虑地震波的时空变化特性,将地震时间t纳入到优化算法中,以寻找边坡安全系数的最小上限解。
2.5 对比与分析
为了验证本文开发的优化程序的有效性,将本章计算结果与已出版文献进行对比,比较中采用c' = 20{\text{ kPa}},{\gamma _{{\text{sat}}}} = 20{\text{ kN/}}{{\text{m}}^3}的典型值。边坡高度可以很容易地从Gao等[17]提供的临界高度({{\gamma H} \mathord{\left/ {\vphantom {{\gamma H} c}} \right. } c})中得到。边坡几何参数对本文提出的水平切片法的精度有一定影响,表 1出了无吸力({q \mathord{\left/ {\vphantom {q {{k_{\text{s}}} = - 1}}} \right. } {{k_{\text{s}}} = - 1}})情况下在选定参数\beta ,\varphi ,{B \mathord{\left/ {\vphantom {B H}} \right. } H}不同层数得到的安全系数。可以发现,安全系数非常接近理论值({F_{\text{s}}} = 1.0),说明本文提出的优化程序是合理有效的。并且当层数达到200时,分层数的增加对水平切片法误差的影响几乎可以忽略不计。因此,本文综合考虑了误差影响及搜索效率的影响,将层数m设置为200。
案例 β/(°) φ/(°) B/H γH/c(a) F_{\text{s}}^{{\text{(b)}}} 分层数m(c) 200 300 400 500 1 60 15 1.0 12.831 1.000 1.006 1.006 1.005 1.004 2 75 15 0.8 11.074 1.000 1.003 1.003 1.003 1.002 3 75 30 3.0 11.120 1.000 1.004 1.002 1.001 1.001 4 90 30 1.0 10.503 1.000 0.998 0.998 0.998 0.997 5 90 30 1.5 8.704 1.000 1.005 1.004 1.003 1.003 注:(a):文献[17]中提供的临界高度;(b):与临界高度值相对应的理论安全系数;(c):本文计算的不同层数下无吸力的安全系数。 3. 参数研究
本节对非饱和土边坡稳定性进行了一系列参数研究。详细讨论了地震波频率、地震系数、三维效应以及吸力对边坡稳定性的影响。由于缺乏相关模型,渗透和地震力对边坡的影响是不耦合的。基本参数设置为:H = 5{\text{ m}},{\text{ z}_0} = 0,\beta = {60^ \circ }。
3.1 地震波动力参数研究
本节详细讨论了修正拟动力法(MPDM)使用参数(剪切波振动周期T和剪切波速{v_{\text{s}}}以及阻尼比\xi )对边坡稳定性的影响。基本参数设置为:T = 0.2{\text{ s}},{k_{\text{h}}} = 0.2,{B \mathord{\left/ {\vphantom {B {H = 2.0}}} \right. } {H = 2.0}},{q \mathord{\left/ {\vphantom {q {{k_{\text{s}}}}}} \right. } {{k_{\text{s}}}}} = 0。选择表 2中列出的两种沿湿润路径的真实土参数进行分析。
表 2 湿润条件下不同土壤参数Table 2. Different soil parameters under wetting conditions土壤类型 c'/
kPaφ'/
(°)Srw αw/
kPa-1nw k_{{\text{sat}}}^{\text{w}}/
(10-7m·s-1)A 8.8 24 0.303 0.13 1.52 0.039 B 6.6 37 0.413 0.17 2.20 0.146 注:数据来自于文献[18]。 已有研究表明:地震加速度的大小受地震波长影响(即波速{v_{\text{s}}}和周期T的乘积),边坡的高度与波长之间的关系就直接决定了加速度的大小[19]。因此,在本节中将对归一化频率{H \mathord{\left/ {\vphantom {H {{\lambda _{\text{s}}}}}} \right. } {{\lambda _{\text{s}}}}}({\lambda _{\text{s}}} = T \cdot {v_{\text{s}}})作为参数研究,并提供了不同放大系数{f_{\text{a}}}下传统拟动力法(CPDM)求得的安全系数的比较。
图 4给出了修正拟动力法和传统拟动力法两种方法的比较。可以看出,传统拟动力法计算得到的安全系数曲线是一条随归一化频率增大而呈非线性增大的曲线,这是因为该方法没有考虑场地效应和阻尼,而是依赖于单一因素{f_{\text{a}}}来考虑土壤放大系数。而对于修正拟动力法,安全系数曲线的变化是非单调的,这是因为当地震波频率接近土体的固有频率时(式(45))就会发生共振现象[19]。
H/(T\cdot {v}_{\text{s}})=\frac{n}{2}-\frac{1}{4}\text{ }。 (45) 式中: n=1,\text{ }2,\text{ }3,\text{ }\cdots ,用于计算边坡共振时土的固有频率。
当发生共振时,土壤中的地震驱动力与边坡滑移方向一致,地震功率全为正功。并且地震力施加在边坡中的能量也会累积,由于土壤的放大效应,振动沿着边坡高度不断加强,地震加速度不断升高,使得地震力功率增大,边坡安全系数迅速降低。同时,安全系数曲线的振幅会随着阻尼比\xi 的减小而迅速增大。这是因为土体阻尼比越小,地震波克服阻尼损失的能量越小,地震引起的外功越大,最终{F_{\text{s}}}减小。对于与土体固有频率相差较大的{H \mathord{\left/ {\vphantom {H {{\lambda _{\text{s}}}}}} \right. } {{\lambda _{\text{s}}}}},阻尼比\xi 对{F_{\text{s}}}的影响很小,为了更好地说明这种现象。图 5中绘出了不同阻尼比下边坡坡顶处归一化加速度振幅与归一化频率{H \mathord{\left/ {\vphantom {H {{\lambda _{\text{s}}}}}} \right. } {{\lambda _{\text{s}}}}}的关系图,可以观察到,当振动频率接近共振频率时,加速度振幅会显著快速增加,最大振幅出现在{H \mathord{\left/ {\vphantom {H {{\lambda _{\text{s}}}}}} \right. } {{\lambda _{\text{s}}}}} = {1 \mathord{\left/ {\vphantom {1 4}} \right. } 4}处,并且随着阻尼比的减小而显著增大。
3.2 地震系数研究
为了体现修正拟动力法与传统拟动力法和拟静力法的差异性,图 6给出了3种方法的安全系数随地震系数和归一化频率变化趋势图。基本参数设置为:{B \mathord{\left/ {\vphantom {B {H = 2.0}}} \right. } {H = 2.0}},{q \mathord{\left/ {\vphantom {q {{k_{\text{s}}}}}} \right. } {{k_{\text{s}}}}} = 0,\xi = 0.1,T = 0.2{\text{ s}}。
当{H \mathord{\left/ {\vphantom {H {{\lambda _{\text{s}}}}}} \right. } {{\lambda _{\text{s}}}}} \to 0时,修正拟动力法可以退化为拟静力法。拟静力法无法体现运动频率、坡顶与坡脚之间波形的相位差以及土的放大效应,因此,在图 6中该方法仅为一点。而对于传统拟动力法,放大因子{f_{\text{a}}} = 1,这里没有考虑土体放大效应。传统拟动力法无法预测边坡安全系数的局部最小值,当地震波频率接近土体固有频率时会明显高估边坡的稳定性[20-21]。当归一化频率较很小时,3种方法得到的安全系数几乎相同,且均随着水平地震系数的增加而显著降低。说明水平地震加速度系数对边坡地震稳定性影响较大。
3.3 三维效应的影响
在本节中,基于修正拟动力法研究了不同{k_{\text{h}}}下边坡三维效应的影响。基本参数设置为:\xi = 0.1,{q \mathord{\left/ {\vphantom {q {{k_{\text{s}}} = }}} \right. } {{k_{\text{s}}} = }} 0,T = 0.2{\text{ s}},{H \mathord{\left/ {\vphantom {H {{\lambda _{\text{s}}} = {1 \mathord{\left/ {\vphantom {1 4}} \right. } 4}}}} \right. } {{\lambda _{\text{s}}} = {1 \mathord{\left/ {\vphantom {1 4}} \right. } 4}}}。
从图 7中可以看出,在同一{B \mathord{\left/ {\vphantom {B H}} \right. } H}下,边坡安全系数随{k_{\text{h}}}的增大呈先快后慢的下降趋势。边坡安全系数随宽高比的增加而逐渐减小,并最终趋近于二维效果。当{B \mathord{\left/ {\vphantom {B H}} \right. } H} \lt 3.0时,边坡三维效应明显,安全系数随{k_{\text{h}}}的增大下降较快。当{B \mathord{\left/ {\vphantom {B H}} \right. } H} \gt 3.0时,边坡三维效应减弱,{F_{\text{s}}}的下降趋势逐渐变缓。这说明边坡的三维效应与地震波的动力特性密切相关。在实际工程中,对于宽高比较小的边坡,应根据实际情况考虑其三维效应。
3.4 基质吸力的影响
非饱和土边坡中水分含量的不同,导致了确定临界面失效的差异性。因此,研究不同入渗条件下边坡中不同位置处土壤的非饱和特性是十分必要的。图 8(a),(b)分别显示了在不同归一化入渗条件下,土壤A非饱和土层(高度H=5 m)中毛细黏聚力c''和有效重度\gamma '的分布。可以看出,毛细黏聚力和有效重度的剖面都沿深度变化呈非线性特征。从图 8(a)中可以看出,土壤A的毛细黏聚力在远离地下水位时总是不断增加,随着归一化入渗率({{ - q} \mathord{\left/ {\vphantom {{ - q} {{k_{\text{s}}}}}} \right. } {{k_{\text{s}}}}})的增加,毛细黏聚力不断降低,当{{ - q} \mathord{\left/ {\vphantom {{ - q} {{k_{\text{s}}}}}} \right. } {{k_{\text{s}}}}} = 1时,毛细黏聚力等于0。土壤A的有效重度剖面(图 8(b))也随着渗透速率的增加,逐渐接近饱和条件下的土壤重度。
图 9描述了土壤A和B的非饱和边坡在不同归一化入渗率下安全系数随{k_{\text{h}}}的变化。参数设置为:B/H = 2.0,\xi = 0.1,T = 0.2{\text{ s}},H/{\lambda _{\text{s}}} = 1/4。
可以发现吸力对边坡稳定性有显著的提升作用,对于土A,当{k_{\text{h}}} = 0.1时, - {q \mathord{\left/ {\vphantom {q {{k_{\text{s}}}}}} \right. } {{k_{\text{s}}}}}从0增加到0.2,{F_{\text{s}}}从2.112降到了1.474,降低了43.3%;当 - {q \mathord{\left/ {\vphantom {q {{k_{\text{s}}}}}} \right. } {{k_{\text{s}}}}}从0.4增加到0.6时,{F_{\text{s}}}从1.202减小到0.992,降低了21.1%。随着{k_{\text{h}}}和 - {q \mathord{\left/ {\vphantom {q {{k_{\text{s}}}}}} \right. } {{k_{\text{s}}}}}的增大,吸力对边坡的影响会逐渐减弱。对比两种土得到的安全系数变化趋势图,可以发现土A的安全系数会显著降低,说明对于细粒土边坡来说,吸力诱导效应较强,而对于粗颗粒土(n \gt 2.0)就会相对较弱。
4. 结论
基于极限分析上限定理和一维稳定入渗模型,研究了非饱和土边坡的三维地震稳定性。提出一种改进的三维水平切片法,有效考虑了非饱和土中水分含量和地震加速度沿边坡高度变化的时空变化特性。与已出版文献进行对比,验证了该方法的有效性。采用修正拟动力法考虑了地震作用的影响,并比较了修正拟动力法、传统拟动力法、拟静力法的特性。分析了各种参数对三维边坡稳定性的影响,主要得到以下3点结论。
(1)当边坡受到与土体固有频率接近的地震波作用时,会发生共振现象,这会导致边坡的稳定性迅速降低,并且随着水平地震系数{k_{\text{h}}}的增加和阻尼比\xi 的降低,效果愈加明显。
(2)忽略边坡的三维效应会低估边坡的稳定性。当{B \mathord{\left/ {\vphantom {B H}} \right. } H} \lt 3时,边坡三维效应明显,当{k_{\text{h}}}增加时,安全系数{F_{\text{s}}}会迅速降低。随着{B \mathord{\left/ {\vphantom {B H}} \right. } H}的增加,三维效应会逐渐过渡到平面应变机制,{F_{\text{s}}}的下降趋势也会变缓。这说明,对于三维效应明显的边坡,地震的影响尤为明显。
(3)吸力会提高边坡的安全系数,并且土壤参数的选择强烈影响着吸力诱导效应。对于地震条件下的边坡,随着归一化入渗率 - {q \mathord{\left/ {\vphantom {q {{k_{\text{s}}}}}} \right. } {{k_{\text{s}}}}}的增加,边坡稳定性对地震烈度的敏感性变弱。
-
表 1 单相热弹性介质的物理力学参数
Table 1 Physical-mechanical parameters of single-phase thermoelastic media
参数 λe/Pa μe/Pa ρe/(kg·m-3) βs/K-1 Ke/(J·s-1·m-1·K-1) 量值 12.0×109 8.0×109 2700 4.0×10-4 2.5 参数 Te/K τqe/s τθe/s T0/K cse/(J·kg-1·K-1) 量值 293.2 2.0×10-7 1.5×10-7 300 1000 表 2 非饱和热弹性介质的物理力学参数
Table 2 Physical-mechanical parameters of unsaturated thermoelastic media
参数 λ/Pa μ/Pa ρs/(kg·m-3) βT/K-1 K/(J·s-1·m-1·K-1) 量值 4.4×109 2.8×109 2650 1.0×10-4 0.6 参数 T/K τq/s τθ/s T/K cs/(J·kg-1·K-1) 量值 293.2 2.0×10-7 1.5×10-7 300 1000 参数 cl/(J·kg-1·K-1) cg/(J·kg-1·K-1) βST/K-1 βWT/K-1 βWP/(Pa-1) 量值 4180 1900 7.8×10-6 2.1×10-4 4.58×10-10 参数 Sl Ssat Sres n k/m2 量值 0.6 1.0 0.05 0.4 1.0×10-10 参数 ρl/(kg·m-3) Ks/kPa Pg/kPa m d 量值 1000 3.5×107 101.3 0.5 2 -
[1] BIOT M A. General solutions of the equations of elasticity and consolidation for a porous material[J]. Journal of Applied Mechanics, 1956, 23(1): 91-96. doi: 10.1115/1.4011213
[2] LIU Z X, LIU J Q, MENG S B, et al. Diffraction of elastic waves by a fluid-filled crack in a fluid-saturated poroelastic half-space[J]. Geophysical Journal International, 2021, 225: 1530-1553. doi: 10.1093/gji/ggab043
[3] 范凯祥, 申玉生, 闻毓民, 等. 平面Rayleigh波入射下饱和土中浅埋隧道复合式衬砌的动力响应[J]. 岩土工程学报, 2022, 44(3): 444-455. doi: 10.11779/CJGE202203006 FAN Kaixiang, SHEN Yusheng, WEN Yumin, et al. Dynamic response of composite linings of shallowly buried tunnels in saturated soils subjected to incidence of plane Rayleigh waves[J]. Chinese Journal of Geotechnical Engineering, 2022, 44(3): 444-455. (in Chinese) doi: 10.11779/CJGE202203006
[4] 禹海涛, 王治坤, 刘中宪. SV波入射下均匀饱和地层渗透系数对深埋隧道的影响机制[J]. 岩土工程学报, 2022, 44(2): 201-211. doi: 10.11779/CJGE202202001 YU Haitao, WANG Zhikun, LIU Zhongxian. Influence mechanism of permeability coefficient in homogeneously saturated strata on responses of deep tunnels under incidence of SV waves[J]. Chinese Journal of Geotechnical Engineering, 2022, 44(2): 201-211. (in Chinese) doi: 10.11779/CJGE202202001
[5] 巴振宁, 梁建文, 梅雄一. 斜入射平面SH波在层状饱和半空间中沉积谷地周围的三维散射[J]. 岩土工程学报, 2013, 35(3): 476-486. http://cge.nhri.cn/article/id/15001 BA Zhenning, LIANG Jianwen, MEI Xiongyi. Three-dimensional scattering of obliquely incident plane SH waves in alluvial valley embedded in layered saturated half-space[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(3): 476-486. (in Chinese) http://cge.nhri.cn/article/id/15001
[6] 邹新军, 杨紫健, 吴文兵. 非饱和土地基中端承桩对SH波的水平地震响应[J]. 岩土工程学报, 2024, 46(1): 72-80. doi: 10.11779/CJGE20221125 ZOU Xinjun, YANG Zijian, WU Wenbing. Horizontal seismic response of end-bearing piles in unsaturated soil foundation under SH waves[J]. Chinese Journal of Geotechnical Engineering, 2024, 46(1): 72-80. (in Chinese) doi: 10.11779/CJGE20221125
[7] WEI C F, MURALEETHARAN K K. A continuum theory of porous media saturated by multiple immiscible fluids: Ⅰ linear poroelasticity[J]. International Journal of Engineering Science, 2002, 40(16): 1807-1833. doi: 10.1016/S0020-7225(02)00068-X
[8] CHEN W Y, XIA T D, HU W T. A mixture theory analysis for the surface-wave propagation in an unsaturated porous medium[J]. International Journal of Solids and Structures, 2011, 48(16/17): 2402-2412.
[9] CHEN W Y, XIA T D, CHEN W, et al. Propagation of plane P-waves at interface between elastic solid and unsaturated poroelastic medium[J]. Applied Mathematics and Mechanics, 2012, 33(7): 829-844.
[10] TOMAR S K, ARORA A. Reflection and transmission of elastic waves at an elastic/porous solid saturated by two immiscible fluids[J]. International Journal of Solids and Structures, 2006, 43(7/8): 1991-2013.
[11] LORD H W, SHULMAN Y. A generalized dynamical theory of thermoelasticity[J]. Journal of the Mechanics and Physics of Solids, 1967, 15(5): 299-309.
[12] PECKER C, DERESIEWICZ H. Thermal effects on wave propagation in liquid-filled porous media[J]. Acta Mechanica, 1973, 16(1/2).
[13] YOUSSEF H M. Theory of generalized porothermoelasticity[J]. International Journal of Rock Mechanics and Mining Sciences, 2007, 44(2): 222-227.
[14] 柳鸿博, 周凤玺, 岳国栋, 等. 非饱和土中热弹性波的传播特性分析[J]. 岩土力学, 2020, 41(5): 1613-1624. LIU Hongbo, ZHOU Fengxi, YUE Guodong, et al. Analysis of propagation characteristics of thermoelastic waves in unsaturated soils[J]. Rock and Soil Mechanics, 2020, 41(5): 1613-1624. (in Chinese)
[15] 郑荣跃, 刘干斌, 邓岳保, 等. SV波在饱和多孔热弹性介质平面界面上的反射[J]. 岩土工程学报, 2013, 35(增刊2): 839-843. http://cge.nhri.cn/article/id/15502 ZHENG Rongyue, LIU Qianbin, DENG Yuebao, et al. Reflection of SV waves on plane interface of saturated porous thermoelastic media[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(S2): 839-843. (in Chinese) http://cge.nhri.cn/article/id/15502
[16] 柳鸿博, 周凤玺, 郝磊超. 平面S波在饱和多孔热弹性介质边界的反射问题研究[J]. 地震工程学报, 2021, 43(1): 105-112. LIU Hongbo, ZHOU Fengxi, HAO Leichao. Reflection of plane S waves at the boundary of saturated porous thermo-elastic media[J]. China Earthquake Engineering Journal, 2021, 43(1): 105-112. (in Chinese)
[17] ZHOU F X, ZHANG R L, LIU H B, et al. Reflection characteristics of plane-S-wave at the free boundary of unsaturated porothermoelastic media[J]. Journal of Thermal Stresses, 2020, 43(5): 579-593.
[18] LIU H B, DAI G L, ZHOU F X, et al. A mixture theory analysis for reflection phenomenon of homogeneous plane-P1-wave at the boundary of unsaturated porothermoelastic media[J]. Geophysical Journal International, 2021, 228(2): 1237-1259.
[19] WEI W, ZHENG R Y, LIU G B, et al. Reflection and refraction of P wave at the interface between thermoelastic and porous thermoelastic medium[J]. Transport in Porous Media, 2016, 113(1): 1-27.
[20] HOU W T, FU L Y, CARCIONE J M, et al. Reflection and transmission of inhomogeneous plane waves in thermoelastic media[J]. Frontiers in Earth Science, 2022, 10: 850331.
[21] LIU H B, DAI G L, ZHOU F X, et al. Propagation behavior of homogeneous plane-P1-wave at the interface between a thermoelastic solid medium and an unsaturated porothermoelastic medium[J]. European Physical Journal Plus, 2021, 136(11): 1163.
[22] 李伟华, 郑洁. 饱和度对平面P波入射下自由场地地震反应的影响分析[J]. 岩土工程学报, 2017, 39(3): 427-435. doi: 10.11779/CJGE201703005 LI Weihua, ZHENG Jie. Effects of saturation on free-field responses of site due to plane P-wave incidence[J]. Chinese Journal of Geotechnical Engineering, 2017, 39(3): 427-435. (in Chinese) doi: 10.11779/CJGE201703005
[23] LI W H, ZHENG J, TRIFUNAC M D. Saturation effects on ground motion of unsaturated soil layer-bedrock system excited by plane P and SV waves[J]. Soil Dynamics and Earthquake Engineering, 2018, 110: 159-172.
[24] LIU H B, JIANG M, ZHOU F, et al. Attenuation characteristics of thermoelastic waves in unsaturated soil[J]. Arabian Journal of Geosciences, 2021, 14(18): 1-11.
[25] 郑洁. 非饱和土自由场地地震地面运动的解析分析[D]. 北京: 北京交通大学, 2015. ZHENG Jie. Analytical Analysis of Earthquake Ground Motion in Unsaturated Soil Free Field[D]. Beijing: Beijing Jiaotong University, 2015. (in Chinese)