宁帅， 庄妍,2， 崔晓艳， 张永攀

(1.东南大学 土木工程学院，江苏 南京 211189；2.东南大学 混凝土及预应力混凝土结构教育部重点实验室，江苏 南京 211189；3.河海大学 岩土力学与堤坝工程教育部重点实验室，江苏 南京 210098)

抗滑桩施工工艺简单且对边坡的扰动较小，因此抗滑桩在加固边坡中得到了广泛应用[1-4]。抗滑桩边坡稳定性分析大多将土体视为强度各向异性和均质性，但是大部分自然沉积土的抗剪强度具有各向异性与非均质性[5-7]。目前为止，强度各向异性条件下的边坡稳定性主要集中在：1) 极限分析理论，Chen等[8]提出了极限分析上限解；Hwang等[9]通过有限元法得出边坡不同位置处水平与竖向应力之比，给出了开挖边坡内平面应变压区和拉区的分布；Zdravkovic等[10]通过有限元软件模拟了足尺路堤破坏试验，研究了土体强度各向异性对软土路堤稳定性的影响;2) 极限平衡理论，AL等[11]采用传统的极限平衡法探讨了二维均质粘性土坡中粘结力的各向异性对土坡稳定安全系数Fs和临界滑动面位置的影响；王建锋等[12]运用普遍条带分法对各向异性边坡进行分析，建立一个能够同时耦合各向异性和非线性的新的强度准则。有限元分析软件虽然解决了手动计算造成的误差大、过程繁琐等缺点，但是模型的建立和优化需要耗费大量的时间调试，一个参数的改变往往需要通过重新建模来实现。国内外的规范大多以极限平衡理论为主，英国加筋路堤设计规范《BS 8006-1:2010+A1:2016》[13](简称英国规范BS8006)考虑了抗滑桩的竖向承载力Fpi和加筋层Tr的抗拉作用，中国规范《GB 50843-2013》[14]与《GB 50330-2013》[15]均采用沿滑动面切线方向的力的平衡求解安全系数，文献[14-15]额外考虑了支护结构的有效抗力R0i，规范《JTG/T D31-02-2013》[16]采用对圆弧滑动面圆心的力矩平衡求解安全系数，其在抗滑力矩中考虑了加筋层的抗拉作用TGCj。英国规范《BS 8006-2-2011》[17]考虑了泥钉轴向设计力Tdj的分解值；另外，美国规范《ANSI/ASCE 1-82·1983》[18]虽然未给出完整的计算式，但是其将边坡稳定性研究的安全因素定义为基于可用强度与外加应力或其他荷载效应的比率，仍是以极限平衡的方法求解安全系数；欧洲规范《EN 1997-1-2004》[19]给出了计算过程中不同荷载分项系数组合；美国规范《ASCE 20-96·1997》[20]在计算边坡稳定安全系数时采用经验系数法等规范方法同样值得参考应用。英国加筋路堤设计规范BS8006[13]可以用于计算抗滑桩边坡的安全系数，但并没有考虑土体强度各向异性和非均质性。

本文将强度各向异性和非均质性引入BS8006，得到修正英国规范BS8006抗滑桩边坡安全系数计算方法。运用ABAQUS有限元软件建立考虑强度各向异性和非均质性的抗滑桩边坡有限元模型，验证计算方法的适用性。

1 考虑土体强度各向异性、非均质性的抗滑桩加固边坡稳定计算方法

英国规范BS8006中的抗滑桩加固路堤的稳定分析方法基于极限平衡法，具有计算思路明确合理、计算过程简便、适用性强等特点。因此，本文在英国规范BS8006方法的基础上，将强度各向异性和非均质性引入边坡稳定计算方法中，同时将边坡中的水平土拱效应考虑在内，对边坡土体以及抗滑桩的抗滑力矩进行修正，以提出一种考虑强度各向异性和非均质性抗滑桩边坡安全系数的计算方法。

1.1 英国路堤整体稳定性设计规范BS8006

英国规范BS8006方法计算简图如图1所示，抗滑桩加固路堤的安全系数为：

图1 计算路堤稳定安全系数简化模型

(1)

式中：MD为土体的下滑力矩；MRS为土的抗滑力矩；MRP为桩的抗滑力矩；MRR为加筋体的抗滑力矩。

土体下滑力矩为：

(2)

式中：Wi为土条自重；Wsi是土条外部荷载；bi是土条宽度；ffs为重力荷载分项系数，取ffs=1.3；fq为外部荷载分项系数，取fq=1.3。

土体抗滑力矩：

(3)

式中：c为土体粘聚力；ru为与孔隙水压力相关的系数；fms为土体参数分项系数，对于tanφ，分项系数取fms=1.0，对于c，分项系数取fms=1.6。

桩体抗滑力矩：

MRP=FP1XP1+FP2XP2

(4)

式中：FP1、FP2为滑动面内桩体承担的竖向承载力；XP1、XP2为桩体到滑动面圆心的水平距离。

加筋体抗滑力矩为：

