Python abaqus实现二维裂纹扩展(XFEM)

自己做的一个简单的案例，用python用于abaqus的模拟比较省时间，留作参考

在abaqus中，点击file——run script——运行脚本就可以了，abaqus工作目录下的.rpy文件是abaqus的日志，进行每一步操作也有相应的代码，可以直接拷贝过来用

#!/usr/bin/env/ python
# -*- coding:utf-8 -*-

# === 导入必要的模块 ===
from abaqus import *
from abaqusConstants import *
import mesh
import regionToolset
import job
#----------------------------------------------------------------

# === 设置参数 ====
modelName = "mixed_mode_2d"

length = 3.0
height = 6.0
width = 1.0

crack_len = 1.5
crack_y = 0.05                          # 裂缝在y方向平移的距离

YM = 2.10E+11                           # 杨氏模量(Pa)
MU = 0.3                                # 泊松比
MAXPS = 2.2E+08                        # 损伤开始的最大主应力
DTOL = 0.05                             # 损伤容限
GI = 42200.0                            # 断裂能
ETA = 1.0                               # Power-law的指数

#----------------------------------------------------------------

# === 创建模型和视图设置 ===

Mdb()

viewportName = session.Viewport(name=modelName)
viewportName.makeCurrent()
viewportName.maximize()

plateModel = mdb.Model(name=modelName)
del mdb.models['Model-1']

#----------------------------------------------------------------

# === 创建部件，首先是平板部件 ===
plateSketch = plateModel.ConstrainedSketch(name='plateProfile',sheetSize=height)
plateSketch.rectangle(point1=(0.0, -height/2.0),point2=(length,height/2.0))
platePart = plateModel.Part(dimensionality=TWO_D_PLANAR,name='plate',type=DEFORMABLE_BODY)
platePart.BaseShell(sketch=plateModel.sketches['plateProfile'])


# === 创建裂缝几何部件 ===
crackSketch = plateModel.ConstrainedSketch(name='crackProfile',sheetSize=height)
crackSketch.Line(point1=(0.0,crack_y),point2=(crack_len,crack_y))
crackPart=plateModel.Part(dimensionality=TWO_D_PLANAR,name='crack',type=DEFORMABLE_BODY)
crackPart.BaseWire(sketch=plateModel.sketches['crackProfile'])

#----------------------------------------------------------------

# === 创建一些集合（Set）===
platePart.Set(faces=platePart.faces[:],name='ALL')

e1 = platePart.edges.findAt(((length/2.0, -height/2.0, 0.0),))
platePart.Set(edges=e1,name='bottom')

e1 = platePart.edges.findAt(((length/2.0, +height/2.0, 0.0),))
platePart.Set(edges=e1,name='top')

platePart.Set(faces=crackPart.faces[:],name='crack')

#----------------------------------------------------------------

# === 定义材料和截面属性 ===
plateMat1 = plateModel.Material(name='elas')

plateMat1.Elastic(table=((YM,MU),))
plateMat1.MaxpsDamageInitiation(table=((MAXPS,),),tolerance=DTOL)
plateMat1.maxpsDamageInitiation.DamageEvolution(
        mixedModeBehavior=POWER_LAW, power=ETA, table=((GI,GI,GI),),type=ENERGY)
plateMat1.maxpsDamageInitiation.DamageStabilizationCohesive(
        cohesiveCoeff=0.001)
plateModel.HomogeneousSolidSection(material='elas',name='solid',thickness=width)

#----------------------------------------------------------------

# === 赋给材料截面属性和方向 ===
platePart.MaterialOrientation(fieldName='', localCsys=None, orientationType=GLOBAL,region=
        platePart.sets['ALL'],stackDirection=STACK_3)
        
platePart.SectionAssignment(region=platePart.sets['ALL'],sectionName='solid')

#----------------------------------------------------------------

# ===网格控制和划分网格===

# 单元类型
elemType1 = mesh.ElemType(elemCode=CPE4R, elemLibrary=STANDARD, 
    secondOrderAccuracy=OFF, hourglassControl=DEFAULT, 
    distortionControl=DEFAULT)
elemType2 = mesh.ElemType(elemCode=CPE3, elemLibrary=STANDARD, 
    secondOrderAccuracy=OFF, distortionControl=DEFAULT)
p1 = mdb.models['mixed_mode_2d'].parts['plate']
f = p1.faces
faces = f.getSequenceFromMask(mask=('[#1 ]', ), )
pickedRegions =(faces, )
p1.setElementType(regions=pickedRegions, elemTypes=(elemType1, elemType2))



# 单元划分技术
pickedRegions = f.getSequenceFromMask(mask=('[#1 ]', ), )
p1.setMeshControls(regions=pickedRegions, elemShape=QUAD, technique=STRUCTURED)

