考虑土体各向异性的抗滑桩边坡稳定性分析

宁 帅 庄 妍 谈云志 王东星

（1.东南大学 苏州联合研究生院, 南京 211189；2.东南大学 土木工程学院, 南京 211189；3.三峡大学三峡库区地质灾害教育部重点实验室, 湖北 宜昌 443002；4.武汉大学 土木建筑工程学院, 武汉 430072）

关于土拱效应的研究一直受到许多学者的关注,1943年,Terzaghi(太沙基)[1]率先验证了土拱效应的真实存在,ITO等[2]通过塑性流动理论和塑性变形理论推导出土拱效应下桩的侧向力公式,Richard等[3]针对小主应力的作用轨迹,提出了土拱形状是一条悬链线的观点；Kellogg等[4]将土拱按大致的形状分为半球形、三角形、圆顶形、抛物线形等.由于土拱形状的不确定性和土体参数的复杂性,导致关于土拱效应的理论计算方法十分有限,并且大多数的计算方法都是以一定的假设为前提条件.通过数值模拟的方法,学者们在土拱效应的研究上取得了许多进展.Ladanyi等[5-10]利用理论与数值模拟相结合的方法,对土拱效应展开了更为详细的研究；董捷等[11]运用数值模拟对三维土拱效应进行了研究并将土拱分为水平、竖向和临空3类.目前比较常见的研究是通过数值模拟得到的数据与理论计算得出的结果进行对比验证,对现有的计算方法进行更准确的修正.同时,部分学者将对土拱的研究结果应用于室内试验与实际工程中,为工程的设计和安全提供了更加充分地保障[12-15].在抗滑桩加固边坡的稳定性研究中,已有许多学者对抗滑桩后土体的土拱效应进行了详细的研究[16-21],通过这些研究可以发现,土体几乎都被视为各向同性,并且简化了数值模拟试验的模型和加载过程,因此无法充分体现出真实情况下抗滑桩加固边坡的复杂性.在抗滑桩边坡的数值模拟试验中,极少有学者能够在有限元软件中考虑到各向异性来模拟实际工况.然而,就笔者调查研究所知,如果在试验过程中将土体视为各向同性,计算结果所得的安全系数会明显高于土体为各向异性时的安全系数,并且差距会随着土体的参数和模型的尺寸而改变.因此,在计算抗滑桩边坡稳定性时,需要将土体的各向异性条件考虑在内.王栋等[22]利用ABAQUS和强度折减法实现了各向异性条件下边坡安全系数的数值分析,但是其模拟的情况为未加固的普通边坡,极少有学者在有限元软件中充分考虑各向异性土体和水平土拱效应来研究抗滑桩加固边坡的稳定性.本文针对以上问题建立了抗滑桩加固边坡的力学模型,通过对ABAQUS子程序的编写及调用,根据迭代计算结果动态地调整了每一增量步上每一个积分点的黏聚力,从而实现了土体强度的各向异性,利用强度折减法求出抗滑桩加固各向异性边坡的安全系数,并通过改变模型的桩间距,分析了土拱效应和安全系数在不同桩间距下的变化情况,得到了考虑了各向异性和水平土拱效应条件下的最优桩间距.

1 计算模型及基本算例

1.1 基本计算模型

本文基本计算模型取自文献[23],为了充分发挥出抗滑桩的加固作用,令桩的埋深贯通土体,本文边坡的宽度取2倍的桩间距以减小边界效应的影响,具体参数如图1所示.

图1 基本模型计算尺寸

土体与抗滑桩的相关参数见表1.在建模的过程中,桩和土之间的切向接触为Penalty,摩擦系数0.3,约为tan(0.75φ) ,法向接触为Hard.整个模型采用底部固定、侧面法向的约束形式,其他平面不设置约束条件.

表1 抗滑桩-边坡体系强度参数

1.2 强度折减有限元法

在有限元软件中,强度折减法根据迭代计算不收敛或者滑动面上的位移突变作为边坡失稳依据,因此安全系数可视为使边坡达到失稳状态时的强度折减系数[24].通常根据下式进行计算：

式中：cm和φm为土体实际所需的抗剪强度；c和φ为土体所能够提供的抗剪强度；Fs为强度折减系数.

1.3 基本算例及其结果验证

基本算例中的土体设为均匀土体,未考虑各向异性,有限元模型网格如图2所示,为了增加计算效率并保证计算结果的精确度,网格在滑动面与桩孔位置处进行了加密处理.

计算结果为：Fs=1.50,与原题1.536基本一致,相对误差不足3%.微小的差别可能来自不同的网格划分方法、网格的密度、增量步长等,通过对基本模型的验证,说明了本文采用的建模过程是正确且可行的.

图2 有限元网格

2 基于土体各向异性完善计算模型

2.1 各向异性在ABAQUS中的实现

在基本算例的基础上,通过改变特定的土体参数,从而实现土体的各向异性.由于ABAQUS主程序只能将土体视为黏聚力与内摩擦角均为固定常数的各向同性土体,因此需要通过调用子程序实现土体强度和内摩擦角的变化.相对于黏聚力的各向异性,内摩擦角φ的各向异性产生的效果并不明显,因此本文主要针对黏聚力的各向异性展开分析.各向异性的黏聚力通常按照下式确定[25]：

式中：cv和ch分别为垂直方向和水平方向的黏聚力；i为大主应力方向与垂直方向的夹角；k为各向异性系数,本文取k=0.5.