MRR=TrY

(5)

式中：Tr为加筋体的最大拉力;Y为最大拉力距滑动面圆心的垂直距离。

1.2 考虑土体强度的各向异性

边坡稳定性分析中，由于固结应力，历史应力以及沉积物类型的不同使得土体表现出明显的各向异性现象，如图2所示。如果采用Mohr-Coulomb理想弹塑性模型描述土的应力-应变关系，内摩擦角φ通常并不表现出显著的各向异性，因此本文假定粘聚力c具有各向异性和非均质性，各向异性的粘聚力为[11,21-22]：

图2 开挖边坡中的强度各向异性

ci=ch+(cv-ch)cos2i

(6)

式中：ch为水平方向上的粘聚力强度；cv为土体竖向粘聚力强度；i为最大主应力与竖直方向的夹角。各向异性系数可以表示为水平与竖向黏聚力强度的比值，即：

(7)

式中粘聚力ci与k和i的关系如图3所示，本文假定在同一种土体材料中各向异性系数相同。将各向异性系数引入式(6)可得：

图3 粘聚力各向异性

(8)

1.3 考虑土体的非均质性

土体非均质性表现在土体粘聚力随着深度增加而增大，土体在正常固结条件下，粘聚力c随着深度的变化一般分为5种模式[23-24]。模式1假设土体的粘聚力为常量，即土体为均质的，其余模式假设土体的粘聚力沿深度而发生变化；模式2以圆弧滑动面为计算模型，将土体视为正常固结不排水，地下水位在地面以下时粘聚力c随着深度的变化，故选取模式3为本文的非均质模型，粘聚力c在边坡中的具体分布如图4所示。图4中ch0为土体层面初始粘聚力，h为土体高度，λ为粘聚力线性变化时的斜率。粘聚力ch为：

图4 粘聚力随深度的变化情况

ch=ch0+λh

(9)

1.4 土体抗滑力矩修正

本文应用极限平衡方法研究强度各向异性非均质边坡的稳定性，修正英国规范BS8006边坡稳定安全系数计算方法。假定强度各向异性非均质边坡圆弧滑动面如图5所示。

图5 强度各向异性非均质边坡的圆弧滑动面

该边坡滑动面半径为R，坡顶滑动面与水平向圆心夹角为θ0，坡顶处与圆心高差为H0，圆心夹角为θi的滑动面处大主应力与竖向夹角为i，滑动面与大主应力面夹角为m，滑动面与水平向夹角为αi。由图5中几何关系可以得到：

H0=R·sinθ0

(10)

(11)

(12)

剪破面与大主应力面的夹角m可以通过莫尔-库伦破坏准则得到：

(13)

结合式(9)和(10)可以得到该点水平粘聚力为：

ch=ch0+λR(sinθi-sinθ0)

(14)

将式(14)代入式(8)可以得到强度各向异性非均质边坡在该点处的粘聚力ci为：

ci=kich0+λRki(sinθi-sinθ0)

(15)

依据BS8006土体抗滑力矩MRS式(3)在各向异性非均质边坡中就求解抗滑力矩时需对式中粘聚力分项进行修正。将第i块土条粘聚力分项定义为Zi,则在三维边坡内有：

Zi=cibi

(16)

式中：bi为第i块土条宽度。为了计算方便将土条所在滑动面弧长近似为土条宽度bi。在条块内将Zi沿滑动面进行积分则有：

Csin3θ+Dsin(2θ)sinθ+λRsin2θ+Ecos2θsin2θ+

Fsin4θ+Gsin(2θ)sin2θ-Esinθ0cos2θsinθ-

Fsinθ0sin3θ-Gsinθ0sin(2θ)sinθ)dθ

(17)

其中：

对式(17)进行积分计算可得：

(18)

其中：

当只考虑土体各向异性时：

(19)

当只考虑土体非均质性时：

(20)

由式(18)最终得到应用英国规范在各向异性非均质边坡中抗滑力矩式为：

(21)

1.5 抗滑桩抗滑力矩修正

Ito[25]基于Tomio塑性变形理论以及土拱效应，得到单位厚度土层作用在桩上的侧向力，并得到文献[26-28]的应用和验证。侧向力为：

(22)

x轴方向上单位厚度土层作用在桩上的侧向力可以表示为：

(23)

式中：Pbc是抗滑桩在土体表层所受的力；Phdm是抗滑桩在滑动面层所受的力；h为抗滑桩所在位置处滑动面深度；l为梯形形心处到滑动面圆形处的垂直距离。

对于只考虑强度各向异性土体边坡，水平向粘聚力ch不沿着深度发生变化。Ito[25]计算抗滑桩极限侧向压力的式采用水平粘聚力ch进行计算，因而不考虑cv对抗滑桩极限侧向压力的影响。求出抗滑桩极限侧向压力ph后，只需考虑抗滑桩所在位置滑动面深度即可得出考虑土体强度各向异性条件下抗滑桩的抗滑力矩MRP。侧向压力ph和抗滑力矩MRP分别为：

