当前位置：服务支持 > 软件文章 > [Abaqus揭秘]Xfem功能大探索

[Abaqus揭秘]Xfem功能大探索

阅读数 7

啥叫XFEM

你可能会想到

还有这个

其实，还有这个

那么，既然世界上有商业软件这个东西，作为纵横工业界与科研界的我们，还是比较喜欢不用编代码的东西。

这是个复杂的问题，官方是这样说的：是一种解决断裂力学问题的新的有限元方法，其理论最早亍 1999 年，由美国西北大学的教授 Belyschko 和 Black 首次提出，主要是采用独立亍网格剖分的思想解决有限元中的裂纹扩展问题，在保留传统有限元所有优点的同时，并不需要对结构内部存在的裂纹等缺陷进行网格划分。

等等这个重要吗？

这个当然不重要！！作为工科生我们只关心两件事

第一，这个玩意怎么搞，有个啥用

第二，这个玩意怎么用

但是，作为一名搞科研的人，还是要把原理大致阐述一下，就是ABAQUS 在6.7版本引入XFEM这种算法可以很好地描述，断裂的特性，应力、裂纹扩展的情况，在机械、岩石有着广阔的应用，那么作为一个工科生，我可以负责人的讲，前方高能预警

要是用还是要干这两档子事，

第一，选择模型中可能出现裂纹的区域，将其单元设置为具有扩展有限元性质的富集

单元。

第二，选择合适的破坏准则，使得单元在达到条件时发生破坏，裂纹得以扩展。

首先，在材料中添加Maxps Damage ，这个是传说中的损伤系数，用来说明破坏

准则，Damage Evolution是损伤演化参数主要用来损伤的发展情况

Damage Stabilization Cohesive 是损伤稳定性系数，改善收敛用的

第二，档子事就是告诉电脑损伤区在哪？在哪，自己定义

在相互作用（interaction）模块，做一下剧透，我一直对这个Special很好奇，所以我就把所以的选项都试了试，这个Special是个神奇的东西，第一个质量特性，可以定义什么质量、重心和转动惯量啊啥的，第二个我一会讲，第三个，Spring/Dashpots是用来定义两点之间的弹簧、阻尼啥的，这个弹簧还可以是什么六个自由度的、非线性的很好用，第四个，就是一个链接单元，机械里经常用，我们很少涉猎。

当然， ABAQUS 允许在裂缝发展区域中插入一个初始裂缝，这个初始裂缝可以在 Assembly里进行组装。

当然这个凡是都有缺陷，所以虽然Abaqus很好很强大，但是这也是万能的，为啥呢。因为

其实单元内部的位移函数（形函数）可以是任意形状的，但大多数的计算软件都采用了多项式戒者揑值多项式作为手段来描述单元内部的位移场，这是因为采用这种方法更加便亍在编程中进行处理。但是这种方法的缺点就是，由亍形函数的连续性，导致单元内部不可能存在间断。直到 Belytschko 提出采用水平集函数作为手段这种扩充形函数能够描述单元内位移场在裂缝两边的跳跃性，

但是由于裂缝存在于单元内部，其扩展独立不其他单元，使得计算变得高效。但是这种方法也存在一些问题，XFEM 采用的形函数模式会导致其求解方程很容易接近线形相关，极大的增加了收敛难度，因此导致 XFEM 方法一直没有办法得到很好的推广。其实 ABAQUS 在集成 XFEM 方法时做了大量的简化，目的都是减小求解的难度， ABAQUS 的帮助在介绍 XFEM 的时候其实说：

The following limitations exist with an enriched feature:

几个限制：

An enriched element cannot be intersected by more than one crack.
不能模拟交叉裂缝
A crack is not allowed to turn more than 90° in one increment during an analysis.
裂缝扩展方向不能偏转超过90
Only asymptotic crack-tip fields in an isotropic elastic material are considered for a stationary crack.
对于制定裂缝的方法来讲，裂尖渐近位移场计算只能使用各项异性材料
Adaptive remeshing is not supported.
自适应网格也不能用
Composite solid elements are not supported.
复合材料单元也不能用

