基于Python的ABAQUS有限元强度折减法程序在边坡稳定性分析中的应用

荣光旭，彭 艳，田 凯

(1. 安徽工业经济职业技术学院 地质与建筑工程学院，安徽 合肥 230051；2. 中国地质调查局成都地质调查中心，四川 成都 610081)

0 引 言

边坡稳定性分析计算的目的是为了判断边坡稳定状态进而对发展趋势做出预测，以防止灾害的发生. 强度折减法是边坡稳定性分析的方法之一，它通过定义当前土体强度与避免结构破坏所需的最小抗剪强度之比来计算安全系数，或者可定义为即将发生破坏时，在Mohr-Coulomb破坏准则下的土体强度最小值.

有限元强度折减法(Strength Reduction Finite Element Method，SRFEM)理论首先由Zienkiewicz O.C[1]提出. SRFEM可应用于非均值材料，如土体破坏的研究[2]，但主要还是应用于边坡稳定性分析[3]. Yuan等[4]根据边坡土体内聚力折减参数和内摩擦角折减参数的比值不同，提出了基于双参数的折减方式，唐芬等[5]等讨论了黏性土和砂土c，φ不同的衰减速度，陈国庆等[6]研究了边坡强度折减的范围，Griffiths D.V等[7]提出了边坡失稳判据，曹先锋等[8]提出了基于温控参数的强度折减有限元法； 孙超伟等[9]利用强度折减法基于Hoek-Brown准则提出了一套边坡稳定性分析的图表法求解思路； 朱文炜等[10]综合分析了强度折减法与重度增加法对稳定性系数的影响； 李宁等[11-12]利用二分法来求解边坡的安全系数，同时利用Fortran编制程序并结合“基于场变量的有限元强度折减法”实现了ABAQUS的二次开发. Wei等[13]完成了单排桩加固的强度折减分析，Shi等[14]以等效塑性应变带连接作为破坏准则，利用SRFEM研究了钢管桩加固前后土坡的安全系数，Mohammad Reza Arvin等[15]利用SRFEM考虑了三维模型的情形下，对采用了土工格栅加固的边坡的安全系数进行了研究，同时对影响边坡稳定的主要因素进行了分析，Mehdipour I等[16]用极限平衡水平切片法分析了土工格栅加筋边坡的稳定性. 以上研究均从不同角度利用强度折减法对稳定性系数的求解过程进行了论述，但是关于如何在建模以及边坡稳定性分析过程中进行控制以利于计算过程的收敛，以及在边坡参数化建模分析等研究方面未做分析.

ABAQUS是一款功能强大的有限元模拟软件. 曹伟等、 张文东等、 秦宇等、 冷伍明等、 窦远明等[17-21]基于ABAQUS平台利用UMAT子程序及Python脚本分别完成了冻土蠕变模型、 材料裂纹模拟、 切削仿真模拟的二次开发，验证了Python脚本程序完成二次开发的可行性. 但是，利用Python编写脚本程序在ABAQUS中实现边坡稳定性分析的研究较少. 本文基于ABAQUS平台，利用具有高效率数据结构的Python语言进行二次开发，基于Mohr-Coulomb破坏准则，利用Python中字符串索引方式修改关键字，完整讨论了程序编写过程及要点，实现了基于场变量边坡稳定性系数的自动提交分析计算. 通过工程实例验证了程序的可行性.

1 有限元强度折减法

1.1 失稳判据

根据边坡失稳的定义，利用有限元强度折减法计算边坡稳定性系数有多种不同的解释. 本文采用常用的三种判据：

1) 塑性区从坡脚至坡顶贯通；

2) 数值求解不收敛；

3) 用户定义的某个节点出现大变形，即突变点，通常用于分析具体问题.

以上三种判据尽管存在不足[22-25]，但是经过不断改进[26-27]，在理论方面和实践方面均取得良好的效果. 塑性区贯通和位移突变点所得稳定性系数相差不大，据文献[27]可知，以位移突变标准作为失稳判据有其合理性，因此，本文以边坡特征点位移突变点作为失稳的评价标准.

1.2 强度折减理论

在外荷载不变的情况下，强度折减法认为边坡失稳主要是岩土体的剪切破坏. 在极限状况下，边坡内岩土体自身抗剪强度与岩土体抵御外荷载所发挥的最低抗剪强度的比值定义为边坡整体稳定性系数. 基本原理是将c，φ同时折减，折减后的抗剪强度参数可表示为

