目前，研究人员在边坡的稳定性问题上通常视土体为各向同性介质，即岩土材料在各个方向上表现出相同的物理、力学性质. 然而，土颗粒受重力影响下，长轴趋于水平方向沉积和排列，这种分布会引起土的强度表现出不同程度的各向异性. 此外，在原位条件下，作用于土体上的应力通常也是不等向的，而且在开挖或筑堤过程中土体的主应力轴会发生偏转，这也会使其天然强度随破坏面的方向不同而变化. Koutsoftas和Ladd^[1]曾建议，对于二维边坡应采用3种不同类型的试验获得土的抗剪强度指标，即与坡顶、坡中、坡底位置对应，分别采用平面应变压缩试验(PSC)、直剪试验(DSS)和平面应变伸长试验(PSE). 在这些试验中，土试样的应力状态能够更加接近坡体不同位置处的实际受力情况，因而所得结果也更加真实准确. 然而，实际应用中大都按照土体固结方向的抗剪强度指标来评价边坡的稳定性. 如果忽略这种各向异性的影响，在某些情况下可能会使分析结果与实际产生较大偏差.

目前已有诸多学者在边坡稳定性分析中进行相关研究. Lo^[2]切取与沉积方向成不同倾角的London黏土试样进行不排水三轴压缩试验，发现其强度随试样倾角的变化与Casagrande公式^[3]比较吻合. 按照常规的极限平衡法研究了黏聚力的各向异性对边坡稳定性的影响. Chen等^[4]视边坡失稳时的滑动体为一个独立单元，将破坏面采用对数螺旋线表示. 在考虑滑面不同位置处土体强度的差异后，利用极限分析上限定理对黏性土边坡进行研究，并得出稳定系数的表达式. 此后，在边坡稳定分析中考虑强度各向异性的极限平衡法^[5-7]、极限分析上限法^[8-10]及有限元法^[11]大多是在Lo^[2]和Chen等^[4]的基础上展开.

对试验结果进行回归分析得出的Casagrande公式虽然考虑了各向异性的影响，但却不具有普遍适用性且不能从理论上给出严格的证明. 本文从微观角度出发，在不同主应力方向下借助材料组构张量，建立了各向异性强度参数表达式. 在此基础上，利用加拿大岩土工程软件GeoStudio中的SLOPE/W模块，分析强度各向异性对边坡稳定的影响.

1 各向异性强度参数公式的建立 1.1 组构张量

组构是用于描述材料颗粒及其相关孔隙空间分布的专门术语. 对于黏土而言，其颗粒既不是球形也不是椭球形，但每一个颗粒有一个较长的主轴是确定的，当该主轴相对沉积面成不同角度时，其组构也将表现出不同程度的各向异性. 组构张量 ${{F}}' $ 是二阶对称张量，可用3个材料主方向分量F₁₁、F₂₂、F₃₃表示. 假定 ${F_{11}} \geqslant {F_{22}} \geqslant {F_{33}}$ ，如果主方向与参考坐标系重合，则该组构张量可表示为

${{F}}' = \left[ {\begin{array}{*{20}{c}} {{F_{11}}}&0&0 \\ 0&{{F_{22}}}&0 \\ 0&0&{{F_{33}}} \end{array}} \right].$

(1)

一般认为土是横观各向同性的材料，所以组构张量中的两个主分量F₂₂和F₃₃相等. 由于组构张量 ${{F}}' $ 具有单位迹，也即F₁₁+F₂₂+F₃₃=1，所以对于轴对称横观各向同性土体，只需一个标量参数 $\varDelta $ 便可定义相应的组构张量. Oda等^[12]给出了如下形式

${{F}}' = \frac{1}{{3 + \varDelta }}\left[ {\begin{array}{*{20}{c}} {1 - \varDelta }&0&0 \\ 0&{1{\rm{ + }}\varDelta }&0 \\ 0&0&{1{\rm{ + }}\varDelta } \end{array}} \right].$

(2)

式(2)中， $\varDelta $ 代表颗粒优势排列的各向异性程度，依赖于颗粒形状和土的沉积状态；当 $\varDelta $ =0时，组构张量的主对角元素值相等，表示材料为各向同性；当 $\varDelta $ =1时，组构张量的主对角元素值一个为0，另外两个相等，表示颗粒的长轴均位于横观各向同性水平面内.

