基于Python的ABAQUS有限元強(qiáng)度折減法程序在邊坡穩(wěn)定性分析中的應(yīng)用

榮光旭，彭 艷，田 凱

(1. 安徽工業(yè)經(jīng)濟(jì)職業(yè)技術(shù)學(xué)院 地質(zhì)與建筑工程學(xué)院，安徽 合肥 230051；2. 中國地質(zhì)調(diào)查局成都地質(zhì)調(diào)查中心，四川 成都 610081)

0 引 言

邊坡穩(wěn)定性分析計(jì)算的目的是為了判斷邊坡穩(wěn)定狀態(tài)進(jìn)而對(duì)發(fā)展趨勢(shì)做出預(yù)測(cè)，以防止災(zāi)害的發(fā)生. 強(qiáng)度折減法是邊坡穩(wěn)定性分析的方法之一，它通過定義當(dāng)前土體強(qiáng)度與避免結(jié)構(gòu)破壞所需的最小抗剪強(qiáng)度之比來計(jì)算安全系數(shù)，或者可定義為即將發(fā)生破壞時(shí)，在Mohr-Coulomb破壞準(zhǔn)則下的土體強(qiáng)度最小值.

有限元強(qiáng)度折減法(Strength Reduction Finite Element Method，SRFEM)理論首先由Zienkiewicz O.C[1]提出. SRFEM可應(yīng)用于非均值材料，如土體破壞的研究[2]，但主要還是應(yīng)用于邊坡穩(wěn)定性分析[3]. Yuan等[4]根據(jù)邊坡土體內(nèi)聚力折減參數(shù)和內(nèi)摩擦角折減參數(shù)的比值不同，提出了基于雙參數(shù)的折減方式，唐芬等[5]等討論了黏性土和砂土c，φ不同的衰減速度，陳國慶等[6]研究了邊坡強(qiáng)度折減的范圍，Griffiths D.V等[7]提出了邊坡失穩(wěn)判據(jù)，曹先鋒等[8]提出了基于溫控參數(shù)的強(qiáng)度折減有限元法； 孫超偉等[9]利用強(qiáng)度折減法基于Hoek-Brown準(zhǔn)則提出了一套邊坡穩(wěn)定性分析的圖表法求解思路； 朱文煒等[10]綜合分析了強(qiáng)度折減法與重度增加法對(duì)穩(wěn)定性系數(shù)的影響； 李寧等[11-12]利用二分法來求解邊坡的安全系數(shù)，同時(shí)利用Fortran編制程序并結(jié)合“基于場(chǎng)變量的有限元強(qiáng)度折減法”實(shí)現(xiàn)了ABAQUS的二次開發(fā). Wei等[13]完成了單排樁加固的強(qiáng)度折減分析，Shi等[14]以等效塑性應(yīng)變帶連接作為破壞準(zhǔn)則，利用SRFEM研究了鋼管樁加固前后土坡的安全系數(shù)，Mohammad Reza Arvin等[15]利用SRFEM考慮了三維模型的情形下，對(duì)采用了土工格柵加固的邊坡的安全系數(shù)進(jìn)行了研究，同時(shí)對(duì)影響邊坡穩(wěn)定的主要因素進(jìn)行了分析，Mehdipour I等[16]用極限平衡水平切片法分析了土工格柵加筋邊坡的穩(wěn)定性. 以上研究均從不同角度利用強(qiáng)度折減法對(duì)穩(wěn)定性系數(shù)的求解過程進(jìn)行了論述，但是關(guān)于如何在建模以及邊坡穩(wěn)定性分析過程中進(jìn)行控制以利于計(jì)算過程的收斂，以及在邊坡參數(shù)化建模分析等研究方面未做分析.

ABAQUS是一款功能強(qiáng)大的有限元模擬軟件. 曹偉等、 張文東等、 秦宇等、 冷伍明等、 竇遠(yuǎn)明等[17-21]基于ABAQUS平臺(tái)利用UMAT子程序及Python腳本分別完成了凍土蠕變模型、 材料裂紋模擬、 切削仿真模擬的二次開發(fā)，驗(yàn)證了Python腳本程序完成二次開發(fā)的可行性. 但是，利用Python編寫腳本程序在ABAQUS中實(shí)現(xiàn)邊坡穩(wěn)定性分析的研究較少. 本文基于ABAQUS平臺(tái)，利用具有高效率數(shù)據(jù)結(jié)構(gòu)的Python語言進(jìn)行二次開發(fā)，基于Mohr-Coulomb破壞準(zhǔn)則，利用Python中字符串索引方式修改關(guān)鍵字，完整討論了程序編寫過程及要點(diǎn)，實(shí)現(xiàn)了基于場(chǎng)變量邊坡穩(wěn)定性系數(shù)的自動(dòng)提交分析計(jì)算. 通過工程實(shí)例驗(yàn)證了程序的可行性.