ABAQUS的编程通常是通过Fortran语言进行子程序的编写而实现的,理想的各向异性需要通过UMAT子程序构建完整的摩尔-库伦本构模型,本构模型的建立需要耗费巨大的工作量,因此本文根据王栋[22]采用的方法,使用主程序自带的本构模型,通过USFLD子程序自定义不同的场变量对应的值(如图3所示).

图3 材料参数设置

将黏聚力按照大主应力方向分为若干小段,每一小段的黏聚力在该段内取均值,通过实时调用每个积分点上的应力分量计算大主应力方向,然后对大主应力方向大小的判断确定场变量的序号,从而确定出该点的黏聚力取值,实现动态更新每个积分点的黏聚力.

2.2 考虑土体各向异性条件下的安全系数

通过强度折减法求解出安全系数,与原文结果对比发现,考虑了土体强度各向异性后,安全系数Fs为1.406,较未考虑土体各向异性时(Fs=1.536)有明显的降低,未考虑土体各向异性时的计算结果会高估抗滑桩边坡的稳定性.这是因为土体在考虑了各向异性条件时,黏聚力不再是某一固定的常数,而是根据i的不断变化而改变,i的取值范围为0°～90°,根据公式(5)可以发现,ci始终小于等于cv,因此当黏聚力整体呈现下降的趋势时,土体的抗滑能力便会下降,滑坡的风险随之上升,从而降低了安全系数.

3 桩间距对水平土拱效应以及安全系数的影响

3.1 桩间距对土拱范围的影响

在土体各向异性的基础上,对桩间距分别为2、4、8、12倍桩径时的土拱范围进行对比分析.由图4可知,桩间土体的主应力方向由于桩内侧与土体的摩擦而产生偏转,在桩间和桩前也呈现出一定范围的水平土拱效应,本文主要针对桩后的水平土拱效应研究桩间距的影响.如图4所示,当桩间距逐渐增加时,桩后拱轴线的曲率逐渐降低,当s/d=12时,桩后水平土拱已经消失,说明随着桩间距的增大,土拱越来越难以形成,当s/d=8以后,虽然能清晰地看到土拱的存在,但是其应力变化幅度已经非常小,可以理解为当s/d=8时,土拱效应已经非常不明显.

图4 水平土拱的变化与桩间距的关系

3.2 桩间距对中轴线上y方向的应力σy的影响

由图5可知：当s/d≤6时,中轴线上的σy的突变较为明显,x在-2～2 m的范围内出现了两个较为明显的拱形,当s/d≥8时,中轴线上的σy为一条平缓的曲线,两侧没有出现明显的拱效应.因此,根据水平土拱效应的发挥程度可以确定出极限桩间距(使桩后能够产生土拱效应的最大桩间距),由此可知模型的极限桩间距为s=8d.

图5 中轴线上y方向应力分量与桩间距的关系

3.3 桩间距对y剖面上应力分量σx的影响

不同桩间距下的y剖面上应力分量σx的变化对比分析如图6所示：离中轴线越近,应力分量越小,当桩间距在较小的范围内(s/d≤3)时,水平土拱范围也较小,并且桩前的水平土拱效应较桩后更加明显；当桩间距逐渐增加时,桩后的土拱范围逐渐增大至饱和的拱形；继续增加桩间距,土拱范围开始逐渐缩小.

图6 桩前、后不同剖面上x方向应力分量

从应力变化曲线可以明显地观察到,当s/d=8时,拱形已经变形,从圆滑的曲线变成扁平的抛物线形状,因此极限桩间距为8倍的桩径.

3.4 桩间距对安全系数的影响

图7 安全系数与桩间距的关系图

图8 最优桩间距的确定

图7给出了不同桩间距下的安全系数曲线,如图7所示：当桩间距逐渐增加时,安全系数呈现出逐渐降低的趋势；当桩间距在较小的范围内时,安全系数降低的速度较快,超过这一范围时,继续增加桩间距,安全系数降低的速度逐渐趋向于零.如果将曲线分为两段,并将这两段曲线各自拟合成一条虚直线,如图8所示：虚线1的安全系数变化较大,而桩间距变化较小,因此虚线1以考虑安全系数为主；虚线2的安全系数变化较小,而桩间距变化较大,因此虚线2以考虑桩间距为主.两条虚线的交点处既充分考虑了桩间距的影响,又充分保证了稳定性,故取两条虚线交点的横坐标为最优桩间距,约为s/d=5.

4 结 论

1)通过对ABAQUS的二次开发实现了土体的各向异性,证明了各向异性对抗滑桩加固边坡有着不可忽视的影响,关于各向异性系数与初始黏聚力的大小对结果的影响还需进一步考究.

2)运用ABAQUS可以合理的模拟三维情况下土拱的形成过程,研究发现水平土拱的影响范围不一定是随着桩间距的增大而减小,在桩间距较小时,水平土拱的影响范围随着桩间距的增大而增大.

3)通过改变桩间距研究了不同路径下的应力变化情况,分析了不同桩间距下安全系数的变化,较为准确地计算出极限桩间距以及最优桩间距的大小.

4)运用强度折减法和二分法可以较为准确地计算出土体在各向异性条件下的安全系数,但是手动输入每一折减系数需要耗费较多的时间,因此在后续的研究中,安全系数的求解方法还需进一步优化.