如果材料主方向与参考坐标系不重合，则该组构张量F的各个分量 ${F_{ij}}$ 可表示为

${F_{ij}} = {Q_{ki}}{Q_{lj}}F_{kl}'.$

(3)

式(3)中， ${Q_{ij}} = {{e}}_i' \cdot {{{e}}_j}$ ，为 ${{e}}_i'$ 和 ${{{e}}_j}$ 夹角的余弦. ${{e}}_i'$ 、 ${{{e}}_j}$ 分别为原坐标系和新坐标系中第i个、第j个基向量.

1.2 各向异性状态参数

由中主应力参数 $b = ({\sigma _2} - {\sigma _3})/({\sigma _1} - {\sigma _3})$ ，可得 ${\sigma _2} = b{\sigma _1} + (1 - b){\sigma _3}$ ；平均主应力 ${\sigma _m} = p =({\sigma _1} + {\sigma _2} + $ $ {\sigma _3})/3$ $ = [(1 + b){\sigma _1} + (2 - b){\sigma _3}]/3$ ；偏应力 ${s_{ij}} = {\sigma _{ij}} - $ $ {\sigma _m}{\delta _{ij}}$ ；归一化后的单位偏应力张量 ${{n}} = {{s}}/\left| s \right|$ . 经数学推导，可得单位偏应力张量n的各个分量，如下式

${{n}} = \frac{1}{{{{\left[ {6({b^2} - b + 1)} \right]}^{1/2}}}}\left[ {\begin{array}{*{20}{c}} {2 - b}&0&0 \\ 0&{2b - 1}&0 \\ 0&0&{ - (1 + b)} \end{array}} \right].$

(4)

由式(4)可知，单位偏应力张量n为二阶无迹张量，即tr n=0且有tr n²=1，因此单位偏应力张量n只表示偏应力的方向而不包括偏应力数值的大小. 根据Pietruszczak等^[13]的研究，可以定义一个组构张量在单位偏应力方向上的投影，记为 ${{F}}:{{n}}$ . 由于 ${\rm{tr}}\; {{n}} = 0$ ，所以 ${{F}}:{{n}}$ 最终只反映了组构张量的偏张量和单位偏应力张量间互相作用的结果，可以定义各向异性状态参数A的表达式为^[14]

$A = {{{F}}:{{n}}} = {{F}}^\pi :{{{n}}}.$

(5)

式中， $F_{ij}^\pi = {F_{ij}} - {F_{kk}}{\delta _{ij}}/3$ ， ${\delta _{ij}} $ 为Kronecker符号.

将式(2)~(4)代入式(5)中，化简后可得各向异性状态参数A具体为

$A = \sqrt {\frac{2}{3}} \frac{\varDelta }{{3 + \varDelta }}\frac{{1 + b - 3{{\cos }^2}\alpha }}{{{{({b^2} - b + 1)}^{1/2}}}} = D \times \frac{{1 + b - 3{{\cos }^2}\alpha }}{{{{({b^2} - b + 1)}^{1/2}}}}.$

(6)

式(6)中， $D = \sqrt {2/3} \varDelta /(3 + \varDelta )$ ；α为大主应力方向与沉积面法线间的夹角. 由此可见，各向异性状态参数A是标量参数 $\varDelta $ 、中主应力参数b和大主应力方向角α的函数.

1.3 各向异性强度参数公式

参考Yang等^[14]、Schweiger等^[15]和张玉军等^[16]的研究，对于Mohr-Coulomb准则，可以将土的抗剪强度指标c、φ看作为各向异性状态参数A的函数，在横观各向同性条件下，表达为

$c = {c_0}(1 + F_{ij}^cn_{ij}^c + F_{ijkl}^cn_{ijkl}^c + F_{ijklmn}^cn_{ijklmn}^c + \cdots ).$

(7)

$\varphi \!=\! \arctan [\tan {\varphi _0}(1 \!+\! F_{ij}^\varphi n_{ij}^\varphi \! + \!F_{ijkl}^\varphi n_{ijkl}^\varphi \!+\! F_{ijklmn}^\varphi n_{ijklmn}^\varphi \!+\! \cdots )].$