1 有限元強(qiáng)度折減法

1.1 失穩(wěn)判據(jù)

根據(jù)邊坡失穩(wěn)的定義，利用有限元強(qiáng)度折減法計(jì)算邊坡穩(wěn)定性系數(shù)有多種不同的解釋. 本文采用常用的三種判據(jù)：

1) 塑性區(qū)從坡腳至坡頂貫通；

2) 數(shù)值求解不收斂；

3) 用戶定義的某個(gè)節(jié)點(diǎn)出現(xiàn)大變形，即突變點(diǎn)，通常用于分析具體問題.

以上三種判據(jù)盡管存在不足[22-25]，但是經(jīng)過不斷改進(jìn)[26-27]，在理論方面和實(shí)踐方面均取得良好的效果. 塑性區(qū)貫通和位移突變點(diǎn)所得穩(wěn)定性系數(shù)相差不大，據(jù)文獻(xiàn)[27]可知，以位移突變標(biāo)準(zhǔn)作為失穩(wěn)判據(jù)有其合理性，因此，本文以邊坡特征點(diǎn)位移突變點(diǎn)作為失穩(wěn)的評(píng)價(jià)標(biāo)準(zhǔn).

1.2 強(qiáng)度折減理論

在外荷載不變的情況下，強(qiáng)度折減法認(rèn)為邊坡失穩(wěn)主要是巖土體的剪切破壞. 在極限狀況下，邊坡內(nèi)巖土體自身抗剪強(qiáng)度與巖土體抵御外荷載所發(fā)揮的最低抗剪強(qiáng)度的比值定義為邊坡整體穩(wěn)定性系數(shù). 基本原理是將c，φ同時(shí)折減，折減后的抗剪強(qiáng)度參數(shù)可表示為

cm=c/Fr,

(1)

(2)

式中：c和φ是土體所能提供的抗剪強(qiáng)度；cm和φm是維持平衡所需要的或者土體實(shí)際發(fā)揮的抗剪強(qiáng)度；Fr是折減系數(shù).計(jì)算過程中不斷增加Fr，當(dāng)邊坡達(dá)到臨界狀態(tài)時(shí)，折減系數(shù)Fr就是邊坡穩(wěn)定性系數(shù)Fs.

2 二次開發(fā)思路

ABAQUS是巖土工程領(lǐng)域常用的有限元軟件之一. Python作為面向?qū)ο蟮母呒?jí)語言，亦是ABAQUS內(nèi)核腳本語言. ABAQUS為用戶提供了Python接口，用戶可以利用Python繞過ABAQUS/CAE GUI，從而直接進(jìn)行內(nèi)核操作.

ABAQUS中未預(yù)置強(qiáng)度折減法，在進(jìn)行邊坡穩(wěn)定性分析時(shí)，每次都需要將巖土體材料參數(shù)值及折減后的值代入進(jìn)行試算，在滿足前述失穩(wěn)判據(jù)時(shí)得到最終邊坡的穩(wěn)定性系數(shù)值，效率較低. 對(duì)于幾何形狀相同但是特征參數(shù)變化的邊坡穩(wěn)定性分析，可以在ABAQUS中預(yù)先定義場(chǎng)變量，將土體材料黏聚力及內(nèi)摩擦角隨場(chǎng)變量折減，只需要一次提交即可完成計(jì)算過程，避免了重復(fù)仿真工作，實(shí)現(xiàn)了參數(shù)化分析，大大提高了工作效率.

當(dāng)通過定義材料參數(shù)值隨場(chǎng)變量變化的程序進(jìn)行分析時(shí)，需要編輯修改關(guān)鍵字，在關(guān)鍵字中添加關(guān)于場(chǎng)變量的描述語句. 本文擬通過字符串索引方式，利用字符串修改添加分析步中場(chǎng)變量變化范圍語句，真正實(shí)現(xiàn)自動(dòng)提交分析. 本文討論邊坡穩(wěn)定性分析中前處理及后處理模塊中所有步驟的二次開發(fā)代碼，采用圖1 所示的結(jié)構(gòu)框架.

圖1 穩(wěn)定性計(jì)算流程圖

邊坡穩(wěn)定性分析二次開發(fā)基本思路是：前處理時(shí)建立集合Set，用于邊界條件的設(shè)定； 在定義材料屬性時(shí)設(shè)置場(chǎng)變量為1，并在0.5～2范圍內(nèi)變化； 修改關(guān)鍵字，將定義場(chǎng)變量語句添加到inp文件中； 計(jì)算完成后得到所需的數(shù)據(jù)及圖件.

3 代碼實(shí)現(xiàn)及應(yīng)用實(shí)例