# 网格种子
ex = 30
ey = 60
e1 = platePart.edges.findAt(( (length/2.0, -height/2.0, 0.0), ))
platePart.seedEdgeByNumber(edges=e1, number=ex)

e1 = platePart.edges.findAt(( (length, 0.0, 0.0), ))
platePart.seedEdgeByNumber(edges=e1, number=ey)

platePart.generateMesh()
# === 部件部分结束 ===

#----------------------------------------------------------------


# === 装配 ===
plateModel.rootAssembly.DatumCsysByDefault(CARTESIAN)
plateModel.rootAssembly.Instance(dependent=ON, name='plate_ 1', part=platePart)
plateModel.rootAssembly.Instance(dependent=ON, name='crack_ 1', part=crackPart)

# 创建参考点，以用于施加位移边界条件
rp_db = plateModel.rootAssembly.ReferencePoint(point=(length, -height/2.0, 0.0))
plateModel.rootAssembly.features.changeKey(fromName='RP-1', toName='db')

#----------------------------------------------------------------

# ===创建装配体下的集台(set) ===
v1 = (plateModel. rootAssembly .referencePoints[rp_db.id], )
plateModel.rootAssembly.Set(name='bdisp', referencePoints=v1)

# ===装配结束===

#----------------------------------------------------------------

# === 创建的束方程，以用于施加边界条件 ===
mdb.models['mixed_mode_2d'].Equation(name='ce_bot', terms=((1.0, 
    'plate_ 1.bottom', 1), (-1.0, 'bdisp', 1)))
#plateModel.Equation(name='ce_bot', terms=((1.0, 'plate_1.bottom',1),(-1.0,'bdisp',1)))

#----------------------------------------------------------------

# === 创建分析步和输出设置 ===
plateModel.StaticStep(initialInc=0.005,maxInc=0.01, maxNumInc=10000,
        minInc=1e-09,name='Static', nlgeom=ON, previous='Initial')
        
plateModel.steps['Static']. control. setValues(
        allowPropagation=OFF,discontinuous=ON,resetDefaultValues=OFF,
        timeIncrementation=(8.0, 10.0, 9.0,16.0, 10.0, 4.0, 12.0, 20.0, 6.0, 3.0,
        50.0))

mdb.models['mixed_mode_2d'].fieldOutputRequests['F-Output-1'].setValues(
    variables=('S', 'LE', 'U', 'RF', 'PHILSM', 'STATUSXFEM'))           

mdb.models['mixed_mode_2d'].HistoryOutputRequest(createStepName='Static',
        name='H-0utput-2', rebar=EXCLUDE, region=
        mdb.models[modelName].rootAssembly.sets['bdisp'],sectionPoints
        =DEFAULT, variables=('U1', 'RF1'))
       
#----------------------------------------------------------------

# === 施加边界条件 ===
plateModel. DisplacementBC( amplitude=UNSET, createStepName=
        'Static', distributionType=UNIFORM, fieldName='', fixed=OFF, localCsys=None
        ,name='rp', region=plateModel.rootAssembly.sets['bdisp'], u1=0.0027,u2=
        UNSET, ur3=UNSET)
        
plateModel. DisplacementBC ( amplitude=UNSET, createStepName=
        'Static', distributionType=UNIFORM, fieldName='', fixed=OFF, localCsys=None
        ,name='bot', region=plateModel.rootAssembly.instances['plate_ 1'].sets['bottom']
        ,u1=UNSET,u2=-0.00162,ur3=UNSET)
        
plateModel. DisplacementBC( amplitude=UNSET, createStepName=
        'Static', distributionType=UNIFORM, fieldName='', fixed=OFF, localCsys=None
        ,name='top',region=plateModel.rootAssembly.instances['plate_ 1'].sets['top']
        ,u1=-0.0027,u2=0.00162, ur3=UNSET)

#----------------------------------------------------------------

# === 定义扩充和初始裂纹 ===
plateModel.ContactProperty('contact')
a = mdb.models['mixed_mode_2d'].rootAssembly
crackDomain = a.instances['plate_ 1'].sets['ALL']
a = mdb.models['mixed_mode_2d'].rootAssembly
e1 = a.instances['crack_ 1'].edges
edges1 = e1.getSequenceFromMask(mask=('[#1 ]', ), )
crackLocation = regionToolset.Region(edges=edges1)
a = mdb.models['mixed_mode_2d'].rootAssembly
a.engineeringFeatures.XFEMCrack(name='enrl', crackDomain=crackDomain, 
    interactionProperty='contact', crackLocation=crackLocation)
#----------------------------------------------------------------

# ===创建Job并提交===
mdb.Job(model=modelName, name='mixed_mode_xfem_cpe4',
        description='Crack propagation in a plate under mixed-mode loading (XFEM)')
mdb.jobs['mixed_mode_xfem_cpe4'].submit()
#----------------------------------------------------------------