此外，计算过程中很容易发现裂缝是丌能停留在单元内部的，这说明了ABAQUS 放弃了单元内部对裂尖的描述。同时，由亍 ABAQUS 在计算 XFEM 的损伤时采用的是基亍能量释放率的 paris 法则，虽然这是基亍弹塑性断裂力学的经典手段，但由亍承认了裂尖位置的塑性效应，使得在模拟损伤时也只能对低周疲劳能有比较好的近似。

总之，尽管有如此之多的限制，但 XFEM 方法依然为有限元领域注入了一剂新鲜的空气，它所提供的是一种广义有限元的实现过程，为有限元理论的发展提供了一种新的思维方式。

这是极官方的评价

考虑大家基础比较薄弱，下面较大家如何添加裂缝:

第一步，做一个Part

第二部建一个裂缝

第三步，加上位移荷载啥的

第四步，编辑裂缝

第五步，定义材料属性

第六步，定义求解步，一般规定个时间啥的

第七步，为了保证能够看清楚裂缝扩展的过程，要把求解增量步变小

第七步，定义损伤稳定系数

第八步，选取输出场变量

第九步，画网格

第十步，求解

如果还是不明白

咱也有法子，你把下边的代码丢到最下面的那个命令框里，就好了。

计算结果：

# === Modules ===

from part import *
from material import *
from section import *
from assembly import *
from step import *
from regionToolset import *
from interaction import *
from load import *
from mesh import *
from job import *
from sketch import *
from visualization import *
from connectorBehavior import *


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

# === Parameters ===

modelName = 'modeII_2d_lefm'

length = 3.0                              # plate half-length
height = 6.0                              # plate height
width  = 1.0                              # plate width

crack_len = 1.5                           # crack length
crack_y   = 0.05                          # crack offset in y-direction

YM    = 2.10e+11                          # Young's modulus
NU    = 0.3                               # Poisson's ratio
MAXPS = 8.44e7                            # Damage initiation
DTOL  = 0.05                              # Damage tolerance
GI    = 42200.0                           # Fracture energy
ETA   = 1.0                               # Power-law exponent


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

# === Create model ===

Mdb()

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

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


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

# === Create parts ===

plateSketch = plateModel.ConstrainedSketch(name='plateProfile',sheetSize=height)

plateSketch.setPrimaryObject(option=STANDALONE)
plateSketch.sketchOptions.setValues(decimalPlaces=3)

plateSketch.Line(point1=(0.0, -height/2.0), point2=(0.0, height/2.0))
plateSketch.Line(point1=(0.0, height/2.0), point2=(length, height/2.0))
plateSketch.Line(point1=(length, height/2.0), point2=(length, -height/2.0))
plateSketch.Line(point1=(length, -height/2.0), point2=(0.0, -height/2.0))

platePart = plateModel.Part(dimensionality=TWO_D_PLANAR, name='plate', 
                type=DEFORMABLE_BODY)
platePart.BaseShell(sketch=plateModel.sketches['plateProfile'])

plateSketch.unsetPrimaryObject()
del plateSketch


# Part for crack geometry

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'])

crackSketch.unsetPrimaryObject()
del crackSketch


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


# === Create geometry sets ===

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


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

# === Define material and section properties ===

plateMatl = plateModel.Material(name='elas')

plateMatl.Elastic(table=((YM, NU), ))
plateMatl.MaxpsDamageInitiation(table=((MAXPS, ), ), tolerance=DTOL)

plateModel.HomogeneousSolidSection(material='elas', name='solid', thickness=width)


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

# === Assign section and orientation ===

#platePart.MaterialOrientation(fieldName='', localCsys=None, orientationType=GLOBAL, region=
#    platePart.sets['All'], stackDirection=STACK_3)

platePart.SectionAssignment(region=platePart.sets['All'], sectionName='solid')


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

# === Assign mesh controls and mesh plate ===

