基于FLAC3D软件的抗滑桩设计推力曲线计算方法与流程

本发明涉及岩土工程数值模拟计算技术领域，特别涉及一种基于FLAC3D软件的抗滑桩设计推力曲线计算方法。

背景技术：

进行抗滑桩加固边坡的支挡结构设计时，抗滑桩的设计推力计算十分重要，抗滑桩设计推力的计算结果直接影响桩身截面尺寸、纵向配筋、锚固段长度、桩间距等设计指标，也会间接影响工程施工的难易程度、项目预算、施工周期等。基于工程设计目标，采用合理的计算理论，提高抗滑桩设计推力计算的合理性和精度，具有非常实用的工程应用价值和理论意义。

在我国滑坡灾害治理中，多采用传递系数法计算抗滑桩的设计推力，该方法采用刚体极限平衡理论，计算过程忽略边坡几何形态的影响，滑动面倾角变化能显著改变抗滑桩设计推力的计算结果。由于刚体极限平衡理论忽略土体的弹塑性本构特性，不能合理反映边坡土体的受力变形特性和破坏过程。针对传递系数法计算理论的局限性，有必要建立一种以弹塑性本构模型为基础的抗滑桩设计推力曲线计算方法，充分考虑土体的受力变形特性和边坡形态，弥补现有传递系数法计算理论的不足，具有重要的工程实践意义。

针对传递系数法计算边坡的抗滑桩设计推力曲线不能考虑土体的弹塑性本构特性、边坡形态，以及抗滑桩设计推力计算结果过于依赖滑动面倾角变化的特点，基于FLAC3D数值模拟计算软件，以土体的弹塑性本构模型为基础，建立横向宽度为1米的边坡三维模型，根据工程设计安全系数FS，按照公式φF＝tan^-1((tanφ)/FS)折减滑动面极限内摩擦角φ，按照公式cF＝c/FS折减滑动面极限粘聚力c，分别计算抗滑桩桩后土体的滑坡推力、弹性力，以及抗滑桩桩前土体抗力，桩后土体滑坡推力减去弹性力和桩前土体抗力的计算结果，即为监测截面的抗滑桩设计推力，将边坡各个监测截面的抗滑桩设计推力值连成曲线，建立基于FLAC3D软件的抗滑桩设计推力曲线计算方法。

技术实现要素：

本发明针对现有技术的缺陷，提供了一种基于FLAC3D软件的抗滑桩设计推力曲线计算方法，解决现有传递系数法忽略土体的弹塑性本构特性、不考虑边坡几何形态、滑动面倾角变化对滑坡推力影响过大的缺陷。

为了实现以上发明目的，本发明采取的技术方案如下：

一种基于FLAC3D软件的抗滑桩设计推力曲线计算方法，其特征在于，包括：横向宽度为1米的边坡三维模型、滑动面极限内摩擦角和极限粘聚力的计算方法、监测截面后侧土体产生弹性力的计算方法、监测截面后侧土体产生滑坡推力的计算方法、监测截面前侧土体抗力的计算方法和监测截面抗滑桩设计推力的计算方法；

根据所述横向宽度为1米的边坡三维模型，基于二分法的计算原理，利用LOOP循环语句折减滑动面内摩擦角和粘聚力的计算方法计算极限值，利用弹性本构模型计算监测截面后侧土体产生的弹性力，根据工程设计安全系数，利用摩尔-库仑本构模型计算滑动面极限内摩擦角和极限粘聚力折减时的桩后滑坡推力和桩前土体抗力，桩后滑坡推力减去弹性力和桩前土体抗力的计算结果，即为各个监测截面的抗滑桩设计推力，将它们连成曲线，得到边坡的抗滑桩设计推力曲线。

进一步地，所述横向宽度为1米的边坡三维模型具体为：边坡模型采用实体单元网格模拟，几何尺寸和原型保持一致，边坡模型的横向宽度均为1米，在边坡模型的前侧边界、后侧边界、底部边界和横向宽度方向上，设置固定速度和位移的边界条件。