cm=c/Fr,

(1)

(2)

式中：c和φ是土体所能提供的抗剪强度；cm和φm是维持平衡所需要的或者土体实际发挥的抗剪强度；Fr是折减系数.计算过程中不断增加Fr，当边坡达到临界状态时，折减系数Fr就是边坡稳定性系数Fs.

2 二次开发思路

ABAQUS是岩土工程领域常用的有限元软件之一. Python作为面向对象的高级语言，亦是ABAQUS内核脚本语言. ABAQUS为用户提供了Python接口，用户可以利用Python绕过ABAQUS/CAE GUI，从而直接进行内核操作.

ABAQUS中未预置强度折减法，在进行边坡稳定性分析时，每次都需要将岩土体材料参数值及折减后的值代入进行试算，在满足前述失稳判据时得到最终边坡的稳定性系数值，效率较低. 对于几何形状相同但是特征参数变化的边坡稳定性分析，可以在ABAQUS中预先定义场变量，将土体材料黏聚力及内摩擦角随场变量折减，只需要一次提交即可完成计算过程，避免了重复仿真工作，实现了参数化分析，大大提高了工作效率.

当通过定义材料参数值随场变量变化的程序进行分析时，需要编辑修改关键字，在关键字中添加关于场变量的描述语句. 本文拟通过字符串索引方式，利用字符串修改添加分析步中场变量变化范围语句，真正实现自动提交分析. 本文讨论边坡稳定性分析中前处理及后处理模块中所有步骤的二次开发代码，采用图1 所示的结构框架.

图1 稳定性计算流程图

边坡稳定性分析二次开发基本思路是：前处理时建立集合Set，用于边界条件的设定； 在定义材料属性时设置场变量为1，并在0.5～2范围内变化； 修改关键字，将定义场变量语句添加到inp文件中； 计算完成后得到所需的数据及图件.

3 代码实现及应用实例

本文以二维边坡建模分析为例，所用环境如下：ABAQUS平台版本为6.13版本，Python为Anaconda自带的Python3.7，编译器采用PyCharm社区版. 因部分代码书写形式在ABAQUS6.13和6.14版本有所区别，且Python2.6和Python3.7版本语法形式亦有所不同，故此处写明代码适用版本. 以文献[28]中燕山集滑坡为例，该滑坡位于重庆市万州区燕山集镇，平面形状呈“长舌状”，高程405 m～406 m，前缘至燕山场镇沿江公路边，剪出口高程为328 m～329 m 左右. 纵向长305 m，横向宽约155 m，如图2 所示. 该滑坡相关参数见表1.

(a) 地理位置

(b) 现场示意图

表1 燕山集滑坡物理力学指标及几何参数

3.1 模块导入

导入abaqus，abaqusConstants，visualization模块以便脚本访问对象并可以使用对象成员和函数. 需要导入的模块较多，此处仅列举几个：

from abaqus import *

import regionToolset

此处可以使用changeKey更改模型名称，但需注意命名空间，代码如下：

myMdb.models.changeKey(fromName=“Model-1”,toName=“slope”)

slopeModel=mdb.models[“slope”]

3.2 建模

模型的建立需导入sketch和part模块，因为是二维模型，所以只需要依据实际几何模型建立点与点之间的连线后成为规则的图形，部分代码如下：

slopeSketch=slopeModel.ConstrainedSketch(name=“2D_slope_sketch”, sheetSize=10.0)

slopeSketch.Line(point1=(0,0),point2=(305.0,0))

slopeSketch.Line(point1=(305.0,0),point2=(305.0,16.0))

…

slopeSketch.Line(point1=(0,78),point2=(0,0))

为了作图方便，也可以通过ConstructionLine和FixedConstraint建立辅助线及约束.

3.3 材料

对材料命名并利用Material对象中的Density和Elastic对土体材料赋值，同时定义材料的c，φ随场变量的变化. 源码如下：

slopeMaterial.Denstity(table=((1.958,),))

slopeMaterialElastic(table=((100E6,0.30),))#弹性模量100 MPa，泊松比0.30

slopeMaterial.MohrCoulombPlasticity(dependencies=1, table=((18.208 4,0,0.5),(12.369 0,0,0.75),(9.34,0,1),(7.495 8,0,1.25),(6.257 4,0,1.5),(5.369 2,0,1.75),(4.701 2,0,2)))

