利用ABAQUS中的rpy文件实现边坡可靠度计算

1.rpy文件

rpy（repeat python）：在Abaqus/CAE中建立模型时，工作目录下将自动生成abaqus.rpy文件。该文件记录了所有Abaqus/CAE操作对应的命令，将其复制，改变扩展名，由.rpy改为.py即可对脚本进行编辑。

利用rpy文件进行脚本编辑有以下两个优点：
1、可避免编写较长的命令。
2、可减少代码出错概率。

2.边坡可靠度计算

对于abaqus中的边坡模型输入不同的参数（c、φ）后进行安全系数(F)的计算。利用蒙特卡洛法原理，使得频率逼近概率，得到边坡的失效概率。

3.实操

实现思路：
1、首先在Abaqus中进行边坡模型的建立，并输出inp文件。（该步骤前文已经详细介绍过，在此不赘述）。
2、在abaqus工作路径中找到rpy文件（该文件中包含刚才建模过程中的所有操作对应的python代码），将其复制出来，并修改其后缀名为.py。当然，此文件中包含许多冗余的不影响建模的语句，需自己识别后删除即可。
3、利用循环语句，将不同的参数输入模型，再输出inp文件：
①利用随机数生成输入参数c、phi值（注意：abaqus中c与phi值必须为正值，所以随机数必须为正值）。
②在循环语句内，将生成的随机数输入相对位置，完成参数赋值。
③在循环语句内，执行inp输出语句，即可输出不同c、phi值的模型inp文件。
4、批量计算、批量提取odb文件（该步骤前文已经详细介绍过，在此不赘述）。

下面为建模后rpy文件的内容：