进一步地，所述滑动面极限内摩擦角和极限粘聚力的计算方法，具体为：建立横向宽度为1米的边坡三维模型，以力不平衡比率不大于1×10^-5为收敛指标，将土体粘聚力和抗拉强度设为1×10¹⁰Pa，采用更改强度参数的摩尔-库仑本构模型生成初始地应力场。设置折减系数的上限值k2和下限值k1，设定折减系数的计算精度，并预定最大计算时步数。根据二分法的计算原理，采用FISH语言编写LOOP循环语句，参照公式和折减边坡滑动面的内摩擦角φ0和粘聚力c0，当力不平衡比率不大于1×10^-5时，则下一次计算中的下限值替换为当力不平衡比率大于1×10^-5时，则下一次计算中上限值替换为当经过数次折减计算，最终折减系数仍为k2时，则应设置更大的上限值；当经过数次折减计算，最终折减系数仍为k1时，则应设置更小的下限值。利用摩尔-库仑本构模型求解相应的地应力场，通过循环折减计算，当上下限值的差值小于计算精度时，计算结束，得到滑动面的极限内摩擦角和极限粘聚力。

进一步地，所述监测截面后侧土体产生弹性力的计算方法，具体为：监测截面沿坡脚方向的岩土体定义为刚体，利用弹性本构模型计算监测截面相邻土体单元网格的水平土压力，再乘以相应监测截面土体单元网格高度，求和汇总得到监测截面后侧土体产生的弹性力。

进一步地，所述监测截面后侧土体产生滑坡推力的计算方法，具体为：根据工程设计安全系数FS，按照公式φF＝tan^-1((tanφ)/FS)和cF＝c/FS，对滑动面的极限内摩擦角φ和极限粘聚力c进行折减，将监测截面沿坡脚方向的岩土体定义为刚体，利用摩尔-库仑本构模型计算监测截面相邻土体单元网格的水平土压力，再乘以相应监测截面土体单元网格高度，求和汇总得到监测截面后侧土体产生的滑坡推力。

进一步地，所述监测截面前侧土体抗力的计算方法，具体为：根据工程设计安全系数FS，按照公式φF＝tan^-1((tanφ)/FS)和cF＝c/FS，对滑动面的极限内摩擦角φ和极限粘聚力c进行折减，将监测截面沿坡顶方向的岩土体定义为刚体，利用摩尔-库仑本构模型计算监测截面相邻土体单元网格的水平土压力，再乘以相应监测截面土体单元网格高度，求和汇总得到监测截面前侧的土体抗力。

进一步地，所述监测截面抗滑桩设计推力的计算方法，具体为：抗滑桩设计推力＝滑坡推力-弹性力-土体抗力。

与现有技术相比，本发明的优点在于：基于FLAC3D数值模拟计算软件，利用摩尔-库仑塑性本构模型、弹性本构模型进行边坡的抗滑桩设计推力计算，能合理反映边坡滑坡推力的变化规律，提高计算精度。

附图说明

图1为本发明基于FLAC3D软件的抗滑桩设计推力曲线计算流程图；

图2为本发明横向宽度1米的边坡三维模型图；

图3为本发明滑动面极限内摩擦角和极限粘聚力的计算流程图；

图4为本发明边坡模型的竖向监测截面划分示意图；

图5为本发明计算监测截面后侧土体弹性力或滑坡推力的三维模型图；

图6为本发明计算监测截面前侧土体抗力的三维模型图；

图7为本发明边坡的抗滑桩设计推力曲线图；

图7a为工程设计安全系数FS＝1.1时的抗滑桩设计推力曲线图；

图7b为工程设计安全系数FS＝1.25时的抗滑桩设计推力曲线图。

具体实施方式

为使本发明的目的、技术方案及优点更加清楚明白，以下根据附图并列举实施例，对本发明做进一步详细说明。

一种基于FLAC3D软件的抗滑桩设计推力曲线计算方法，如图1所示，该方法包括：

步骤1：建立横向宽度1米的边坡三维模型

根据边坡原型几何尺寸，利用FLAC3D数值模拟软件建立实体单元网格的边坡三维模型，横向宽度均为1米，如图2所示，在边坡模型的前侧边界、后侧边界、底部边界、横向宽度方向均固定节点速度和位移。