slopeMaterial.MohrCoulombPlasticity.MohrCoulombHardening(dependencies=1,table=((19.62,0,0.5),(13.08,0,0.75),(9.81,0,1),(7.848,0,1.25),(6.54,0,1.5),(5.606,0,1.75),(4.905,0,2)))

后两句代码是定义变化范围0.5～2的场变量，并根据式(1)、 式(2)计算得出不同场变量下的c，φ值. 本部分代码在书写时需注意所有数据为元组(tuple)形式，当仅有一个数据时，后面一定要有逗号.

3.4 截面

通过HomogeneousSolidSection对截面命名，并指定材料名称为前述名称，对整个区域赋予截面属性，部分代码如下：

slope_region=(slope.faces,)

slopePart.SectionAssignment(region=slope_region,sectionName=“slopesection”,offset=0.0,offsetType=MIDDLE_SURFACE,offsetField=“”,thicknessAssignment=From_Section)

3.5 装配及分析步

利用根装配(rootAssembly)对象中的Instance方法将多个部件装配为一个实体，并命名，代码如下：

slopeInstance=slopeModel.rootAssembly.Instance(name=“slope instance”,part = slopePart,dependent = ON)

在定义分析步步骤中，定义两个分析步，分别命名为“load”，“reduce”； 增加输出结果场变量“FV”：

slopeModel.StaticStep(name=“load”,previous=“Initial”, initialInc=0.1, matrixSolver=DIRECT, matrixStorage = UNSYMMETRIC)#定义分析步“load”

slopeModel.StaticStep(name=“reduce”,previous=“load”,initialInc=0.1,matrixSolver=DIRECT,matrixStorage = UNSYMMETRIC)#定义分析步“load”

slopeModel.fieldOutputRequests.changeKey(fromName=“F-Output-1”,toName=“selected field outputs”)

slopeModel.fieldOutputRequests[“selected field outputs”].setValuesInStep(stepName=“reduce”,variables=(“S”,“PE”, “PEEQ”,“PEMAG”,“LE”,“U”,“RF”,“CF”,“CSTRESS”,“CDISP”,“FV”))#增加“FV”

ABAQUS中已经自定义initial分析步，本步骤中设置初始增量步为0.1，将Matrix storage选择为非对称分析(UNSYMMETRIC).

3.6 边界条件定义

在定义边界条件时，需要用到区域(regions), 它可以是集合、 表面对象或临时区域对象[29]. 本节使用findAt方法查找边(edge)，findAt方法的参数可以是边上的任意一点，返回值为包含该点ID的对象边，然后即可对边施加位移约束. 在边坡稳定性分析中，一般左、 右边界均限制水平位移，边坡底部限制水平和纵向两个方向的位移. 以下仅列举左边界限制水平位移的代码[30-31]：

left_edge=slopeModel.rootAssembly.instances[“slope instance”].edges.findAt(((0,40.0,0),))#查找左边界上的点

left_region=regionToolset.Region(edges=left_edge)#确定左边界

slopeModel.DisplacementBC(name=“leftBC-1”,createStepName=“reduce”,region=left_region,u1=0.0,u2=UNSET,ur3=UNSET, amplitude=UNSET,fixed=OFF,distributionType=UNIFORM,fieldName=“”,localCsys=None)#u1=0.0即为限制水平位移，底部边界为u1=0.0，u2=0.0

边坡的重力载荷通过施加体力的方式来实现. 对所有区域施加体力，需首先利用findAt方法确定面(face)，同样地，参数可以是面上任意一点，然后通过regionToolset赋给相应区域，施加体力. 代码如下：

Body_Force_faces=slopeModel.rootAssembly.instances[“slope instance”].faces.findAt(((120,43,0),))

Body_Force_region=regionToolset.Region(faces=Body_Force_faces)

slopeModel.BodyForce(name=“bodyforce”,createStepName=“reduce”,region=Body_Force_region,comp2=-20.0)

3.7 网格划分

首先利用3.6节所述的findAt方法确定面以后，赋予指定的区域，然后进行网格的划分. 本文中网格类型采用自由网格划分形式，单元为四节点平面应变减缩积分单元. 部分代码如下：

slope_mesh_elemType=mesh.ElemType(elemCode=CPE4R,elemLibrary=STANDARD)

slopeMeshFace=slopeInstance.faces.findAt(((120,35,0),))