# -*- coding: mbcs -*-
#
# Abaqus/CAE Release 2020 replay file
# Internal Version: 2019_09_14-01.49.31 163176
# Run by pc on Sun Jan 14 10:58:07 2024
## from driverUtils import executeOnCaeGraphicsStartup
# executeOnCaeGraphicsStartup()
#: Executing "onCaeGraphicsStartup()" in the site directory ...
from abaqus import *
from abaqusConstants import *
session.Viewport(name='Viewport: 1', origin=(0.0, 0.0), width=192.67707824707, height=155.555557250977)
session.viewports['Viewport: 1'].makeCurrent()
session.viewports['Viewport: 1'].maximize()
from caeModules import *
from driverUtils import executeOnCaeStartup
executeOnCaeStartup()
session.viewports['Viewport: 1'].partDisplay.geometryOptions.setValues(referenceRepresentation=ON)
s = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=200.0)
g, v, d, c = s.geometry, s.vertices, s.dimensions, s.constraints
#以下代码为模型创建代码
s.setPrimaryObject(option=STANDALONE)
s.Line(point1=(0.0, 0.0), point2=(20.0, 0.0))
s.HorizontalConstraint(entity=g[2], addUndoState=False)
s.Line(point1=(20.0, 0.0), point2=(20.0, 13.0))
s.VerticalConstraint(entity=g[3], addUndoState=False)
s.PerpendicularConstraint(entity1=g[2], entity2=g[3], addUndoState=False)
s.Line(point1=(20.0, 13.0), point2=(12.0, 13.0))
s.HorizontalConstraint(entity=g[4], addUndoState=False)
s.PerpendicularConstraint(entity1=g[3], entity2=g[4], addUndoState=False)
s.Line(point1=(12.0, 13.0), point2=(2.0, 3.0))
s.Line(point1=(2.0, 3.0), point2=(0.0, 3.0))
s.HorizontalConstraint(entity=g[6], addUndoState=False)
s.Line(point1=(0.0, 3.0), point2=(0.0, 0.0))
s.VerticalConstraint(entity=g[7], addUndoState=False)
s.PerpendicularConstraint(entity1=g[6], entity2=g[7], addUndoState=False)
p = mdb.models['Model-1'].Part(name='slope', dimensionality=TWO_D_PLANAR, type=DEFORMABLE_BODY)
p = mdb.models['Model-1'].parts['slope']
p.BaseShell(sketch=s)
s.unsetPrimaryObject()
p = mdb.models['Model-1'].parts['slope']
session.viewports['Viewport: 1'].setValues(displayedObject=p)
del mdb.models['Model-1'].sketches['__profile__']
p = mdb.models['Model-1'].parts['slope']
f, e, d1 = p.faces, p.edges, p.datums
t = p.MakeSketchTransform(sketchPlane=f[0], sketchPlaneSide=SIDE1, origin=(12.175439, 5.508772, 0.0))
s1 = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=47.7, gridSpacing=1.19, transform=t)
g, v, d, c = s1.geometry, s1.vertices, s1.dimensions, s1.constraints
s1.setPrimaryObject(option=SUPERIMPOSE)
p = mdb.models['Model-1'].parts['slope']
p.projectReferencesOntoSketch(sketch=s1, filter=COPLANAR_EDGES)
s1.Line(point1=(-10.175439, -2.508772), point2=(7.82456100005157, -2.508772))
s1.HorizontalConstraint(entity=g[8], addUndoState=False)
s1.CoincidentConstraint(entity1=v[6], entity2=g[4], addUndoState=False)
s1.Line(point1=(-10.175439, -2.508772), point2=(-10.175439, -5.50877199998597))
s1.VerticalConstraint(entity=g[9], addUndoState=False)
s1.CoincidentConstraint(entity1=v[7], entity2=g[3], addUndoState=False)
session.viewports['Viewport: 1'].view.setValues(width=19.74, height=14.774, cameraPosition=(10.1437, 6.22066, 47.7074), cameraTarget=(10.1437, 6.22066, 0))
p = mdb.models['Model-1'].parts['slope']
f = p.faces
pickedFaces = f.getSequenceFromMask(mask=('[#1 ]', ), )
e1, d2 = p.edges, p.datums
p.PartitionFaceBySketch(faces=pickedFaces, sketch=s1)
s1.unsetPrimaryObject()
del mdb.models['Model-1'].sketches['__profile__']
session.viewports['Viewport: 1'].partDisplay.setValues(sectionAssignments=ON, engineeringFeatures=ON)
session.viewports['Viewport: 1'].partDisplay.geometryOptions.setValues(referenceRepresentation=OFF)
mdb.models['Model-1'].Material(name='soil')
mdb.models['Model-1'].materials['soil'].Elastic(table=((100000.0, 0.33), ))
mdb.models['Model-1'].materials['soil'].MohrCoulombPlasticity(dependencies=1, table=((36.0, 0.0, 0.5), (25.887, 0.0, 0.75), (20.0, 0.0, 1.0), (16.2343, 0.0, 1.25), (13.639, 0.0, 1.5), (11.749, 0.0, 1.75), (10.3141, 0.0, 2.0)))
mdb.models['Model-1'].materials['soil'].mohrCoulombPlasticity.MohrCoulombHardening(dependencies=1, table=((24.76, 0.0, 0.5), (16.5067, 0.0, 0.75), (12.38, 0.0, 1.0), (9.904, 0.0, 1.25), (8.25333, 0.0, 1.5), (7.07429, 0.0, 1.75), (6.19, 0.0, 2.0)))
mdb.models['Model-1'].materials['soil'].mohrCoulombPlasticity.TensionCutOff(temperatureDependency=OFF, dependencies=0, table=((0.0, 0.0), ))
mdb.models['Model-1'].HomogeneousSolidSection(name='Section-1', material='soil', thickness=None)
session.viewports['Viewport: 1'].view.setValues(nearPlane=41.8992, farPlane=53.5157, width=28.4409, height=20.0339, viewOffsetX=0.702957, viewOffsetY=0.352482)
p = mdb.models['Model-1'].parts['slope']
f = p.faces
faces = f.getSequenceFromMask(mask=('[#7 ]', ), )
region = regionToolset.Region(faces=faces)
p = mdb.models['Model-1'].parts['slope']
p.SectionAssignment(region=region, sectionName='Section-1', offset=0.0, offsetType=MIDDLE_SURFACE, offsetField='', thicknessAssignment=FROM_SECTION)
a = mdb.models['Model-1'].rootAssembly
session.viewports['Viewport: 1'].setValues(displayedObject=a)
session.viewports['Viewport: 1'].assemblyDisplay.setValues(optimizationTasks=OFF, geometricRestrictions=OFF, stopConditions=OFF)
a = mdb.models['Model-1'].rootAssembly
a.DatumCsysByDefault(CARTESIAN)
p = mdb.models['Model-1'].parts['slope']
a.Instance(name='slope-1', part=p, dependent=ON)
session.viewports['Viewport: 1'].assemblyDisplay.setValues(adaptiveMeshConstraints=ON)
mdb.models['Model-1'].StaticStep(name='load', previous='Initial', initialInc=0.1, matrixSolver=DIRECT, matrixStorage=UNSYMMETRIC)
session.viewports['Viewport: 1'].assemblyDisplay.setValues(step='load')
mdb.models['Model-1'].StaticStep(name='reduce', previous='load', initialInc=0.1, matrixSolver=DIRECT, matrixStorage=UNSYMMETRIC)
session.viewports['Viewport: 1'].assemblyDisplay.setValues(step='reduce')
session.viewports['Viewport: 1'].assemblyDisplay.setValues(step='load')
mdb.models['Model-1'].fieldOutputRequests['F-Output-1'].setValues(variables=('S', 'PE', 'PEEQ', 'PEMAG', 'LE', 'U', 'RF', 'CF', 'CSTRESS', 'CDISP', 'FV'))
del mdb.models['Model-1'].fieldOutputRequests['F-Output-1']
mdb.models['Model-1'].FieldOutputRequest(name='F-Output-1', createStepName='load', variables=('FV', ))
session.viewports['Viewport: 1'].assemblyDisplay.setValues(interactions=ON, constraints=ON, connectors=ON, engineeringFeatures=ON, adaptiveMeshConstraints=OFF)
session.viewports['Viewport: 1'].assemblyDisplay.setValues(loads=ON, bcs=ON, predefinedFields=ON, interactions=OFF, constraints=OFF, engineeringFeatures=OFF)
session.viewports['Viewport: 1'].view.setValues(nearPlane=43.0546, farPlane=52.3603, width=22.8176, height=16.0728, viewOffsetX=0.249526, viewOffsetY=0.464813)
a = mdb.models['Model-1'].rootAssembly
f1 = a.instances['slope-1'].faces
faces1 = f1.getSequenceFromMask(mask=('[#7 ]', ), )
region = regionToolset.Region(faces=faces1)
mdb.models['Model-1'].BodyForce(name='Load-1', createStepName='load', region=region, comp2=-20.0)
session.viewports['Viewport: 1'].assemblyDisplay.setValues(step='Initial')
a = mdb.models['Model-1'].rootAssembly
e1 = a.instances['slope-1'].edges
edges1 = e1.getSequenceFromMask(mask=('[#c4 ]', ), )
region = regionToolset.Region(edges=edges1)
mdb.models['Model-1'].DisplacementBC(name='BC-1', createStepName='Initial', region=region, u1=SET, u2=UNSET, ur3=UNSET, amplitude=UNSET, distributionType=UNIFORM, fieldName='', localCsys=None)
a = mdb.models['Model-1'].rootAssembly
e1 = a.instances['slope-1'].edges
edges1 = e1.getSequenceFromMask(mask=('[#28 ]', ), )
region = regionToolset.Region(edges=edges1)
mdb.models['Model-1'].DisplacementBC(name='BC-2', createStepName='Initial', region=region, u1=SET, u2=SET, ur3=UNSET, amplitude=UNSET, distributionType=UNIFORM, fieldName='', localCsys=None)
session.viewports['Viewport: 1'].assemblyDisplay.setValues(mesh=ON, loads=OFF, bcs=OFF, predefinedFields=OFF, connectors=OFF)
session.viewports['Viewport: 1'].assemblyDisplay.meshOptions.setValues(meshTechnique=ON)
p = mdb.models['Model-1'].parts['slope']
session.viewports['Viewport: 1'].setValues(displayedObject=p)
session.viewports['Viewport: 1'].partDisplay.setValues(sectionAssignments=OFF, engineeringFeatures=OFF, mesh=ON)
session.viewports['Viewport: 1'].partDisplay.meshOptions.setValues(meshTechnique=ON)
session.viewports['Viewport: 1'].view.setValues(width=26.5023, height=18.7857, viewOffsetX=0.771327, viewOffsetY=0.424364)
p = mdb.models['Model-1'].parts['slope']
f = p.faces
pickedRegions = f.getSequenceFromMask(mask=('[#7 ]', ), )
p.setMeshControls(regions=pickedRegions, elemShape=QUAD, technique=SWEEP)
elemType1 = mesh.ElemType(elemCode=CPE4, elemLibrary=STANDARD)
elemType2 = mesh.ElemType(elemCode=CPE3, elemLibrary=STANDARD)
p = mdb.models['Model-1'].parts['slope']
f = p.faces
faces = f.getSequenceFromMask(mask=('[#7 ]', ), )
pickedRegions =(faces, )
p.setElementType(regions=pickedRegions, elemTypes=(elemType1, elemType2))
p = mdb.models['Model-1'].parts['slope']
p.seedPart(size=1.0, deviationFactor=0.1, minSizeFactor=0.1)
p = mdb.models['Model-1'].parts['slope']
p.seedPart(size=0.5, deviationFactor=0.1, minSizeFactor=0.1)
p = mdb.models['Model-1'].parts['slope']
p.generateMesh()
session.viewports['Viewport: 1'].partDisplay.setValues(mesh=OFF)
session.viewports['Viewport: 1'].partDisplay.meshOptions.setValues(meshTechnique=OFF)
session.viewports['Viewport: 1'].partDisplay.geometryOptions.setValues(referenceRepresentation=ON)
p = mdb.models['Model-1'].parts['slope']
f = p.faces
faces = f.getSequenceFromMask(mask=('[#7 ]', ), )
p.Set(faces=faces, name='slope')
#: The set 'slope' has been created (3 faces).
import job
#一下代码为修改关键字代码
mdb.models['Model-1'].keywordBlock.synchVersions(storeNodesAndElements=False)
mdb.models['Model-1'].keywordBlock.insert(33, """
*initial conditions, type=field, variable=1
slope-1.slope,0.5""")
mdb.models['Model-1'].keywordBlock.insert(49, """*field,variable=1
slope-1.slope,2""")
a = mdb.models['Model-1'].rootAssembly
session.viewports['Viewport: 1'].setValues(displayedObject=a)
session.viewports['Viewport: 1'].assemblyDisplay.setValues(mesh=OFF)
session.viewports['Viewport: 1'].assemblyDisplay.meshOptions.setValues(meshTechnique=OFF)
mdb.Job(name='Job-1', model='Model-1', description='', type=ANALYSIS, atTime=None, waitMinutes=0, waitHours=0, queue=None, memory=90, memoryUnits=PERCENTAGE, getMemoryFromAnalysis=True, explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE, echoPrint=OFF, modelPrint=OFF, contactPrint=OFF, historyPrint=OFF, userSubroutine='', scratch='', resultsFormat=ODB, multiprocessingMode=DEFAULT, numCpus=1, numGPUs=0) #创建Job-1
mdb.jobs['Job-1'].writeInput(consistencyChecking=OFF) #输出inp文件