# Element type
platePart.setElementType(elemTypes=(ElemType(elemCode=CPE4, elemLibrary=STANDARD), 
    ElemType(elemCode=CPE3, elemLibrary=STANDARD)), regions=platePart.sets['All'])

# Mesh technique
platePart.setMeshControls(elemShape=QUAD, regions=platePart.faces[:], technique=STRUCTURED)


# Seed mesh
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()

# === End part ===


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

# === Assemble ===

plateModel.rootAssembly.DatumCsysByDefault(CARTESIAN)
plateModel.rootAssembly.Instance(dependent=ON, name='plate_1', part=platePart)
plateModel.rootAssembly.Instance(dependent=ON, name='crack_1', part=crackPart)

# Reference points for displacement application
rp_db = plateModel.rootAssembly.ReferencePoint(point=(length, -height/2.0, 
    0.0))
plateModel.rootAssembly.features.changeKey(fromName='RP-1', toName='db')


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

# === Create assembly sets ===

v1 = (plateModel.rootAssembly.referencePoints[rp_db.id], )
plateModel.rootAssembly.Set(name='bdisp', referencePoints=v1)

# === End assembly ===


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

# === Create constraint equation ===

plateModel.Equation(name='ce_bot', terms=((1.0, 'plate_1.bottom', 1), (-1.0, 'bdisp', 1)))

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

# === Create step ===

plateModel.StaticStep(timePeriod=0.483,initialInc=0.01, maxInc=0.05, maxNumInc=10000,
    minInc=1e-012, name='Static', nlgeom=ON, previous='Initial')

plateModel.steps['Static'].control.setValues(
    allowPropagation=OFF,  resetDefaultValues=OFF, 
    timeIncrementation=(4.0, 8.0, 9.0, 16.0, 10.0, 4.0, 12.0, 20.0, 6.0, 3.0, 
    50.0))

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

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


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

# === Apply boundary conditions ===

plateModel.TabularAmplitude(name='Amp-1', timeSpan=STEP, 
    smooth=SOLVER_DEFAULT, data=((0.0, 0.0), (1, 1.0)))
plateModel.DisplacementBC( createStepName=
    'Static', distributionType=UNIFORM, fieldName='', fixed=OFF, localCsys=None
    , name='rp',amplitude='Amp-1',  region=plateModel.rootAssembly.sets['bdisp'], u1=0.01, 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, ur3=UNSET)

plateModel.DisplacementBC( createStepName=
    'Static', distributionType=UNIFORM, fieldName='', fixed=OFF, localCsys=None
    , name='top',amplitude='Amp-1',  region=plateModel.rootAssembly.instances['plate_1'].sets['top']
    , u1=-0.01, u2=0, ur3=UNSET)


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

# === Define enrichment and initial crack ===

edges = plateModel.rootAssembly.instances['crack_1'].edges
e1 = edges.findAt(( (crack_len/2.0, crack_y, 0.0), ))
crackLocation = Region(edges=e1)

plateModel.ContactProperty('contact')

plateModel.rootAssembly.engineeringFeatures.XFEMCrack( name='enr1',
    crackDomain=plateModel.rootAssembly.instances['plate_1'].sets['All']
    , interactionProperty='contact', crackLocation=crackLocation)

 
plateModel.interactionProperties['contact'].FractureCriterion(
    initTable=((GI, GI, GI, 1.0, 1.0, 1.0), ), 
    mixedModeBehavior=POWER, tolerance=0.2, viscosity=5e-05)

plateModel.rootAssembly.engineeringFeatures.cracks['enr1'].setValues(interactionProperty='contact')
#---------------------------------------------------------------------------------

# === Create Job === 

mdb.Job(model=modelName, name='modeII_lefm_xfem_cpe4',
    description='Crack propagation in a plate under mixed-mode loading (XFEM)')


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

搞定。。。

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删

返回上级列表

联系我们

，获取更多内容

Abaqus实战经验大放送：论坛精华整理（一）

解锁ABAQUS高效实用技巧集锦

SIMULIA Abaqus R2023x Explicit 新功能揭秘 | 达索系统百世慧分享