本文以二維邊坡建模分析為例，所用環(huán)境如下：ABAQUS平臺(tái)版本為6.13版本，Python為Anaconda自帶的Python3.7，編譯器采用PyCharm社區(qū)版. 因部分代碼書寫形式在ABAQUS6.13和6.14版本有所區(qū)別，且Python2.6和Python3.7版本語法形式亦有所不同，故此處寫明代碼適用版本. 以文獻(xiàn)[28]中燕山集滑坡為例，該滑坡位于重慶市萬州區(qū)燕山集鎮(zhèn)，平面形狀呈“長舌狀”，高程405 m～406 m，前緣至燕山場(chǎng)鎮(zhèn)沿江公路邊，剪出口高程為328 m～329 m 左右. 縱向長305 m，橫向?qū)捈s155 m，如圖2 所示. 該滑坡相關(guān)參數(shù)見表1.

(a) 地理位置

(b) 現(xiàn)場(chǎng)示意圖

表1 燕山集滑坡物理力學(xué)指標(biāo)及幾何參數(shù)

3.1 模塊導(dǎo)入

導(dǎo)入abaqus，abaqusConstants，visualization模塊以便腳本訪問對(duì)象并可以使用對(duì)象成員和函數(shù). 需要導(dǎo)入的模塊較多，此處僅列舉幾個(gè)：

from abaqus import *

import regionToolset

此處可以使用changeKey更改模型名稱，但需注意命名空間，代碼如下：

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

slopeModel=mdb.models[“slope”]

3.2 建模

模型的建立需導(dǎo)入sketch和part模塊，因?yàn)槭嵌S模型，所以只需要依據(jù)實(shí)際幾何模型建立點(diǎn)與點(diǎn)之間的連線后成為規(guī)則的圖形，部分代碼如下：

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 材料

對(duì)材料命名并利用Material對(duì)象中的Density和Elastic對(duì)土體材料賦值，同時(shí)定義材料的c，φ隨場(chǎng)變量的變化. 源碼如下：

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的場(chǎng)變量，并根據(jù)式(1)、 式(2)計(jì)算得出不同場(chǎng)變量下的c，φ值. 本部分代碼在書寫時(shí)需注意所有數(shù)據(jù)為元組(tuple)形式，當(dāng)僅有一個(gè)數(shù)據(jù)時(shí)，后面一定要有逗號(hào).

3.4 截面

通過HomogeneousSolidSection對(duì)截面命名，并指定材料名稱為前述名稱，對(duì)整個(gè)區(qū)域賦予截面屬性，部分代碼如下：

slope_region=(slope.faces,)

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

3.5 裝配及分析步

利用根裝配(rootAssembly)對(duì)象中的Instance方法將多個(gè)部件裝配為一個(gè)實(shí)體，并命名，代碼如下：

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

在定義分析步步驟中，定義兩個(gè)分析步，分別命名為“l(fā)oad”，“reduce”； 增加輸出結(jié)果場(chǎng)變量“FV”：

slopeModel.StaticStep(name=“l(fā)oad”,previous=“Initial”, initialInc=0.1, matrixSolver=DIRECT, matrixStorage = UNSYMMETRIC)#定義分析步“l(fā)oad”

slopeModel.StaticStep(name=“reduce”,previous=“l(fā)oad”,initialInc=0.1,matrixSolver=DIRECT,matrixStorage = UNSYMMETRIC)#定義分析步“l(fā)oad”

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中已經(jīng)自定義initial分析步，本步驟中設(shè)置初始增量步為0.1，將Matrix storage選擇為非對(duì)稱分析(UNSYMMETRIC).

3.6 邊界條件定義

在定義邊界條件時(shí)，需要用到區(qū)域(regions), 它可以是集合、 表面對(duì)象或臨時(shí)區(qū)域?qū)ο骩29]. 本節(jié)使用findAt方法查找邊(edge)，findAt方法的參數(shù)可以是邊上的任意一點(diǎn)，返回值為包含該點(diǎn)ID的對(duì)象邊，然后即可對(duì)邊施加位移約束. 在邊坡穩(wěn)定性分析中，一般左、 右邊界均限制水平位移，邊坡底部限制水平和縱向兩個(gè)方向的位移. 以下僅列舉左邊界限制水平位移的代碼[30-31]：

left_edge=slopeModel.rootAssembly.instances[“slope instance”].edges.findAt(((0,40.0,0),))#查找左邊界上的點(diǎn)

left_region=regionToolset.Region(edges=left_edge)#確定左邊界

