Theoretical basis and application verification of scale effects of deformation characteristics of rockfill
-
摘要: 针对筑坝料变形特性的尺度效应,以接触变形理论为基础,从粒径尺度、试件尺度及外部影响因素角度,结合试验的方法开展了研究。建立了粒径尺度与归一化压缩模量的指数相关模型,论证了试件尺度径径比大于5时压缩试验的变形可靠性,揭示了侧限条件、干密度、边缘破碎及荷载等级5个外部因素对尺度效应的影响规律和内在机理,并分别建立了各因素与变形的相关模型。基于试验回溯现场过程,建立了综合考虑可破碎材料（堆石料）与不破碎材料（砂砾石）的多尺度统一修正模型。通过对大坝堆石料EB模型计算参数的修正及变形计算,验证了多尺度统一修正模型对变形计算的适用性与可靠性,较直接采用缩尺试验结果准确度提升了约20%,为解决坝体变形预测不准确提供了新的思考方法与解决路径。Abstract: In view of the scale effects of deformation characteristics of rockfill, based on the contact deformation theory, the research is carried out from the perspective of particle sizes, specimen sizes and external factors, combined with the test method. The exponential correlation model between the particle sizes and the normalized compression modulus is established. The deformation reliability of compression tests is demonstrated when the particle size ratio is greater than 5. The influence rules and internal mechanism of five external factors on the scale effects, namely, the lateral limit condition, the dry density, the edge breakage and the load level, are revealed, and the relevant models for each factor and deformation are established respectively. Based on the field experiments, a multi-scale unified correction model is established, which takes into account the crushed material (rockfill) and the non-crushed materials (gravel). The applicability and reliability of the multi-scale unified correction model for deformation calculation are verified by modifying the parameters and deformation calculation of the EB model for dam rockfill. The accuracy of the model is about 20% higher than that of the direct scale tests. This provides a new thinking method and solution path to solve the inaccuracy of dam deformation prediction.
-
Keywords:
- scale effect /
- compression modulus /
- influencing factor /
- correction model /
- dam deformation
-
0. 引言
中国在公路、铁路、水利及城市建设等多个行业中经常会遇到岩质边坡工程，工程上常采用二维极限平衡法进行岩质边坡的稳定性分析。事实上，岩质边坡失稳呈现显著的三维空间效应，二维极限平衡法已不能满足一些重要工程的精度要求。因此，开展岩质边坡三维极限平衡法的理论研究具有重要的理论意义和工程应用前景。
经典的三维极限平衡法是基于二维条分法的拓展, 陈祖煜[1]对早期研究代表性成果做了总结，相关研究奠定了三维极限平衡法的理论基础。朱大勇等[2-3]建立了基于滑面正应力修正的三维极限平衡法；Zheng[4]将三维极限平衡法归结为代数特征值问题，解决了安全系数的不收敛问题；近年来，Zhou等[5]提出了考虑条间力和基于位移的严格极限平衡法；邓东平等[6]提出了三维边坡稳定性分析的非严格法和严格法；高玉峰等[7]建立了基于一个方向力（力矩）平衡的三维非对称边坡稳定性分析方法；卢坤林等[8]建立了适用任意形态滑面的三维极限平衡法。上述研究成果都是基于Mohr-Coulomb（简称M-C）强度准则，适用于土质边坡。
朱合华等[9]认为M-C强度准则不能完全反映岩石力学特性。对于岩质边坡稳定性研究应选用Hoek-Brown（简称H-B）强度准则。孙超伟等[10]采用H-B强度折减法，编制了一套岩质边坡安全系数的稳定性图表；林杭等[11]实现了H-B强度折减法的三维岩质边坡稳定性分析方法；Deng等[12]采用非线性H-B强度准则的泰勒级数线性展开式实现了岩质边坡安全系数的计算；Kumar等[13]将H-B破坏准则与等效的M-C强度参数相结合，采用PSO分析岩质边坡的稳定性；Yang等[14]基于H-B强度准则与极限分析上限法，实现了动力和静力荷载下两级三维岩质边坡稳定性分析方法；Michalowski等[15]基于H-B强度准则，采用极限分析运动学方法对弧形滑动的岩质边坡进行稳定性评估。
基于M-C强度准则的土质边坡三维极限平衡法已有了长足的进展，但基于H-B强度准则的岩质边坡三维极限平衡法的研究尚不够充分。鉴于此，本文将H-B强度准则与构造滑面正应力分布的极限平衡法[8]相结合，合理构造滑面正应力分布函数，建立满足三个力平衡和一个绕垂直滑动方向的力矩平衡的岩质边坡三维极限平衡法。
1. 滑面正应力的构造
恰当合理地构造滑面正应力分布模式是边坡三维稳定性分析正确解答的前提。朱大勇等[2-3]假设了滑面正应力分布合理形式。卢坤林等[8]研究表明，边坡稳定性系数对合理构造的滑面正应力分布不敏感，建议构造的滑面正应力分布函数（坐标原点选在滑体中心）：
$$ \sigma = {\sigma _0}{\text{(}}1 + {\lambda _1} + {\lambda _2} \cdot x + {\lambda _3} \cdot y)。 $$ (1) 式中：${\lambda _1}$，${\lambda _2}$，${\lambda _3}$为待定参数；${\sigma _0}$为初始滑面正应力，采用Hovland[16]的瑞典法3D扩展的滑面正应力。
合理构造的滑面正应力必须满足：①滑面正应力非负，滑面正应力连续光滑且单值；②初始滑面正应力${\sigma _0}$为滑面正应力构造的核心。
2. 逐点等效M-C强度参数
2.1 H-B强度准则
H-B强度准则是Hoek等通过对大量岩石三轴试验资料和岩体现场试验成果的统计分析，并综合考虑了岩体结构、岩块强度、应力状态等多方面的影响，而得出的岩体破坏经验公式。在岩质边坡稳定性分析中，通常需要将强度准则表达成法向应力和切应力的关系，为此Hoek等[17]建立的表达式为
$$ {\tau _{\text{f}}} = A{\sigma _{\text{c}}}{\left( {\frac{{\sigma - {\sigma _{{\text{tm}}}}}}{{{\sigma _{\text{c}}}}}} \right)^B}。 $$ (2) 式中：A，B均为与岩体材料相关的拟合参数，通过对三轴应力点拟合确定[18]；${\sigma _{\text{c}}}$为岩石单轴抗压强度；$\sigma $为滑面的法向应力；${\sigma _{{\text{tm}}}}$为节理岩体的抗拉强度。
2.2 逐点等效M-C强度参数
如图 1所示，选取H-B强度包络线上的任意一点（σi，τi）做该点切线，其斜率的反正切值即为该点的φi值，其与y轴的截距即为该点的ci值。拓展到整条曲线，各点的φ值和c值可表示为
$$ \varphi = \arctan \left[ {AB{{\left( {\frac{{\sigma - {\sigma _{{\text{tm}}}}}}{{{\sigma _{\text{c}}}}}} \right)}^{B - 1}}} \right] \text{，} $$ (3) $$ c = A{\sigma _{\text{c}}}{\left( {\frac{{\sigma - {\sigma _{{\text{tm}}}}}}{{{\sigma _{\text{c}}}}}} \right)^B} - \sigma \tan \varphi 。 $$ (4) ![]()
图 1 H-B逐点等效M-C等效方式示意Figure 1. Schematic of H-B strength parameters by point-by-point equivalent M-C strength parameter式中：φ为逐点等效内摩擦角；c为逐点等效黏聚力，均为滑面正应力$\sigma $的函数，并与坐标位置密切关联。由于正应力$\sigma $的分布未知，逐点等效内摩擦角φ和逐点等效黏聚力c的值也不能直接获取。
引入构造滑面正应力的分布（式（1）），得到逐点等效内摩擦角$ \varphi $和逐点等效黏聚力c的表达式为
$$ \varphi =\mathrm{arctan}\left[AB{\left(\frac{{\sigma }_{0}（1+{\lambda }_{1}+{\lambda }_{2}\cdot x+{\lambda }_{3}\cdot y)-{\sigma }_{\text{tm}}}{{\sigma }_{\text{c}}}\right)}^{B-1}\right] \text{，} $$ (5) $$\begin{array}{c} c=A{\sigma }_{\text{c}}{\left(\frac{{\sigma }_{0}（1+{\lambda }_{1}+{\lambda }_{2}\cdot x+{\lambda }_{3}\cdot y)-{\sigma }_{\text{tm}}}{{\sigma }_{\text{c}}}\right)}^{B}- \\ {\sigma }_{0}（1+{\lambda }_{1}+{\lambda }_{2}\cdot x+{\lambda }_{3}\cdot y)\mathrm{tan}\varphi 。 \end{array}$$ (6) 由式（5），（6）能够获得整个滑面上每个点的等效等效内摩擦角φ和等效黏聚力c，在空间上分布不恒定，这与采用Hoek等[19]提出的常规等效M-C参数为固定值有着显著的区别。
3. 岩质边坡三维极限平衡法
3.1 基本假定
（1）滑动方向与oxz面平行，即滑面剪应力与滑动方向平行。
（2）图 2（a）所示的边坡与三维滑面，坡面方程z=g，滑面方程为z=S，Sx和Sy为滑面S沿x，y轴的斜率。
（3）图 2（b）为典型条柱，作用其上的重力W=wdxdy，w为滑体单位面积上的条柱重量，岩体重度为γ，$ w = \gamma \cdot (g - S) $。条块底面法向力N=$ \sigma \mathit{\Delta }{\text{d}}x{\text{d}}y $和切向力T=$ \tau \mathit{\Delta }{\text{d}}x{\text{d}}y $，其中$ {\mathit{\Delta } = }\sqrt {1 + S_x^2 + S_y^2} $。滑面正应力分布采用式（1），其中${\sigma _0} = w/(1 + S_x^2 + S_y^2)$，为条柱体积力对滑面正应力的贡献分量，也就是采用Hovland[16]的瑞典法3D扩展的滑面正应力。
3.2 平衡方程
考虑整个滑体的3个方向的力平衡和绕垂直方向（y轴）的力矩平衡，忽略另两个力矩平衡条件可得：
$$ \iint \sigma \cdot n_\sigma ^y \cdot \mathit{\Delta }{\text{d}}x{\text{d}}y = 0 \text{，} $$ (7a) $$ {\displaystyle \iint （\sigma }\cdot {n}_{\sigma }^{x}+\tau \cdot {n}_{\tau }^{x})\cdot \mathit{\Delta }\text{d}x\text{d}y=0 \text{，} $$ (7b) $$ {\displaystyle \iint （\sigma }\cdot {n}_{\sigma }^{z}+\tau \cdot {n}_{\tau }^{z})\cdot \mathit{\Delta }\text{d}x\text{d}y-W=0 \text{，} $$ (7c) $$\begin{array}{l} {\displaystyle \iint （\sigma }\cdot {n}_{\sigma }^{z}+\tau \cdot {n}_{\tau }^{z})\cdot x\cdot \mathit{\Delta }\text{d}x\text{d}y- \\ \iint {(\sigma \cdot n_\sigma ^x + \tau \cdot n_\tau ^x)} \cdot S \cdot \mathit{\Delta }{\text{d}}x{\text{d}}y = M 。 \end{array}$$ (7d) 式中：$ W = \iint w \cdot {\text{d}}x{\text{d}}y $；$ M = \iint w \cdot x \cdot {\text{d}}x{\text{d}}y $；$ n_\sigma ^x = - \frac{{{S_x}}}{\mathit{\Delta }} $；$ n_\sigma ^y = - \frac{{{S_y}}}{\mathit{\Delta }} $；$ n_\sigma ^z = \frac{1}{\mathit{\Delta }} $; $ n_\tau ^x = \frac{1}{{{\mathit{\Delta } '}}} $；$ n_\tau ^y = 0 $；$ n_\tau ^z = \frac{{{S_x}}}{{{\mathit{\Delta } '}}} $；$ {\mathit{\Delta } '} = \sqrt {1 + S_x^2} $。
根据H-B准则（式（2））有：
$$ \tau = \frac{{{\tau _{\text{f}}}}}{{{F_{\text{s}}}}} = \frac{{A{\sigma _{\text{c}}}{{\left( {\frac{{\sigma - {\sigma _{{\text{tm}}}}}}{{{\sigma _{\text{c}}}}}} \right)}^B}}}{{{F_{\text{s}}}}}。 $$ (8) 由式（8）表示的τ在代入式（7a）~（7d）后，由于其指数B的取值为小于1的非整数，使得与τ相关的未知量无法提出而全部堆叠在方程组右侧。经多次演算发现，在后续迭代过程中均出现了难以收敛的问题。Deng等[12]也指出在计算过程中由于要嵌套多个公式使得迭代过程变得困难。故采用逐点等效M-C强度参数（式（5），（6））：
$$ \tau = \frac{{c + \sigma \tan \varphi }}{{{F_{\text{s}}}}} 。 $$ (9) 式中：c，φ为逐点等效黏聚力和逐点等效内摩擦角，是滑面正应力的变量，随滑面位置和滑面正应力逐点变化，并非是固定值。
将式（9）代入式（7a）~（7d），进一步整理得：
$$\begin{array}{c} {\lambda _1}\iint {( - {S_y}){\sigma _0}}{\text{d}}x{\text{d}}y{\text{ + }}{\lambda _2}\iint {( - {S_y})}{\sigma _0}x{\text{d}}x{\text{d}}y + \\{\lambda _3}\iint {( - {S_y})}{\sigma _0}y{\text{d}}x{\text{d}}y = \iint {{S_y}{\sigma _0}}{\text{d}}x{\text{d}}y \text{，} \end{array}$$ (10a) $$\begin{array}{c} {\lambda _1}\iint {\left( { - {S_x} + \frac{{\mathit{\Delta }\tan \varphi }}{{{\mathit{\Delta } '}{F_{\text{s}}}}}} \right){\sigma _0}}{\text{d}}x{\text{d}}y{\text{ + }}{\lambda _2}\iint {\left( { - {S_x} + \frac{{\mathit{\Delta }\tan \varphi }}{{{\mathit{\Delta } '}{F_{\text{s}}}}}} \right){\sigma _0}}x{\text{d}}x{\text{d}}y + \\ {\lambda _3}\iint {\left( { - {S_x} + \frac{{\mathit{\Delta }\tan \varphi }}{{{\mathit{\Delta } '}{F_{\text{s}}}}}} \right){\sigma _0}y{\text{d}}x{\text{d}}y} = - \iint c \cdot \frac{\mathit{\Delta }}{{{\mathit{\Delta } '}{F_{\text{s}}}}}{\text{d}}x{\text{d}}y - \\ \iint {\left( { - {S_x} + \frac{{\mathit{\Delta }\tan \varphi }}{{{\mathit{\Delta } '}{F_{\text{s}}}}}} \right){\sigma _0}}{\text{d}}x{\text{d}}y \text{，} \end{array}$$ (10b) $$\begin{array}{c} {\lambda _1}\iint {\left( {1 + {S_x}\frac{{\mathit{\Delta }\tan \varphi }}{{{\mathit{\Delta } '}{F_{\text{s}}}}}} \right){\sigma _0}}{\text{d}}x{\text{d}}y{\text{ + }}{\lambda _2}\iint {\left( {1 + {S_x}\frac{{\mathit{\Delta }\tan \varphi }}{{{\mathit{\Delta } '}{F_{\text{s}}}}}} \right){\sigma _0}}x{\text{d}}x{\text{d}}y + \\ {\lambda _3}\iint {\left( {1 + {S_x}\frac{{\mathit{\Delta }\tan \varphi }}{{{\mathit{\Delta } '}{F_{\text{s}}}}}} \right){\sigma _0}y{\text{d}}x{\text{d}}y} = W - \iint c \cdot {S_x}\frac{\mathit{\Delta }}{{{\mathit{\Delta } '}{F_{\text{s}}}}}{\text{d}}x{\text{d}}y - \\ \iint {\left( {1 + {S_x}\frac{{\mathit{\Delta }\tan \varphi }}{{{\mathit{\Delta } '}{F_{\text{s}}}}}} \right){\sigma _0}}{\text{d}}x{\text{d}}y \text{，} \end{array}$$ (10c) $$ \begin{array}{c} {F_{\text{s}}} = \left[ {\iint {{\sigma _0}(1 + {\lambda _1}} + {\lambda _2}x + {\lambda _3}y) \cdot ({S_x} \cdot x - S) \cdot \frac{\mathit{\Delta }}{{{\mathit{\Delta } '}}}\tan \varphi {\text{d}}x{\text{d}}y + } \right. \hfill \\ \left. {\iint c \cdot ({S_x} \cdot x - S)\frac{\mathit{\Delta }}{{{\mathit{\Delta } '}}}{\text{d}}x{\text{d}}y} \right]/\left[ {M - \iint {{\sigma _0}(1 + {\lambda _1}} + {\lambda _2}x + {\lambda _3}y) \cdot } \right. \hfill \\ \left. {\mathop {}\limits^{} (x + {S_x} \cdot S){\text{d}}x{\text{d}}y} \right] 。 \end{array} $$ (10d) 将式（10a）~（10d）简记为
$$ {\lambda _1}{D_1} + {\lambda _2}{D_2} + {\lambda _3}{D_3} = {D_4} \text{，} $$ (11a) $$\begin{array}{c} {\lambda _1}({F_{\text{s}}}{A_{11}} + {A'_{11}}) + {\lambda _2}({F_{\text{s}}}{A_{12}} + {A'_{12}}) + {\lambda _3}({F_{\text{s}}}{A_{13}} + {A'_{13}})\\ = {F_{\text{s}}}{A_{14}} + {A'_{14}} \text{，} \end{array}$$ (11b) $$\begin{array}{c} {\lambda _1}({F_{\text{s}}}{A_{21}} + {A'_{21}}) + {\lambda _2}({F_{\text{s}}}{A_{22}} + {A'_{22}}) + {\lambda _3}({F_{\text{s}}}{A_{23}} + {A'_{23}})\\ = {F_{\text{s}}}{A_{24}} + {A'_{24}} \text{，} \end{array}$$ (11c) $$\begin{array}{c} {\lambda _1}({F_{\text{s}}}{A_{31}} + {A'_{31}}) + {\lambda _2}({F_{\text{s}}}{A_{32}} + {A'_{32}}) + {\lambda _3}({F_{\text{s}}}{A_{33}} + {A'_{33}})\\ = {F_{\text{s}}}{A_{34}} + {A'_{34}} 。 \end{array}$$ (11d) 式（11a）~（11d）中的简写符号的表达式见表 1。
表 1 符号对应具体表达式Table 1. Symbols corresponding to specific expressions简写符号 表达式 简写符号 表达式 简写符号 表达式 $ {D_1} $ $ \iint { - {S_y}{\sigma _0}{\text{d}}x{\text{d}}y} $ $ {A_{14}} $ $ \iint {{S_x}}{\sigma _0}{\text{d}}x{\text{d}}y $ $ {A_{31}} $ $ \iint {{\sigma _0}(x + }S \cdot {S_x}){\text{d}}x{\text{d}}y $ $ {D_2} $ $ \iint { - {S_y}{\sigma _0}x{\text{d}}x{\text{d}}y} $ $ {A'_{14}} $ $ \iint { - (c + {\sigma _0}\tan \varphi )\frac{\mathit{\Delta }}{{{\mathit{\Delta } '}}}}{\text{d}}x{\text{d}}y $ $ {A'_{31}} $ $ \iint {(x{S_x} - S)\frac{\mathit{\Delta }}{{{\mathit{\Delta } '}}}\tan \varphi {\sigma _0}}{\text{d}}x{\text{d}}y $ $ {D_3} $ $ \iint { - {S_y}{\sigma _0}y{\text{d}}x{\text{d}}y} $ $ {A_{21}} $ $ \iint {{\sigma _0}{\text{d}}x{\text{d}}y} $ $ {A_{32}} $ $ \iint {{\sigma _0}(x + }S \cdot {S_x})x{\text{d}}x{\text{d}}y $ $ {D_4} $ $ \iint {{S_y}{\sigma _0}{\text{d}}x{\text{d}}y} $ $ {A'_{21}} $ $ \iint {{S_x}\frac{\mathit{\Delta }}{{{\mathit{\Delta } '}}}\tan \varphi }{\sigma _0}{\text{d}}x{\text{d}}y $ $ {A'_{32}} $ $ \iint {(x{S_x} - S)\frac{\mathit{\Delta }}{{{\mathit{\Delta } '}}}\tan \varphi {\sigma _0}}x{\text{d}}x{\text{d}}y $ $ {A_{11}} $ $ \iint { - {S_x}}{\sigma _0}{\text{d}}x{\text{d}}y $ $ {A_{22}} $ $ \iint {{\sigma _0}x{\text{d}}x{\text{d}}y} $ $ {A_{33}} $ $ \iint {{\sigma _0}(x + }S \cdot {S_x})y{\text{d}}x{\text{d}}y $ $ {A'_{11}} $ $ \iint {\frac{\mathit{\Delta }}{{{\mathit{\Delta } '}}}\tan \varphi }{\sigma _0}{\text{d}}x{\text{d}}y $ $ {A'_{22}} $ $ \iint {{S_x}\frac{\mathit{\Delta }}{{{\mathit{\Delta } '}}}\tan \varphi }{\sigma _0}x{\text{d}}x{\text{d}}y $ $ {A'_{33}} $ $ \iint {(x{S_x} - S)\frac{\mathit{\Delta }}{{{\mathit{\Delta } '}}}\tan \varphi {\sigma _0}}y{\text{d}}x{\text{d}}y $ $ {A_{12}} $ $ \iint { - {S_x}}{\sigma _0}x{\text{d}}x{\text{d}}y $ $ {A_{23}} $ $ \iint {{\sigma _0}y{\text{d}}x{\text{d}}y} $ $ {A_{34}} $ $ M - \iint {{\sigma _0}(x + }S \cdot {S_x}){\text{d}}x{\text{d}}y $ $ {A'_{12}} $ $ \iint {\frac{\mathit{\Delta }}{{{\mathit{\Delta } '}}}\tan \varphi }{\sigma _0}x{\text{d}}x{\text{d}}y $ $ {A'_{23}} $ $ \iint {{S_x}\frac{\mathit{\Delta }}{{{\mathit{\Delta } '}}}\tan \varphi }{\sigma _0}y{\text{d}}x{\text{d}}y $ $ {A'_{34}} $ $ \iint { - (c + {\sigma _0}\tan \varphi )(x{S_x} - S)\frac{\mathit{\Delta }}{{{\mathit{\Delta } '}}}}{\text{d}}x{\text{d}}y $ $ {A_{13}} $ $ \iint { - {S_x}}{\sigma _0}y{\text{d}}x{\text{d}}y $ $ {A_{24}} $ $ W - \iint {{\sigma _0}{\text{d}}x{\text{d}}y} $ $ {A'_{13}} $ $ \iint {\frac{\mathit{\Delta }}{{{\mathit{\Delta } '}}}\tan \varphi }{\sigma _0}y{\text{d}}x{\text{d}}y $ $ {A'_{24}} $ $ \iint { - (c + {\sigma _0}\tan \varphi ){S_x}\frac{\mathit{\Delta }}{{{\mathit{\Delta } '}}}}{\text{d}}x{\text{d}}y $ 3.3 稳定性系数迭代求解
求解稳定性系数步骤如下：
（1）假定初始的λ1(0)，λ2(0)，λ3(0)，代入式（1）得到初始滑面正应力函数。
（2）将式（1）的滑面正应力函数代入式（5）、（6）得到逐点等效黏聚力和逐点等效内摩擦角的初始表达式，再将其代入式（11a）~（11d）。由克拉默法则通过式（11b）~（11d）解得：
$$ \lambda _i^{(1)} = \frac{{{\mathit{\Delta }_i}^{(0)}}}{{{\mathit{\Delta }_0}^{(0)}}}{\text{ }}(i = 1,2,3) 。 $$ (12) 式中：
$$ {\mathit{\Delta }_0}^{(0)} = \left| {\begin{array}{*{20}{c}} {{F_{\text{s}}}{A_{11}} + {{A'}_{11}}}&{{F_{\text{s}}}{A_{12}} + {{A'}_{12}}}&{{F_{\text{s}}}{A_{13}} + {{A'}_{13}}} \\ {{F_{\text{s}}}{A_{21}} + {{A'}_{21}}}&{{F_{\text{s}}}{A_{22}} + {{A'}_{22}}}&{{F_{\text{s}}}{A_{23}} + {{A'}_{23}}} \\ {{F_{\text{s}}}{A_{31}} + {{A'}_{31}}}&{{F_{\text{s}}}{A_{32}} + {{A'}_{32}}}&{{F_{\text{s}}}{A_{33}} + {{A'}_{33}}} \end{array}} \right| \text{，} $$ (13a) $$ {\mathit{\Delta }_1}^{(0)} = \left| {\begin{array}{*{20}{c}} {{F_{\text{s}}}{A_{11}} + {{A'}_{14}}}&{{F_{\text{s}}}{A_{12}} + {{A'}_{12}}}&{{F_{\text{s}}}{A_{13}} + {{A'}_{13}}} \\ {{F_{\text{s}}}{A_{21}} + {{A'}_{24}}}&{{F_{\text{s}}}{A_{22}} + {{A'}_{22}}}&{{F_{\text{s}}}{A_{23}} + {{A'}_{23}}} \\ {{F_{\text{s}}}{A_{31}} + {{A'}_{34}}}&{{F_{\text{s}}}{A_{32}} + {{A'}_{32}}}&{{F_{\text{s}}}{A_{33}} + {{A'}_{33}}} \end{array}} \right| \text{，} $$ (13b) $$ {\mathit{\Delta }_2}^{(0)} = \left| {\begin{array}{*{20}{c}} {{F_{\text{s}}}{A_{11}} + {{A'}_{11}}}&{{F_{\text{s}}}{A_{12}} + {{A'}_{14}}}&{{F_{\text{s}}}{A_{13}} + {{A'}_{13}}} \\ {{F_{\text{s}}}{A_{21}} + {{A'}_{21}}}&{{F_{\text{s}}}{A_{22}} + {{A'}_{24}}}&{{F_{\text{s}}}{A_{23}} + {{A'}_{23}}} \\ {{F_{\text{s}}}{A_{31}} + {{A'}_{31}}}&{{F_{\text{s}}}{A_{32}} + {{A'}_{34}}}&{{F_{\text{s}}}{A_{33}} + {{A'}_{33}}} \end{array}} \right| \text{，} $$ (13c) $$ {\mathit{\Delta }_3}^{(0)} = \left| {\begin{array}{*{20}{c}} {{F_{\text{s}}}{A_{11}} + {{A'}_{11}}}&{{F_{\text{s}}}{A_{12}} + {{A'}_{12}}}&{{F_{\text{s}}}{A_{13}} + {{A'}_{14}}} \\ {{F_{\text{s}}}{A_{21}} + {{A'}_{21}}}&{{F_{\text{s}}}{A_{22}} + {{A'}_{22}}}&{{F_{\text{s}}}{A_{23}} + {{A'}_{24}}} \\ {{F_{\text{s}}}{A_{31}} + {{A'}_{31}}}&{{F_{\text{s}}}{A_{32}} + {{A'}_{32}}}&{{F_{\text{s}}}{A_{33}} + {{A'}_{34}}} \end{array}} \right| 。 $$ (13d) 将式（12）代入式（11a），经整理可得关于稳定性系数的一元三次方程：
$$ {t_3}F_{\text{s}}^3 + {t_2}{F_{\text{s}}}^2 + {t_1}{F_{\text{s}}} + {t_0} = 0 。 $$ (14) 式中：t0，t1，t2，t3为系数，已知。其具体表达式可参考文献[3]，其中简写表达式不同处已在表 1中给出。求解式（14）的最大实根[3]为稳定性系数Fs(1)，将其代入式（12）即可获得迭代后的λ1(1)，λ2(1)，λ3(1)。
（3）将λ1(1)，λ2(1)，λ3(1)代回式（1），得到迭代后的滑面正应力分布函数，并重复（2）的步骤可得稳定性系数Fs(2)，λ1(2)，λ2(2)，λ3(2)。
（4）比较Fs(1)和Fs(2)，λ1(1)和λ1(2)，λ2(1)和λ2(2)，λ3(1)和λ3(2)，若二者间的差值均满足精度要求（ε=0.001），则Fs=（Fs(1)+Fs(2)）/2，否则令λ1(1)=λ1(2)，λ2(1)=λ2(2)，λ3(1)=λ3(2)，重复步骤（2），直至前次和本次稳定性系数与待定参数的差值满足精度要求（ε=0.001）。
具体计算流程如图 3所示。
4. 算例验证
4.1 算例1
如图 4所示，算例1的岩质边坡来自文献[20]，破坏形式为楔体破坏。
采用直角坐标表示的三维滑面数学方程为
$$ \left.\begin{array}{l}{S}_{ABC}=\frac{349.36}{312.23}x+\frac{303.69}{312.23}y+\frac{4706.4}{312.23}\text{，}\\ {S}_{ABD}=\frac{246.62}{205.65}x-\frac{365.62}{205.65}y+\frac{422.3}{205.65}。\end{array}\right\} $$ (15) 岩体强度参数如表 2所示，mb，s，a均为反映岩体质量和岩体类型的常数，具体的表达式见文献[19]。
表 2 岩体参数Table 2. Material parameters of rock mass参数 滑面ABC 滑面ABD 重度γ/(kN·m-3) 25.0 25.0 单轴抗压强度σci/MPa 0.818 0.682 完整岩石材料参数mi 20 15 地质强度指标GSI 100 75 扰动因子D 0 0 mb 20 6.142 s 1 6.22×10-2 a 0.5 0.501 根据Hoek[19]方法计算滑面ABC、滑面ABD的岩体抗拉强度σtm分别为40.82，6.89 kPa。依据陈祖煜等[18]的数据拟合方法计算滑面ABC、滑面ABD的参数A分别为1.1719，0.8208；参数B分别为0.7151，0.7098。采用本文方法计算相应的稳定性系数，并与其他学者成果相比较，结果列于表 3中。
表 3 算例1稳定性系数计算结果Table 3. Calculated results of stability coefficient of example 1由表 3可知，采用本文计算方法得到的稳定性系数为3.562与广义H-B（Deng[20]）的结果基本一致，误差为0.86%，验证了本文方法的正确性与合理性。
4.2 算例2
如图 5所示，有一高度为30 m，坡比为1.5的岩质边坡，岩体的强度参数如表 4所示。
表 4 岩体参数及H-B计算参数Table 4. Parameters of rock mass and Hoek-Brown criterion参数 取值 重度γ/(kN·m-3) 25.0 单轴抗压强度σci/MPa 0.4 完整岩石材料参数mi 8 地质强度指标GSI 60 扰动因子D 0 mb 1.917 s 1.17×10-2 a 0.503 σtm/kPa 2.44 A 0.5630 B 0.6933 三维滑裂面的形态已知，采用直角坐标表示的三维滑面数学方程为
$$ \left.\begin{array}{l}{S}_{ABC}=\frac{3}{4}x+\frac{15}{14}y+15\text{，}\\ {S}_{ABD}=\frac{3}{4}x-\frac{5}{4}y+15。\end{array}\right\} $$ (16) 根据Hoek等[19]方法得到常规等效M-C强度准则的黏聚力c=54.77 kPa和内摩擦角φ=20.23°。分别采用逐点等效M-C和常规等效M-C计算该边坡的稳定性系数，结果列于表 5。
表 5 算例2稳定性系数计算结果Table 5. Calculated results of stability coefficient of example 2计算方法 计算结果 误差/% 逐点等效M-C(本文方法) 1.614 常规等效M-C(卢坤林等[8]) 1.913 15.63 常规等效M-C（三维楔形体法） 1.921 15.98 由表 5可知，本文方法得到的稳定性系数为1.614，低于三维楔形体法计算的稳定性系数1.921，偏低达15.98%；低于卢坤林等[8]的1.913，偏低15.63%。这是由于本文方法采用逐点等效M-C的强度参数，而三维楔形体法和卢坤林等[8]采用常规等效M-C的强度参数。这与Deng等[20]所得到的结论较为一致。
图 6为算例2的计算迭代过程，收敛速度快且收敛效果理想。图 7~9展示了收敛后的构造滑面正应力、逐点等效黏聚力和逐点等效内摩擦角的空间分布。
![]()
图 9 逐点等效内摩擦角分布Figure 9. Distribution of point-by-point equivalent internal friction angle4.3 算例3
如图 10所示，岩质边坡高度为15 m，坡比为0.5。岩体的强度参数同算例2。
三维滑裂面可采用直角坐标表示为
$$ \left.\begin{array}{l}{S}_{ABFG}=\frac{5}{22}x+\frac{75}{11}\text{，}\\ {S}_{ACF}=\frac{25}{38}x-\frac{18}{19}y+\frac{195}{19}\text{，}\\ {S}_{CEF}=\frac{10}{13}x+\frac{10}{13}y+\frac{145}{13}\text{，}\\ {S}_{EFG}=\frac{5}{9}x+\frac{85}{9}\text{，}\\ {S}_{DEG}=\frac{60}{73}x-\frac{50}{73}y+\frac{245}{73}\text{，}\\ {S}_{BDG}=\frac{35}{58}x-\frac{30}{29}y-\frac{75}{29}。\end{array}\right\} $$ (17) 根据Hoek等[19]方法得到常规等效M-C强度的黏聚力c=37.21kPa、内摩擦角φ=24.78°。采用本文方法计算的稳定性系数为2.547，常规等效M-C的三维极限平衡法为2.922，本文方法偏低12.83%，计算结果列于表 6。
表 6 算例3稳定性系数计算结果Table 6. Calculated results of stability coefficient of example 3计算方法 计算结果 误差/% 逐点等效M-C(本文方法) 2.547 常规等效M-C（朱大勇等[3]） 2.922 12.83 结合表 5，6可知，与常规等效M-C强度参数相比，采用逐点等效M-C强度参数进行岩质边坡三维稳定性分析会显著偏低一些。
4.4 算例4
如图 11所示，岩质边坡高20 m，坡比为1，岩体强度参数如表 7所示，滑面为椭球状，采用直角坐标表示的三维滑面数学方程为
$$ \frac{（x+6.9881{)}^{2}}{{35}^{2}}+\frac{{y}^{2}}{{(35*n)}^{2}}+\frac{{(z-34.2953)}^{2}}{{35}^{2}}=1 。 $$ (18) 表 7 岩体参数Table 7. Material parameters of rock mass参数 取值 重度γ/(kN·m-3) 23.0 单轴抗压强度σci/MPa 0.081 完整岩石材料参数mi 15 地质强度指标GSI 70 扰动因子D 0 mb 5.138 s 3.57×10-2 a 0.501 σtm/kPa 0.842 A 0.7771 B 0.7101 式中：n为该滑体y轴轴长与x轴轴长的比值，通过n的变化，可以表征不同形态的椭球状滑面。
根据Hoek等[19]方法得到常规等效M-C强度准则的黏聚力c=46.25 kPa和内摩擦角φ=19.78°。分别采用逐点等效M-C和常规等效M-C计算该边坡的稳定性系数，结果列于表 8。
表 8 算例4稳定性系数计算结果Table 8. Calculated results of stability coefficient of example 4由表 8可知，本文方法得到的稳定性系数低于常规等效M-C所得到的稳定性系数，偏低8%左右；在滑面趋近于二维形态时，本文方法计算得到的稳定性系数更接近二维稳定性系数[21]，再次验证了本文方法的合理性。
5. 讨论
5.1 逐点等效M-C与常规等效M-C的参数差异
逐点等效M-C强度参数是滑面正应力的变量，常规等效M-C强度参数为固定值。从算例2和算例3可看出，本文方法得到的等效黏聚力明显低于常规等效的黏聚力，平均值分别偏低56.9%，46.6%；而本文方法得到的等效内摩擦角明显高于常规等效的内摩擦角，平均值分别偏高75.1%，52.8%。
由于逐点等效M-C强度参数与常规等效M-C强度参数有较大的差异，当采用同一个三维稳定性分析方法计算稳定性系数时，本文方法所得稳定性系数显著偏低。
5.2 滑面正应力分布对稳定性系数的影响
滑面正应力σ的分布对逐点等效M-C强度参数有着直接的影响，进而影响到稳定性系数的大小。为探究这一影响幅度，借鉴刘华丽等[22]的研究思路，在式（1）的基础上引入变量k，此时的正应力表达式为
$$ \sigma ={\sigma }_{0}（1+{\lambda }_{1}+{\lambda }_{2}\cdot x+{\lambda }_{3}\cdot y+k\cdot xy) 。 $$ (19) 以算例2为例，参数k的取值需保证滑面正应力处处为正，经试算k的取值范围为-0.03~0.03。k在不同取值下相对应的稳定性系数如图 12，由图 12可看出滑面正应力分布的不同对边坡稳定性系数的影响不显著，最大差别仅为1.80%。
![]()
图 12 不同正应力分布对应的稳定性系数Figure 12. Stability coefficient corresponding to different normal stress distributions6. 工程实例
6.1 工程概况
如图 13所示，该工程为一矿区边坡，中部冲沟处为不稳定坡体。不稳定坡体后缘以后缘共轭光滑节理面为界，不稳定坡体左右两侧主要以出露的碎裂状凝灰岩结构为界；前缘主要根据地形陡缓交界及开挖的矿坑边坡综合推断确定。不稳定坡体纵长约190 m，最大横宽约55 m，平面面积约1.7×104 m2，根据不稳定坡体前、中部探槽揭示，滑体厚度5.5~30.5 m，平均厚度约18.8 m，总体积约1.9×105 m3。
6.2 计算模型
根据现场实测数据，建立三维边坡计算模型如图 14所示，图 14（a）为三维边坡形态，图 14（b）为潜在三维滑面的形态。由于三维滑面的空间分布不规律，难以准确地采用数学方程来表示，通过Matlab读取三维滑面模型来获取滑面上各点的坐标，采用坐标矩阵代替数学方程来获取三维滑面。
岩体的强度参数如表 9所示。
表 9 岩体参数Table 9. Material parameters of rock mass参数 取值 重度γ/(kN·m-3) 28.0 单轴抗压强度σci /MPa 10 完整岩石材料参数mi 6 地质强度指标GSI 26 扰动因子D 0.8 mb 0.0733 s 1.4×10-5 a 0.529 根据Hoek等[19]方法计算得岩体抗拉强度σtm= 1.842 kPa。采用陈祖煜等[18]的数据拟合方法计算参数A，B得出A=0.1406，B=0.6734。
6.3 计算结果
采用本文方法对该不稳定坡体进行稳定性分析，并与其他三维极限平衡法相比较，列于表 10中。
表 10 不稳定坡体稳定性系数计算结果Table 10. Calculated results of stability coefficient of landslide area计算方法 计算结果 误差/% 逐点等效M-C(本文方法) 0.942 常规等效M-C（朱大勇等[3]） 1.154 18.37 由表 10可知，本文计算方法得到该不稳定坡体的稳定性系数为0.942，低于朱大勇等[3]的1.154，偏低18.37%。实际工程中，该边坡已经临近失稳（稳定系数预计约为1.0），从计算结果来看，朱大勇等[3]方法得到的稳定系数偏高，结果偏于激进；而本文方法得到的稳定系数偏低，结果偏于保守一些。最后，业主采纳了本文方法结果。本文方法偏低的原因可能是H-B模型参数取值有偏差造成的。
计算迭代过程如图 15，图 16~18为迭代稳定后的构造滑面正应力、逐点等效黏聚力和逐点等效内摩擦角空间分布。
![]()
图 18 逐点等效内摩擦角分布Figure 18. Distribution of point-by-point equivalent internal friction angel7. 结论
针对岩质边坡中岩体强度包线呈非线性分布的特征，提出了一种基于构造滑面正应力分布的岩质边坡三维极限平衡法，该方法收敛性较好，收敛速度快，易于编程实现，适用于任意滑面形态。
（1）通过构造滑面正应力分布，提出了逐点等效M-C强度参数，将等效黏聚力和等效内摩擦角转化为滑面正应力的变量，获得了逐点等效黏聚力和逐点等效内摩擦角的空间分布。与常规等效M-C强度参数相比，本文得到的等效参数不再是固定值，逐点等效黏聚力低于常规等效的黏聚力，而逐点等效内摩擦角高于常规等效的内摩擦角。
（2）将H-B强度准则的逐点等效M-C强度参数与构造滑面正应力相结合，根据整个滑体的平衡条件建立了4个平衡方程（3个力平衡方程与1个力矩平衡方程），建立了岩质三维边坡稳定性分析方法。4个算例和1个工程实例验证了本方法的正确性和合理性，本文方法得到的三维岩质边坡稳定性系数可能会偏于保守些。对于考虑整个滑体6个平衡方程的严格三维极限平衡法和最新版的H-B模型将是未来研究方向。
-
![]()
图 7 压缩模量与径径比的关系
Figure 7. Relationship between compression modulus and diameter diameter ratio
![]()
图 9 不同粒径下试样破碎变形计算
Figure 9. Calculation of crushing deformation of samples with different particle sizes
![]()
图 10 粒径与破碎变形量的关系
Figure 10. Relationship between particle size and crushing deformation
![]()
图 11 荷载等级与压缩模量的倍数关系
Figure 11. Relationship between load grade and compression modulus
表 1 修正前后河口村主堆石材料E-B模型参数
Table 1 E-B model parameters of main rockfill materials in Hekou Village
分类 干重度/(kN·m-3) K n Rf Kur c/kPa φ0/(°) Δφ/(°) Kb m 修正前 20.7 1660 0.21 0.85 3320 0 54 10.6 380 0.14 修正后 20.7 1312 0.25 0.82 3320 0 54 10.6 581 0.14 
下载: 导出CSV
表 2 河口村主堆石材料修正计算
Table 2 Correction calculation of main rockfill materials in Hekou Village
项目 f1(ξ) f2(μ) f3(ρd) f4(D) f5(p) 综合影响F(·) 修正值 1.19 1.0 1.47 0.56 1.0 0.84
下载: 导出CSV
表 3 河口村坝体沉降计算对比
Table 3 Comparison of settlement calculation of Hekou Village Dam
项目 计算最大沉降/cm 提升准确度/% 竣工 蓄水 室内缩尺试验 61.2 68.1 监测数据 78.9 109.7 17.43~28.89 统一修正模型 75.6 99.7
下载: 导出CSV
-
[1] IJHD. 2018 World Atlas & Industry Guide[R]. Wallington: International Journal of Hydropower and Dams, 2018.
[2] MA H Q, CHI F D. Technical progress on researches for the safety of high CFRDs[J]. Engineering, 2016, 2(4): 332-339. doi: 10.3969/j.issn.1009-444X.2016.04.011
[3] 秦红玉, 刘汉龙, 高玉峰, 等. 粗粒料强度和变形的大型三轴试验研究[J]. 岩土力学, 2004, 25(10): 1575-1580. doi: 10.3969/j.issn.1000-7598.2004.10.013 QIN Hong-yu, LIU Han-long, GAO Yu-feng, et al. Research on strength and deformation behavior of coarse aggregates based on large-scale triaxial tests[J]. Rock and Soil Mechanics, 2004, 25(10): 1575-1580. (in Chinese) doi: 10.3969/j.issn.1000-7598.2004.10.013
[4] 徐志华, 孙大伟, 张国栋. 堆石料应力-应变特性大型三轴试验研究[J]. 岩土力学, 2017, 38(6): 1565-1572. XU Zhi-hua, SUN Da-wei, ZHANG Guo-dong. Study on stress-strain behavior of rockfill using large-scale triaxial tests[J]. Rock and Soil Mechanics, 2017, 38(6): 1565-1572. (in Chinese)
[5] KONG X, LIU J, ZOU D, et al. Stress-dilatancy relationship of Zipingpu gravel under cyclic loading in triaxial stress states[J]. International Journal of Geomechanics, 2016, 16(4): 4016001. doi: 10.1061/(ASCE)GM.1943-5622.0000584
[6] MARACHI N D, CHAN C K, SEED H B. Evaluation of properties of rockfill materials[J]. Journal of Soil Mechanics & Foundations Div, 1972, 98(1): 95-114.
[7] HUNTER G, FELL R. Rockfill modulus and settlement of concrete face rockfill dams[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2003, 129(10): 909-917. doi: 10.1061/(ASCE)1090-0241(2003)129:10(909)
[8] VARADARAJAN A, SHARMA K G, ABBAS S M, et al. Constitutive model for rockfill materials and determination of material constants[J]. International Journal of Geomechanics, 2006, 6(4): 226-237. doi: 10.1061/(ASCE)1532-3641(2006)6:4(226)
[9] VARADARAJAN A, SHARMA K G, VENKATACHALAM K, et al. Testing and modeling two rockfill materials[J]. Journal of Geotechnical &Geoenvironmental Engineering, 2003, 129(3): 206-218.
[10] 汪小刚. 高土石坝几个问题探讨[J]. 岩土工程学报, 2018, 40(2): 203-222. WANG Xiao-gang. Discussion on some problems observed in high earth-rockfill dams[J]. Chinese Journal of Geotechnical Engineering, 2018, 40(2): 203-222. (in Chinese)
计量
- 文章访问数: 204
- HTML全文浏览量: 22
- PDF下载量: 72