slopeMeshFaceRegion=regionToolset.Region(faces=slopeMeshFace)

slopePart.setElementType(regions=slopeMeshFaceRegion,elemTypes=(slope_mesh_elemType,))

slopeModel.parts[“slope”].generateMesh()

燕山集滑坡二维有限元模型经过划分网格后，如图3 所示.

图3 划分网格后的模型

3.8 关键字修改及工作提交

ABAQUS中通过定义场变量实现强度折减法时需要定义“initial conditions”，因此，需要修改关键字(keywords)[32]. 本文通过文件处理的方式，利用Python中的字符串索引来确定需要修改的关键字的行数，然后插入替代的内容. 部分代码如下：

originstr = “**STEP:Reduce ”

newstr = “**STEP:Reduce Field,VARIABLE=1 slope-1.slope,2; ”#替代的关键字字符串

strindex = lines.index(originstr)#索引确定需要修改关键字的行数

lines[strindex]= newstr

lines.insert(strindex+1,newstr)#插入替代的内容

3.9 后处理

导入visualization模块，打开结果数据文件，查看稳定性系数及滑动面位置. 本算例在第二个分析步t=0.437 73时无法收敛，计算终止. 此时等效塑性应变云图如图4 所示，绘制场变量-位移曲线，即燕山集滑坡特征点的水平位移随折减系数变化的曲线，如图5 所示. 依据前述失稳判据理论，位移突变点所对应的FV1为0.99，即为燕山集滑坡的稳定性系数，这与剩余推力法所得的计算结果 1.012 基本一致. 从输出数据库读取数据部分代码如下：

from odbAcesss import *

from visualization import XYData

odb=visualization.openOdb(path=yanshanji.odb)#打开燕山集odb文件

xyData=session.XYData(plotName,data)

SafetyCurve=session.Curve(xyData)#绘制曲线

xyp=session.XYPlot(plotName)

chart = xyp.charts.values()[0]

chart.setValues(curvesToPlot=(SatyCurve, ), )

vp=session.viewports[“Viewport:1”]

vp.odbDisplay.setPrimaryVariable(variableLabel=“PEMAG”,outputPosition=INTEGRATION_POINT)#显示PEMAG云图

vp.odbDisplay.setPrimaryVariable(variableLabel=“U”)#显示计算终止时增量位移等值线云图，如图6 所示.

图4 t=0.437 73时的塑性区分布图

图5 特征点水平位移随FV变化的曲线

图6 t=0.437 73时的增量位移等值线云图

上述为获取常规结果数据输出的文件，一般边坡稳定性分析中只需要分析得到以上结果. 在某些特定分析中还可以利用getSubset获取边坡不同位置，如坡顶或坡中土体的位移场数据以判断边坡破坏类型，或者通过Viewport函数创建一个新的视图窗口，代码举例如下：

myOdb.step[“Step-1”].frames[1].fieldOutputs[“RF”].getSubset(region=myodb.rootAssembly.nodeSets[“setsname”]).values

myOdb=visualization.openOdb(path=jobName.job)

myViewport.setValues(displayedObject=myOdb)

myDisplacement=frame.fieldOutputs[“U”]

myStress=frame.fieldOutputs[“S”]

将上述代码保存为后缀名.py文件后从编译器退出. 运行该脚本程序的方法有多种，可以在ABAQUS/CAE中直接在菜单中选择从保存路径下Run Script.

4 结 论

本文以燕山集滑坡为例，完整讨论了利用Python语言进行有限元软件ABAQUS的二次开发，通过定义材料参数随场变量的变化以及利用字符串索引修改关键字的方式实现了边坡稳定性自动提交分析，使用方便，得到结果与极限平衡法基本一致. 主要结论有以下几点：

1) 基于场变量的边坡稳定性分析，避免了需多次修改基于cae产生的inp文件中的c和φ值，以及每次均需要提交工作才能获得最终的稳定性系数的问题，分析过程得以简化，达到参数化分析的目的.

2) 通过编写Python脚本程序，调用ABAQUS内核分析实现稳定性系数计算过程的自动化. 脚本程序文件便于研究人员更为深入地理解响应机理，建模以及分析过程中对模型的调试和诊断更有利于计算过程的收敛.

3) 本文通过Python语言中字符串索引方式修改关键字，顺利完成了模型输入文件的修改，实现了自动提交分析.