以上代码在abaqus中运行后即可得到一个完整的边坡模型，并输出inp文件。运行操作如下：

注意：需要将.rpy文件修改为.py才可以运行！

以上为脚本建立的边坡模型。
以下为利用

接下来利用循环语句循环建立模型并输出inp文件

实现思路：

inpNumbers=10   #需要生成的inp文件数量
for i in range(inpNumbers):  #利用循环语句循环输出inp文件参数生成     #利用随机数生成c、phi参数，参数需大于零模型建立（参数替换）#将生成的随机参数输入模型INP输出  #利用INP输出语句输出inp文件

接下来对参数生成详细说明：
①随机数的生成：

python中包含很多随机数生成函数（random、uniform、randint、choice、shuffle、randrange、sample等等，函数功能自行搜索）。本文用random函数进行简单的随机参数生成。

random()：该函数会随机生成一个0.0~1.0之间的浮点数

对于c、phi值，首先确定一个均值
c=15 phi=20
那么将随机数与均值结合后便可生成随机的c、phi值，代码如下：

inpNumbers=10
c=15
phi=20
for i in range(inpNumbers):  #每次循环后都将生成一个新的参数rand_c=c+random.random()rand_phi=phi+random.random()

当然，现实的c、phi值一般会围绕均值上下浮动，且服从对数正态分布（对数正态分布的参数均大于零）。本文仅实现简单的随机参数。