slopeModel.DisplacementBC(name=“l(fā)eftBC-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

邊坡的重力載荷通過施加體力的方式來實(shí)現(xiàn). 對(duì)所有區(qū)域施加體力，需首先利用findAt方法確定面(face)，同樣地，參數(shù)可以是面上任意一點(diǎn)，然后通過regionToolset賦給相應(yīng)區(qū)域，施加體力. 代碼如下：

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 網(wǎng)格劃分

首先利用3.6節(jié)所述的findAt方法確定面以后，賦予指定的區(qū)域，然后進(jìn)行網(wǎng)格的劃分. 本文中網(wǎng)格類型采用自由網(wǎng)格劃分形式，單元為四節(jié)點(diǎn)平面應(yīng)變減縮積分單元. 部分代碼如下：

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()

燕山集滑坡二維有限元模型經(jīng)過劃分網(wǎng)格后，如圖3 所示.

圖3 劃分網(wǎng)格后的模型

3.8 關(guān)鍵字修改及工作提交

ABAQUS中通過定義場(chǎng)變量實(shí)現(xiàn)強(qiáng)度折減法時(shí)需要定義“initial conditions”，因此，需要修改關(guān)鍵字(keywords)[32]. 本文通過文件處理的方式，利用Python中的字符串索引來確定需要修改的關(guān)鍵字的行數(shù)，然后插入替代的內(nèi)容. 部分代碼如下：

originstr = “**STEP:Reduce ”

newstr = “**STEP:Reduce Field,VARIABLE=1 slope-1.slope,2; ”#替代的關(guān)鍵字字符串

strindex = lines.index(originstr)#索引確定需要修改關(guān)鍵字的行數(shù)

lines[strindex]= newstr

lines.insert(strindex+1,newstr)#插入替代的內(nèi)容

3.9 后處理

導(dǎo)入visualization模塊，打開結(jié)果數(shù)據(jù)文件，查看穩(wěn)定性系數(shù)及滑動(dòng)面位置. 本算例在第二個(gè)分析步t=0.437 73時(shí)無法收斂，計(jì)算終止. 此時(shí)等效塑性應(yīng)變?cè)茍D如圖4 所示，繪制場(chǎng)變量-位移曲線，即燕山集滑坡特征點(diǎn)的水平位移隨折減系數(shù)變化的曲線，如圖5 所示. 依據(jù)前述失穩(wěn)判據(jù)理論，位移突變點(diǎn)所對(duì)應(yīng)的FV1為0.99，即為燕山集滑坡的穩(wěn)定性系數(shù)，這與剩余推力法所得的計(jì)算結(jié)果 1.012 基本一致. 從輸出數(shù)據(jù)庫讀取數(shù)據(jù)部分代碼如下：

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”)#顯示計(jì)算終止時(shí)增量位移等值線云圖，如圖6 所示.

圖4 t=0.437 73時(shí)的塑性區(qū)分布圖

圖5 特征點(diǎn)水平位移隨FV變化的曲線

圖6 t=0.437 73時(shí)的增量位移等值線云圖

上述為獲取常規(guī)結(jié)果數(shù)據(jù)輸出的文件，一般邊坡穩(wěn)定性分析中只需要分析得到以上結(jié)果. 在某些特定分析中還可以利用getSubset獲取邊坡不同位置，如坡頂或坡中土體的位移場(chǎng)數(shù)據(jù)以判斷邊坡破壞類型，或者通過Viewport函數(shù)創(chuàng)建一個(gè)新的視圖窗口，代碼舉例如下：

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文件后從編譯器退出. 運(yùn)行該腳本程序的方法有多種，可以在ABAQUS/CAE中直接在菜單中選擇從保存路徑下Run Script.

4 結(jié) 論

本文以燕山集滑坡為例，完整討論了利用Python語言進(jìn)行有限元軟件ABAQUS的二次開發(fā)，通過定義材料參數(shù)隨場(chǎng)變量的變化以及利用字符串索引修改關(guān)鍵字的方式實(shí)現(xiàn)了邊坡穩(wěn)定性自動(dòng)提交分析，使用方便，得到結(jié)果與極限平衡法基本一致. 主要結(jié)論有以下幾點(diǎn)：

1) 基于場(chǎng)變量的邊坡穩(wěn)定性分析，避免了需多次修改基于cae產(chǎn)生的inp文件中的c和φ值，以及每次均需要提交工作才能獲得最終的穩(wěn)定性系數(shù)的問題，分析過程得以簡(jiǎn)化，達(dá)到參數(shù)化分析的目的.

2) 通過編寫Python腳本程序，調(diào)用ABAQUS內(nèi)核分析實(shí)現(xiàn)穩(wěn)定性系數(shù)計(jì)算過程的自動(dòng)化. 腳本程序文件便于研究人員更為深入地理解響應(yīng)機(jī)理，建模以及分析過程中對(duì)模型的調(diào)試和診斷更有利于計(jì)算過程的收斂.

3) 本文通過Python語言中字符串索引方式修改關(guān)鍵字，順利完成了模型輸入文件的修改，實(shí)現(xiàn)了自動(dòng)提交分析.