步骤2：计算边坡滑动面的极限内摩擦角和极限粘聚力

如图3所示，建立横向宽度为1米的边坡三维模型，以力不平衡比率不大于1×10^-5为收敛指标，将土体粘聚力和抗拉强度设为1×10¹⁰Pa，采用更改强度参数的摩尔-库仑本构模型生成初始地应力场。设置滑动面内摩擦角和粘聚力的折减系数上限值k2(如4)和下限值k1(如0.5)，设定折减系数的计算精度为0.001，将计算时步的最大值设为10000步。根据二分法的计算原理，采用FISH语言编写LOOP循环语句，参照公式和折减边坡滑动面的内摩擦角φ0和粘聚力c0，当力不平衡比率不大于1×10^-5时，则下一次计算中的下限值替换为当力不平衡比率大于1×10^-5时，则下一次计算中上限值替换为当经过数次折减计算，最终折减系数仍为k2时，则应设置更大的上限值；当经过数次折减计算，最终折减系数仍为k1时，则应设置更小的下限值。利用摩尔-库仑本构模型求解相应的地应力场，通过循环折减计算，当上下限值的差值小于计算精度0.001时，计算结束，得到滑动面的极限内摩擦角φ和极限粘聚力c。

步骤3：划分边坡模型的监测截面

根据工程设计要求，沿水平方向在边坡模型内部划分一定数量的竖向监测截面，如图4所示，监测截面即抗滑桩的考察桩位。

步骤4：计算监测截面后侧土体产生的弹性力

监测截面沿坡脚方向的岩土体定义为刚体，如图5所示，利用弹性本构模型进行计算，提取监测截面相邻土体单元网格的水平土压力b，再乘以相应监测截面土体单元网格高度h，监测截面后侧土体产生的弹性力为(b1×h1+…+bn×hn)。

步骤5：计算监测截面后侧土体产生的滑坡推力

监测截面沿坡脚方向的岩土体定义为刚体，如图5所示，以力不平衡比率不大于1×10^-5为收敛指标，将土体粘聚力和抗拉强度设为1×10¹⁰Pa，利用更改强度参数的摩尔-库仑本构模型生成初始地应力场。根据工程设计安全系数FS，按照公式φF＝tan^-1((tanφ)/FS)和cF＝c/FS，对滑动面的极限内摩擦角φ和极限粘聚力c进行折减，利用摩尔-库仑本构模型进行计算，提取监测截面相邻土体单元网格的水平土压力a，再乘以相应监测截面土体单元网格高度h，监测截面后侧土体产生的滑坡推力为(a1×h1+…+an×hn)。

步骤6：计算监测截面前侧土体抗力

监测截面沿坡顶方向的岩土体定义为刚体，如图6所示，以力不平衡比率不大于1×10^-5为收敛指标，将土体粘聚力和抗拉强度设为1×10¹⁰Pa，利用更改强度参数的摩尔-库仑本构模型生成初始地应力场。根据工程设计安全系数FS，按照公式φF＝tan^-1((tanφ)/FS)和cF＝c/FS，对滑动面的极限内摩擦角φ和极限粘聚力c进行折减，利用摩尔-库仑本构模型进行计算，当力不平衡比率大于1×10^-5时，不计入抗滑桩前侧土体抗力，当力不平衡比率不大于1×10^-5时，提取监测截面相邻土体单元网格的水平土压力c，再乘以相应监测截面土体单元网格高度h，监测截面前侧的土体抗力为(c1×h1+…+cn×hn)。

步骤7：根据不同的工程设计安全系数，计算监测截面的抗滑桩设计推力，对于每一个监测截面，抗滑桩设计推力＝(a1×h1+…+an×hn)-(b1×h1+…+bn×hn)-(c1×h1+…+cn×hn)，绘制边坡的抗滑桩设计推力曲线，如图7所示。

基于FLAC3D软件的抗滑桩设计推力曲线计算方法的计算精度和边坡三维模型单元体网格划分的精细程度相关，两者呈正相关性。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的实施方法，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。