②将参数输入到模型中

mdb.models['Model-1'].materials['soil'].MohrCoulombPlasticity(dependencies=1, table=((36.0, 0.0, 0.5), (25.887, 0.0, 0.75), (20.0, 0.0, 1.0), (16.2343, 0.0, 1.25), (13.639, 0.0, 1.5), (11.749, 0.0, 1.75), (10.3141, 0.0, 2.0)))
mdb.models['Model-1'].materials['soil'].mohrCoulombPlasticity.MohrCoulombHardening(dependencies=1, table=((24.76, 0.0, 0.5), (16.5067, 0.0, 0.75), (12.38, 0.0, 1.0), (9.904, 0.0, 1.25), (8.25333, 0.0, 1.5), (7.07429, 0.0, 1.75), (6.19, 0.0, 2.0)))

上述代码为原模型中参数设定的代码，需要进行以下修改

#对phi进行赋值mdb.models['Model-1'].materials['soil'].MohrCoulombPlasticity(dependencies=1,table=((float(math.atan((math.tan(rand_phi*(math.pi/180))/0.5))*(180/math.pi)), 0.0, 0.5),(float(math.atan((math.tan(rand_phi*(math.pi/180))/0.75))*(180/math.pi)), 0.0, 0.75), (float(math.atan((math.tan(rand_phi*(math.pi/180))/1.0))*(180/math.pi)), 0.0, 1.0), (float(math.atan((math.tan(rand_phi*(math.pi/180))/1.25))*(180/math.pi)), 0.0, 1.25),(float(math.atan((math.tan(rand_phi*(math.pi/180))/1.25))*(180/math.pi)), 0.0, 1.5),(float(math.atan((math.tan(rand_phi*(math.pi/180))/1.75))*(180/math.pi)), 0.0, 1.75),(float(math.atan((math.tan(rand_phi*(math.pi/180))/2.0))*(180/math.pi)), 0.0, 2.0)))#对c进行赋值mdb.models['Model-1'].materials['soil'].mohrCoulombPlasticity.MohrCoulombHardening(dependencies=1, table=((float(rand_c/0.5), 0.0, 0.5), (float(rand_c/0.75), 0.0, 0.75), (float(rand_c/1.0), 0.0, 1.0), (float(rand_c/1.25), 0.0, 1.25), (float(rand_c/1.5), 0.0, 1.5), (float(rand_c/1.75), 0.0, 1.75), (float(rand_c/2.0), 0.0, 2.0)))