(8)

一般情况下，可写成 $F_{ijkl}^cn_{ijkl}^c = {d_1}F_{ij}^cF_{kl}^cn_{ij}^cn_{kl}^c$ ， $F_{ijklmn}^cn_{ijklmn}^c = {d_2}F_{ij}^cF_{kl}^cF_{mn}^cn_{ij}^cn_{kl}^cn_{mn}^c,\cdots$ ；同理可得 $F_{ijkl}^\varphi n_{ijkl}^\varphi = {b_1}F_{ij}^\varphi F_{kl}^\varphi n_{ij}^\varphi n_{kl}^\varphi $ ， $F_{ijklmn}^\varphi n_{ijklmn}^\varphi = {b_2}F_{ij}^\varphi F_{kl}^\varphi F_{mn}^\varphi \cdot$ $ n_{ij}^\varphi n_{kl}^\varphi n_{mn}^\varphi , \cdots$ . 式(7)和式(8)可进一步展开为

$c = {c_0}[1 + F_{ij}^cn_{ij}^c + {d_1}{(F_{ij}^cn_{ij}^c)^2} + {d_2}{(F_{ij}^cn_{ij}^c)^3} + \cdots ].$

(9)

$\varphi \! =\! \arctan \{ \tan {\varphi _0}[1 \!+\! F_{ij}^\varphi n_{ij}^\varphi \!+\! {b_1}{(F_{ij}^\varphi n_{ij}^\varphi )^2} \!+ \!{b_2}{(F_{ij}^\varphi n_{ij}^\varphi )^3} \!+\! \cdots ]\}. $

(10)

根据式(5)，可将式(9)与式(10)简化为