ph=pBB′-D2{σx}x=0=

(24)

(25)

对于强度各向异性的非均质土体边坡，由于水平粘聚力ch随深度的变化为线性变化，且式中与ch相乘的参数在已知工况下均为常数，因此文献[25]计算抗滑桩极限侧向压力的式同样适用。

2 各向异性边坡有限元数值模型

由于ABAQUS主程序默认土体为各项同性，本文采用文献[29]实现各向异性在ABAQUS中的应用，步骤为：

1) 编写子程序，调取每个积分点上的应力值，计算出各点的i值，并将i划分成若干区段，每一区段的i对应不同的场变量：

(26)

2) 在主程序中添加子程序，结合式(6)计算出每一区段上的平均粘聚力，将每一区段上的平均粘聚力分别存放于对应的场变量中，通过调用子程序选取相应场变量下的粘聚力值，从而实时更新每个积分点上的粘聚力；

3) 为了实现非均质性，将土体水平分成若干层，如图6所示。根据式(9)计算出每层土粘聚力的上下限并取平均值，每一层按照步骤1)、2)实现各向异性；

图6 土体非均质性的抗滑桩加固边坡简化模型

4) 由于ABAQUS无法直接对各向异性和非均质性土体进行强度折减，因此本文利用二分法对各向异性和非均质性土体进行手动折减并得到安全系数。

本文有限元数值模型取自文献[30]，为验证建立的抗滑桩边坡安全系数计算方法的正确性，将计算结果与抗滑桩加固软土边坡的数值模型进行对比分析，并采用该模型分析强度各向异性和非均质性对抗滑桩边坡稳定性的影响。材料具体参数和模型尺寸分别由表1和图7所示。

表1 有限元模型材料参数

图7 有限元模型具体尺寸

在建模的过程中，材料选取摩尔-库伦弹塑性本构模型，为减小边界效应，模型宽度取2倍桩间距，桩径为1 m，桩体打入至模型底部，桩和土之间的切向接触为罚函数，摩擦系数0.3，约为tan(0.75φ)，法向接触为硬接触，模型底部为固定约束，侧面采用法向约束形式。为近似实现非均质性，本文根据不同的分层方案以及结果对比，考虑网格划分难易和计算效率影响，选取层高均为1 m的分层方案。模型的网格划分如图8所示，为了提高计算精度与效率，在斜坡位置与桩孔处对网格进行了加密处理，桩与土的单元类型均为C3D8。

图8 有限元模型网格划分

3 有限元与解析解对比分析

当只考虑各向异性条件时(非均质系数n=1)，三轴和直剪试验中观察到0.5≤k≤1.0[21]，根据修正BS8006抗滑桩边坡稳定性计算方法和数值模拟计算，得到了不同各向异性系数下均质边坡安全系数，由表2所示。

表2 不同各向异性系数下的安全系数

根据各向异性系数k=0.5和k=1时分析可得，未考虑各向异性条件时，抗滑桩边坡的安全系数将明显高于考虑各向异性条件下时抗滑桩边坡的安全系数，约高出9.25%。

其次，在考虑各向异性条件时(k=0.5)，引入非均质性条件，根据修正BS8006抗滑桩边坡稳定性计算方法和数值模拟计算，分别得到了不同非均质性系数n下各向异性边坡安全系数，计算结果由表3所示。

表3 不同非均质性系数下的安全系数

根据各向异性系数n=0.5、k=0.5和n=1、k=1时分析可得，当抗滑桩边坡未考虑各向异性和非均质性条件时，抗滑桩边坡的安全系数将显著提高，约提高21.42%。实际工程中的边坡通常都有不同程度的强度各向异性和非均质性，但是目前评价抗滑桩边坡稳定性时大多按各向同性处理和均质性处理，此时如果按照强度各向同性和均质性去计算边坡的稳定性，就会严重高估边坡的稳定性。因此，各向异性与非均质性在抗滑桩边坡稳定性计算中起着不可忽视的作用，实际计算中应当将各向异性和非均质性考虑在内。

4 结论

1)通过对现行规范的对比研究，发现英国规范BS8006中抗滑桩路堤稳定性的计算方法计算过程简单、较好的模拟了抗滑桩边坡的工况。通过对BS8006中平衡力矩的修正，加入抗滑桩边坡水平土拱效应的影响，可以实现考虑土体强度各向异性与非均质性的抗滑桩边坡稳定性计算方法。

2)强度各向异性均质边坡的安全系数较各项同性均质边坡降低约9.25%；强度各向异性非均质边坡的安全系数较各项同性均质边坡计算所得安全系数降低了20%以上，因此，土体强度各向异性和非均质性在抗滑桩边坡稳定性分析中应予以考虑。

3)建立了粘聚力强度各向异性和非均质性的抗滑桩边坡有限元模型，通过有限元分析软件的二次开发对修正BS8006抗滑桩边坡稳定性计算方法进行对比验证，为后续进行各向异性的抗滑桩边坡稳定性影响因素分析奠定基础。