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

利用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 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__']
    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.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=(-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.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'].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.faces
    faces = 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'].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')
    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'].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')
    session.viewports['Viewport: 1'].view.setValues(width=24.1229, height=16.9922, 
        viewOffsetX=0.561247, viewOffsetY=-0.0205853)
    a = mdb.models['Model-1'].rootAssembly
    e1 = a.instances['slope-1'].edges
    edges1 = 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'].rootAssembly
    e1 = a.instances['slope-1'].edges
    edges1 = 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.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.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(51, """
    *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-'+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文件已放在资源中，同本文中代码相同。