$\left\{\begin{array}{l}c = {c_0}(1 + {A_c} + {d_1}A_c^2 + {d_2}A_c^3 + \cdots );\\\varphi = \arctan [\tan {\varphi _0}(1 + {A_\varphi } + {b_1}A_\varphi ^2 + {b_2}A_\varphi ^3 + \cdots )].\end{array} \right.$

(11)

式(11)中，A_c和A_φ分别为黏聚力和内摩擦角的各向异性状态参数，可由式(6)分别确定；c₀、φ₀、d_i和b_j分别为待定的材料参数.

1.4 参数确定与验证

选取赵红华等^[17]对上海褐黄色粉质黏土的试验结果. 由于该试验采用与沉积面(水平面)成不同角度切取土样，然后进行直剪试验，经过角度转换后可得大主应力偏离沉积面法线 $\alpha $ 角时土的抗剪强度参数c、φ试验值，如图1和图2所示.

直剪试验属于平面应变状态的范畴，而平面应变状态下的罗德角 ${\theta _\sigma }$ 一般在–7°左右，可用于求解平面应变状态下的边坡稳定性问题^[18-20]. 现取 ${\theta _\sigma } = - 7.5^\circ $ ，由 $\tan {\theta _\sigma } = (2b - 1)/\sqrt 3 $ 可得，中主应力参数 $b = 0.386$ . 将式(11)中A_c和A_φ的幂次分别取至1、2、3、…，对图1、图2所示的抗剪强度参试验值数进行最小二乘法拟合，可以得出不同各向异性程度的抗剪强度参数表达式.

A_c和A_φ的幂次取至1时，可得

$\left\{ \begin{array}{l}c = 19.44(1 + {A_c});\\\varphi = \arctan [0.366(1 + {A_\varphi })].\end{array} \right.$

(12)

式(12)中， ${A_c} = - 0.072\;9 \times \displaystyle\frac{{1 + b - 3{{\cos }^2}\alpha }}{{{{({b^2} - b + 1)}^{1/2}}}}$ ， $ {A_\varphi } = - 0.038\;3 \times $ $\displaystyle\frac{{1 + b - 3{{\cos }^2}\alpha }}{{{{({b^2} - b + 1)}^{1/2}}}}$ ， $b = 0.386$ .

A_c和A_φ的幂次取至2时，可得

$\left\{\begin{array}{l}c = 16.77(1 + {A_c} + 46.385A_c^2);\\\varphi = \arctan [0.324(1 + {A_\varphi } + 546.75A_\varphi ^2)].\end{array} \right.$

(13)

式(13)中， ${A_c} = - 0.047\;2 \times \displaystyle\frac{{1 + b - 3{{\cos }^2}\alpha }}{{{{({b^2} - b + 1)}^{1/2}}}}$ ， ${A_\varphi } = - 0.012\;3 \times $ $ \displaystyle\frac{{1 + b - 3{{\cos }^2}\alpha }}{{{{({b^2} - b + 1)}^{1/2}}}}$ ， $b = 0.386$ .

A_c和A_φ的幂次取至3时，可得

$\left\{ \begin{array}{l}c = 16.52(1 + {A_c} + 5.012A_c^2 - 11.865A_c^3);\\\varphi = \arctan [0.323(1 + {A_\varphi } + 62.805A_\varphi ^2 - 181.126A_\varphi ^3)].\end{array} \right.$

(14)

式(14)中， ${A_c} = - 0.155\;6 \times \displaystyle\frac{{1 + b - 3{{\cos }^2}\alpha }}{{{{({b^2} - b + 1)}^{1/2}}}}$ ， $ {A_\varphi } = - 0.037\;2 \times$ $\displaystyle\frac{{1 + b - 3{{\cos }^2}\alpha }}{{{{({b^2} - b + 1)}^{1/2}}}}$ ， $b = 0.386$ .

式(12)~(14)对抗剪强度参数的拟合曲线与其试验值之间的对比如图1和图2所示，其中1~3次幂预测值分别与式(12)~(14)对应. 由图可见，在式(11)中随着各向异性状态参数A_c和A_φ最高幂次的增加，式(12)~(14)的计算结果对试验值的拟合效果越来越好. 当A_c和A_φ的幂次取至3时，误差最小，式(14)能够较好地反映大主应力方向对该黏土抗剪强度指标的影响，同时也验证了本文所建各向异性强度参数公式(11)的有效性.

2 边坡稳定性分析 2.1 典型的破坏模式

图3给出了边坡典型的圆弧破坏模式及潜在滑动面上土条单元的应力状态变化. 自然条件下，通常可认为土体的沉积面为水平面. 图中折线ABC为边坡坡面线，圆弧AC为潜在滑动面，θ_i为土条i的滑面底边倾角，α_i为土条i底面单元所受大主应力σ₁的方向与土层沉积面法线方向(竖直方向)的夹角，ψ为大主应力σ₁的方向与土条滑裂面切线的夹角，β为坡角，H为坡高，其他符号见图3.

2.2 算例分析与讨论 2.2.1 算例1

某一均质黏土边坡模型几何参数为：坡高H=20 m，坡比(竖直:水平)为1:1.5，即坡角β=33.7°. 土体抗剪强度参数选用赵红华等^[17]对上海褐黄色粉质黏土的试验结果，天然重度 $\gamma $ =19.4 kN/m³. 设边坡的潜在滑动面如图3所示，由图可知，采用条分法分析边坡稳定性时，大主应力σ₁偏离竖直方向的角度α_i = 90° – θ_i – ψ. 大主应力σ₁的方向与土条滑裂面切线的夹角一般由ψ = 45° – φ/2确定，根据Lo^[2]的试验研究，发现ψ通常在30°~40°之间. 而由图2中内摩擦角φ的变化，可得ψ的值在33.4°~36.1°之间，据此本算例计算时可取平均值ψ=35°，其对α_i的影响仅在1°左右，因而可以忽略不计.

各向同性条件下，取α=0°时的抗剪强度参数c=23.8 kPa，φ=23.2°，利用GeoStudio中的SLOPE/W模块通过Spencer法计算得到的安全系数F_sI=1.378；各向异性条件下，基于式(14)，利用GeoStudio中的SLOPE/W模块搜索出本例边坡模型的临界滑动面，并根据滑面上不同土条的底面倾角θ_i确定不同的α_i角后，代入不同土条的抗剪强度参数，通过Spencer法可求得边坡安全系数F_sA=1.158. 比较可知，考虑强度各向异性的影响后，本例边坡的安全系数降低约16.0%.

2.2.2 算例2

选取王栋等^[11]对某均质土坡的分析计算模型：坡高 ${{H}} = 8.12$ m，坡角 $\beta = 70^\circ $ ；土体参数为重度 $\gamma = $ 20 kN/m³，内摩擦角 $\varphi = 20^\circ $ 且不考虑其各向异性的影响，垂直方向的黏聚力c_v=20 kPa，各向异性黏聚力c采用Casagrande法表达且各向异性比K=0.7，结果如图4所示.

对于二维边坡稳定性问题，可认为各向异性黏聚力c为平面应变状态下的结果，与算例1的做法相同，可取中主应力参数 $b = 0.386$ . 对不同主应力方向下的黏聚力c进行最小二乘法拟合，可得

$c = 16.77(1 + {A_c}).$

(15)

式(15)中， ${A_c} = - 0.104\;2 \times \displaystyle\frac{{1 + b - 3{{\cos }^2}\alpha }}{{{{({b^2} - b + 1)}^{1/2}}}}$ ； $b = 0.386$ . 式(15)对Casagrande法计算所得黏聚力 $c$ 的拟合效果如图4所示，可见二者结果完全吻合.

与算例1类似，采用Spencer法时的计算结果如表1所示，相对差值为(各向同性值−各向异性值)/各向同性值. 各向同性情况下，取黏聚力c=20 kPa，内摩擦角φ=20°.

由表1中的结果可见，各向异性情况下，本文解法得到的边坡安全系数与文献[11]采用Casagrande法得到的结果差别较小. 各向同性与各向异性情况下的结果也基本相同，相对差值仅为1.34%. 这是因为对于坡角β=70°的这类陡坡，临界滑动面为靠近坡面的浅层滑动，使得条块底面倾角θ_i较大，大主应力σ₁偏离竖直方向的角度α_i = 90° – θ_i – ψ较小. 由图4可见，当α减小时，黏聚力c的各向异性程度减弱，最终使两种情况下的边坡安全系数差别较小.

3 坡角对各向异性边坡的影响

本文前两个算例分别代表了典型的缓坡与陡坡，初步得出如下结论：缓坡角时，各向同性与各向异性情况下的边坡安全系数差别较大；而陡坡角时，两种情况下的边坡安全系数差别较小. 为了不失一般性，本节将建立不同坡角β和坡高H的边坡模型，通过大量计算进一步验证该结论的真实性.

与图3所示的边坡几何形状相似，建立如下不同几何尺寸的边坡模型：坡高H分别取10 m、15 m、20 m；坡角β分别为14°、21.8°、26.6°、33.7°、45°、55°. 边坡土体的物理力学参数与算例1相同，各向同性情况下的抗剪强度参数c=23.8 kPa、φ=23.2°，各向异性情况下的抗剪强度参数采用式(14). 参照算例1的求解方式，利用Spencer法对上述36个边坡模型进行稳定分析，可得边坡安全系数F_s的计算结果如图5和图6所示.

由图5可见，在不同坡高下，各向同性与各向异性情况下的安全系数均随坡角的增加而减小；且随着坡角的增加，两种情况下的安全系数值越来越接近. 这种变化可由图6更加直观表现：对于坡角为14°的缓坡，考虑强度各向异性的影响后，边坡安全系数降低可达22.0%；随着坡角的增加，各向同性与各向异性安全系数相对差值逐渐减小；当坡角达到55°时，两种情况下的安全系数相对差值略大于8.0%.

4 结论

通过以上分析，本文得出的主要结论如下：

(1) 对于黏土材料，将各向异性状态参数定义为材料组构张量与标准化偏应力张量之间的联合不变量，提出了考虑抗剪强度指标 $c$ 和 $\varphi $ 各向异性的表达式，并结合已有的试验结果验证了该公式的有效性.

(2) 算例分析表明，不考虑土体强度的各向异性，可能会高估边坡的稳定性，尤其是对于高而平缓的边坡.

(3) 通过对36个边坡模型计算结果的分析，得出结论为：坡角较小时，各向同性与各向异性安全系数相对差值较大；随着坡角的增大，其相对差值逐渐降低，边坡受强度各向异性的影响减弱，即表明土体强度的各向异性对缓坡安全系数的影响比对陡坡的大. 实际应用中，应当重视强度各向异性对高而平缓边坡的影响.