经过以上操作即可完成将任意输入的c、phi值赋值在模型上

③INP文件的输出

mdb.jobs['Job-1'].writeInput(consistencyChecking=OFF) #输出inp文件

以上为单一INP文件输出语句，若想批量生成INP，需要进行如下修改

for i in rang(inpNumbers):#输出inp文件inp文件名递增mdb.Job(name='Job-'+str(i+1), model='Model-1', description='', type=ANALYSIS, atTime=None, waitMinutes=0, waitHours=0, queue=None, memory=90, memoryUnits=PERCENTAGE, getMemoryFromAnalysis=True, explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE, echoPrint=OFF, modelPrint=OFF, contactPrint=OFF, historyPrint=OFF, userSubroutine='', scratch='', resultsFormat=ODB, multiprocessingMode=DEFAULT, numCpus=1, numGPUs=0)mdb.jobs['Job-'+str(i+1)].writeInput(consistencyChecking=OFF) 输出

通过以上介绍的①②③三个步骤，便可得到批量生成模型的.py文件，代码如下：

# -*- coding: mbcs -*-
#
# Abaqus/CAE Release 2020 replay file
# Internal Version: 2019_09_14-01.49.31 163176
# Run by pc on Sun Jan 14 14:58:55 2024
## from driverUtils import executeOnCaeGraphicsStartup
# executeOnCaeGraphicsStartup()
#: Executing "onCaeGraphicsStartup()" in the site directory ...
from abaqus import *
from abaqusConstants import *
import random
import numpy
import math
mdb.Model(name='Model-2', modelType=STANDARD_EXPLICIT)
#: The model "Model-2" has been created.
session.viewports['Viewport: 1'].setValues(displayedObject=None)
session.viewports['Viewport: 1'].setValues(displayedObject=None)
del mdb.models['Model-1']inpNumbers=10
c_input=float(15)
phi_input=float(20)for i in range(inpNumbers):mdb.Model(name='Model-1', modelType=STANDARD_EXPLICIT)session.Viewport(name='Viewport: 1', origin=(0.0, 0.0), width=192.67707824707, height=155.555557250977)session.viewports['Viewport: 1'].makeCurrent()session.viewports['Viewport: 1'].maximize()from caeModules import *from driverUtils import executeOnCaeStartupexecuteOnCaeStartup()session.viewports['Viewport: 1'].partDisplay.geometryOptions.setValues(referenceRepresentation=ON)s = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=200.0)g, v, d, c = s.geometry, s.vertices, s.dimensions, s.constraintss.setPrimaryObject(option=STANDALONE)s.Line(point1=(0.0, 0.0), point2=(20.0, 0.0))s.HorizontalConstraint(entity=g[2], addUndoState=False)s.Line(point1=(20.0, 0.0), point2=(20.0, 13.0))s.VerticalConstraint(entity=g[3], addUndoState=False)s.PerpendicularConstraint(entity1=g[2], entity2=g[3], addUndoState=False)s.Line(point1=(20.0, 13.0), point2=(12.0, 13.0))s.HorizontalConstraint(entity=g[4], addUndoState=False)s.PerpendicularConstraint(entity1=g[3], entity2=g[4], addUndoState=False)s.Line(point1=(12.0, 13.0), point2=(2.0, 3.0))s.Line(point1=(2.0, 3.0), point2=(0.0, 3.0))s.HorizontalConstraint(entity=g[6], addUndoState=False)s.Line(point1=(0.0, 3.0), point2=(0.0, 0.0))s.VerticalConstraint(entity=g[7], addUndoState=False)s.PerpendicularConstraint(entity1=g[6], entity2=g[7], addUndoState=False)p = mdb.models['Model-1'].Part(name='slope', dimensionality=TWO_D_PLANAR, type=DEFORMABLE_BODY)p = mdb.models['Model-1'].parts['slope']p.BaseShell(sketch=s)s.unsetPrimaryObject()p = mdb.models['Model-1'].parts['slope']session.viewports['Viewport: 1'].setValues(displayedObject=p)del mdb.models['Model-1'].sketches['__profile__']session.viewports['Viewport: 1'].view.setValues(nearPlane=43.5104, farPlane=51.9045, width=20.4584, height=14.4109, viewOffsetX=-0.00816162, viewOffsetY=0.149886)p = mdb.models['Model-1'].parts['slope']f, e, d1 = p.faces, p.edges, p.datumst = p.MakeSketchTransform(sketchPlane=f[0], sketchPlaneSide=SIDE1, origin=(12.175439, 5.508772, 0.0))s1 = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=47.7, gridSpacing=1.19, transform=t)g, v, d, c = s1.geometry, s1.vertices, s1.dimensions, s1.constraintss1.setPrimaryObject(option=SUPERIMPOSE)p = mdb.models['Model-1'].parts['slope']p.projectReferencesOntoSketch(sketch=s1, filter=COPLANAR_EDGES)s1.Line(point1=(-10.175439, -2.508772), point2=(-10.175439, -5.50877199998597))s1.VerticalConstraint(entity=g[8], addUndoState=False)s1.CoincidentConstraint(entity1=v[6], entity2=g[3], addUndoState=False)session.viewports['Viewport: 1'].view.setValues(width=22.3404, height=16.7202, cameraPosition=(10.3414, 6.87094, 47.7074), cameraTarget=(10.3414, 6.87094, 0))s1.Line(point1=(-10.175439, -2.508772), point2=(7.82456100005157, -2.508772))s1.HorizontalConstraint(entity=g[9], addUndoState=False)s1.CoincidentConstraint(entity1=v[7], entity2=g[4], addUndoState=False)p = mdb.models['Model-1'].parts['slope']f = p.facespickedFaces = f.getSequenceFromMask(mask=('[#1 ]', ), )e1, d2 = p.edges, p.datumsp.PartitionFaceBySketch(faces=pickedFaces, sketch=s1)s1.unsetPrimaryObject()del mdb.models['Model-1'].sketches['__profile__']session.viewports['Viewport: 1'].view.setValues(nearPlane=43.0377, farPlane=52.3772, width=22.9019, height=16.1322, viewOffsetX=0.035477, viewOffsetY=0.709127)session.viewports['Viewport: 1'].partDisplay.setValues(sectionAssignments=ON, engineeringFeatures=ON)session.viewports['Viewport: 1'].partDisplay.geometryOptions.setValues(referenceRepresentation=OFF)rand_c=float(c_input+random.random())rand_phi=float(phi_input+random.random())mdb.models['Model-1'].Material(name='soil')mdb.models['Model-1'].materials['soil'].Elastic(table=((100000.0, 0.33), ))mdb.models['Model-1'].materials['soil'].MohrCoulombPlasticity(dependencies=1,table=((float(math.atan((math.tan(rand_phi*(math.pi/180))/0.5))*(180/math.pi)), 0.0, 0.5),(float(math.atan((math.tan(rand_phi*(math.pi/180))/0.75))*(180/math.pi)), 0.0, 0.75), (float(math.atan((math.tan(rand_phi*(math.pi/180))/1.0))*(180/math.pi)), 0.0, 1.0), (float(math.atan((math.tan(rand_phi*(math.pi/180))/1.25))*(180/math.pi)), 0.0, 1.25),(float(math.atan((math.tan(rand_phi*(math.pi/180))/1.25))*(180/math.pi)), 0.0, 1.5),(float(math.atan((math.tan(rand_phi*(math.pi/180))/1.75))*(180/math.pi)), 0.0, 1.75),(float(math.atan((math.tan(rand_phi*(math.pi/180))/2.0))*(180/math.pi)), 0.0, 2.0)))mdb.models['Model-1'].materials['soil'].mohrCoulombPlasticity.MohrCoulombHardening(dependencies=1, table=((float(rand_c/0.5), 0.0, 0.5), (float(rand_c/0.75), 0.0, 0.75), (float(rand_c/1.0), 0.0, 1.0), (float(rand_c/1.25), 0.0, 1.25), (float(rand_c/1.5), 0.0, 1.5), (float(rand_c/1.75), 0.0, 1.75), (float(rand_c/2.0), 0.0, 2.0)))mdb.models['Model-1'].materials['soil'].mohrCoulombPlasticity.TensionCutOff(temperatureDependency=OFF, dependencies=0, table=((0.0, 0.0), ))mdb.models['Model-1'].HomogeneousSolidSection(name='soil', material='soil', thickness=None)p = mdb.models['Model-1'].parts['slope']f = p.facesfaces = f.getSequenceFromMask(mask=('[#7 ]', ), )region = regionToolset.Region(faces=faces)p = mdb.models['Model-1'].parts['slope']p.SectionAssignment(region=region, sectionName='soil', offset=0.0, offsetType=MIDDLE_SURFACE, offsetField='', thicknessAssignment=FROM_SECTION)a = mdb.models['Model-1'].rootAssemblysession.viewports['Viewport: 1'].setValues(displayedObject=a)session.viewports['Viewport: 1'].assemblyDisplay.setValues(optimizationTasks=OFF, geometricRestrictions=OFF, stopConditions=OFF)a = mdb.models['Model-1'].rootAssemblya.DatumCsysByDefault(CARTESIAN)p = mdb.models['Model-1'].parts['slope']a.Instance(name='slope-1', part=p, dependent=ON)session.viewports['Viewport: 1'].assemblyDisplay.setValues(adaptiveMeshConstraints=ON)mdb.models['Model-1'].StaticStep(name='load', previous='Initial', initialInc=0.1, matrixSolver=DIRECT, matrixStorage=UNSYMMETRIC)session.viewports['Viewport: 1'].assemblyDisplay.setValues(step='load')mdb.models['Model-1'].StaticStep(name='reduce', previous='load', initialInc=0.1, matrixSolver=DIRECT, matrixStorage=UNSYMMETRIC)session.viewports['Viewport: 1'].assemblyDisplay.setValues(step='reduce')mdb.models['Model-1'].fieldOutputRequests['F-Output-1'].setValues(variables=('S', 'PE', 'PEEQ', 'PEMAG', 'LE', 'U', 'RF', 'CF', 'CSTRESS', 'CDISP', 'FV'))session.viewports['Viewport: 1'].assemblyDisplay.setValues(loads=ON, bcs=ON, predefinedFields=ON, connectors=ON, adaptiveMeshConstraints=OFF)session.viewports['Viewport: 1'].view.setValues(width=18.803, height=13.2449, viewOffsetX=-0.354432, viewOffsetY=0.115113)session.viewports['Viewport: 1'].assemblyDisplay.setValues(step='load')session.viewports['Viewport: 1'].view.setValues(nearPlane=43.0545, farPlane=52.3604, width=22.8175, height=16.0727, viewOffsetX=1.14506, viewOffsetY=-0.247424)a = mdb.models['Model-1'].rootAssemblyf1 = a.instances['slope-1'].facesfaces1 = f1.getSequenceFromMask(mask=('[#7 ]', ), )region = regionToolset.Region(faces=faces1)mdb.models['Model-1'].BodyForce(name='Load-1', createStepName='load', region=region, comp2=-20.0)session.viewports['Viewport: 1'].assemblyDisplay.setValues(step='Initial')session.viewports['Viewport: 1'].view.setValues(width=24.1229, height=16.9922, viewOffsetX=0.561247, viewOffsetY=-0.0205853)a = mdb.models['Model-1'].rootAssemblye1 = a.instances['slope-1'].edgesedges1 = e1.getSequenceFromMask(mask=('[#242 ]', ), )region = regionToolset.Region(edges=edges1)mdb.models['Model-1'].DisplacementBC(name='BC-1', createStepName='Initial', region=region, u1=SET, u2=UNSET, ur3=UNSET, amplitude=UNSET, distributionType=UNIFORM, fieldName='', localCsys=None)a = mdb.models['Model-1'].rootAssemblye1 = a.instances['slope-1'].edgesedges1 = e1.getSequenceFromMask(mask=('[#180 ]', ), )region = regionToolset.Region(edges=edges1)mdb.models['Model-1'].DisplacementBC(name='BC-2', createStepName='Initial', region=region, u1=SET, u2=SET, ur3=UNSET, amplitude=UNSET, distributionType=UNIFORM, fieldName='', localCsys=None)session.viewports['Viewport: 1'].assemblyDisplay.setValues(mesh=ON, loads=OFF, bcs=OFF, predefinedFields=OFF, connectors=OFF)session.viewports['Viewport: 1'].assemblyDisplay.meshOptions.setValues(meshTechnique=ON)p = mdb.models['Model-1'].parts['slope']session.viewports['Viewport: 1'].setValues(displayedObject=p)session.viewports['Viewport: 1'].partDisplay.setValues(sectionAssignments=OFF, engineeringFeatures=OFF, mesh=ON)session.viewports['Viewport: 1'].partDisplay.meshOptions.setValues(meshTechnique=ON)p = mdb.models['Model-1'].parts['slope']p.seedPart(size=0.5, deviationFactor=0.1, minSizeFactor=0.1)p = mdb.models['Model-1'].parts['slope']f = p.facespickedRegions = f.getSequenceFromMask(mask=('[#7 ]', ), )p.setMeshControls(regions=pickedRegions, elemShape=QUAD, technique=SWEEP)elemType1 = mesh.ElemType(elemCode=CPE4, elemLibrary=STANDARD)elemType2 = mesh.ElemType(elemCode=CPE3, elemLibrary=STANDARD)p = mdb.models['Model-1'].parts['slope']f = p.facesfaces = f.getSequenceFromMask(mask=('[#7 ]', ), )pickedRegions =(faces, )p.setElementType(regions=pickedRegions, elemTypes=(elemType1, elemType2))p = mdb.models['Model-1'].parts['slope']p.generateMesh()session.viewports['Viewport: 1'].partDisplay.setValues(mesh=OFF)session.viewports['Viewport: 1'].partDisplay.meshOptions.setValues(meshTechnique=OFF)session.viewports['Viewport: 1'].partDisplay.geometryOptions.setValues(referenceRepresentation=ON)p = mdb.models['Model-1'].parts['slope']f = p.facesfaces = f.getSequenceFromMask(mask=('[#7 ]', ), )p.Set(faces=faces, name='slope')#: The set 'slope' has been created (3 faces).import jobmdb.models['Model-1'].keywordBlock.synchVersions(storeNodesAndElements=False)mdb.models['Model-1'].keywordBlock.insert(33, """*initial conditions,type=field,variable=1slope-1.slope,0.5""")mdb.models['Model-1'].keywordBlock.insert(51, """*field,variable=1slope-1.slope,2""")a = mdb.models['Model-1'].rootAssemblysession.viewports['Viewport: 1'].setValues(displayedObject=a)session.viewports['Viewport: 1'].assemblyDisplay.setValues(mesh=OFF)session.viewports['Viewport: 1'].assemblyDisplay.meshOptions.setValues(meshTechnique=OFF)mdb.Job(name='Job-'+str(i+1), model='Model-1', description='', type=ANALYSIS, atTime=None, waitMinutes=0, waitHours=0, queue=None, memory=90, memoryUnits=PERCENTAGE, getMemoryFromAnalysis=True, explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE, echoPrint=OFF, modelPrint=OFF, contactPrint=OFF, historyPrint=OFF, userSubroutine='', scratch='', resultsFormat=ODB, multiprocessingMode=DEFAULT, numCpus=1, numGPUs=0)mdb.jobs['Job-'+str(i+1)].writeInput(consistencyChecking=OFF)del mdb.models['Model-1']

4.注意事项

在编写脚本的过程中遇到了以下问题：
1、abaqus中执行python代码时，若使用了第三方库，需要令第三方库与abaqus中的python版本保持一致，否则不可使用，本文中原本使用numpy中的numpy.pi、numpy.arctan、numpy.tan进行折减参数的计算，但由于以上原因，不得不将numpy库变为math库（内置库），相应的函数变为math.pi、math.atan、math.tan。若后续必须使用第三方库才可完成脚本，切记需要保持第三方库和abaqus中的python版本一致。
解决方案：

https://zhuanlan.zhihu.com/p/59440515

以上方案未经验证，如不能实现，自行百度。

2、在循环生成INP文件时，如果不删除上一个Model，可能会出现报错情况，但Abaqus必须存在至少一个Model，所以本文先建立了一个Model-2

mdb.Model(name='Model-2', modelType=STANDARD_EXPLICIT)

此模型并未使用，仅用来占位。
然后在每次生成INP后删除循环中的Model

del mdb.models['Model-1']

这种方法较为粗糙，时间有限，后续会优化该部分代码。

3、脚本中存在许多冗余代码，若有兴趣，可以逐行核对，删除冗余代码后，程序效率会相应提高。

4、对应py文件已放在资源中，同本文中代码相同。