星辰技文|一步步教你如何使用35行代码生成ABAQUS二维随机颗粒模型

前面推荐了一些ABAQUS二次开发小工具，不知道大家是否已经安装使用。

后面以一些小案例带大家熟悉ABAQUS前后处理相关的Python库，以及使用技巧。

星哥开发的插件大多集中在非均质相关断裂问题，相信关注公众号的很多朋友也都是做这方面，那么我们从最初始的非均质几何模型的案例出发，来演示一个随机颗粒模型的代码编写的全过程，效果如下所示：

在这个案例中，最大的帮手是PythonReader,它能让初学者能轻松了解GUI界面中的每个操作对应的代码段是什么，比如，点击新建文件按钮，会弹出以下代码：

Mdb()
#: A new model database has been created.
#: The model "Model-1" has been created.
session.viewports['Viewport: 1'].setValues(displayedObject=None)

其中第一行就是创建新数据模型的命令，第二行和第三行均为注释，描述模型的信息，第4行则是视图设置当前显示对象为None，即空。大多数的前处理操作均可以用这种方式进行录制，只需要了解一些Python基础知识和接口代码的风格和结构，小白也能轻松上手。

为实现随机颗粒模型的代码编写，我们分3步进行：

• • 首先通过录制获得ABAQUS中建模的相关代码。

• • 然后修改相应代码，删除无用代码，实现代码的参数化。

• • 最后在代码中添加循环和干涉判断，实现多个颗粒的随机投递。

第一步：录制

界面中操作 ：

1）新建部件，Modeling Space定义为2D Planar，其它保持默认，进入草图

s = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=200.0)

2）草图内绘制一个矩形，输入坐标[0,0],[50,50]

g, v, d, c = s.geometry, s.vertices, s.dimensions, s.constraints
s.setPrimaryObject(option=STANDALONE)
s.rectangle(point1=(0.0, 0.0), point2=(50.0, 50.0))

3）完成草图获得部件

p = mdb.models['Model-1'].Part(name='Part-1', dimensionality=TWO_D_PLANAR, 
    type=DEFORMABLE_BODY)
p = mdb.models['Model-1'].parts['Part-1']
p.BaseShell(sketch=s)
s.unsetPrimaryObject()
p = mdb.models['Model-1'].parts['Part-1']
session.viewports['Viewport: 1'].setValues(displayedObject=p)
del mdb.models['Model-1'].sketches['__profile__']

4）使用面刨分工具

p = mdb.models['Model-1'].parts['Part-1']
f, e, d1 = p.faces, p.edges, p.datums
t = p.MakeSketchTransform(sketchPlane=f[0], sketchPlaneSide=SIDE1, origin=(
    25.0, 25.0, 0.0))
s1 = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', 
    sheetSize=141.42, gridSpacing=3.53, transform=t)
g, v, d, c = s1.geometry, s1.vertices, s1.dimensions, s1.constraints
s1.setPrimaryObject(option=SUPERIMPOSE)
p = mdb.models['Model-1'].parts['Part-1']
p.projectReferencesOntoSketch(sketch=s1, filter=COPLANAR_EDGES)

5）任意位置绘制一个圆

s1.CircleByCenterPerimeter(center=(-7.06, 10.59), point1=(-4.4125, 10.59))

6）完成草图，实现颗粒的绘制

p = mdb.models['Model-1'].parts['Part-1']
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__']

第二步：参数化

有一定英语基础，并了解基本的Python类语法的概念，我们就能轻松读懂每一行代码，从而了解哪些代码可以被简化或删除：

from abaqus import *
from abaqusConstants import *
# 创建草图，保留
s = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=200.0)
# 定义的变量g,v,d,c均未调用，该行可删除
g, v, d, c = s.geometry, s.vertices, s.dimensions, s.constraints
# 窗口显示草图，删除不影响
s.setPrimaryObject(option=STANDALONE)
# 创建矩形，保留
s.rectangle(point1=(0.0, 0.0), point2=(50.0, 50.0))
# 创建部件，保留
p = mdb.models['Model-1'].Part(name='Part-1', dimensionality=TWO_D_PLANAR, type=DEFORMABLE_BODY)
# p变量已经定义，并指向同一个对象，该行可删除
p = mdb.models['Model-1'].parts['Part-1']
# 部件添加基础面，保留
p.BaseShell(sketch=s)
# 窗口隐藏草图，与显示行对应，删除不影响
s.unsetPrimaryObject()
# p变量已经定义，并指向同一个对象，该行可删除
p = mdb.models['Model-1'].parts['Part-1']
# 视图显示部件，可移到最后
session.viewports['Viewport: 1'].setValues(displayedObject=p)
# 删除临时草图，保留
del mdb.models['Model-1'].sketches['__profile__']
# p变量已经定义，并指向同一个对象，该行可删除
p = mdb.models['Model-1'].parts['Part-1']
# 定义的变量e,d1均未调用，该行可简化
f, e, d1 = p.faces, p.edges, p.datums
# 创建草图投影面，默认草图中心会在f[0]的几何中心，需与全局坐标一致，origin可修改为0，0，保留
t = p.MakeSketchTransform(sketchPlane=f[0], sketchPlaneSide=SIDE1, origin=(
    25.0, 25.0, 0.0))
# 创建草图
s1 = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', 
    sheetSize=141.42, gridSpacing=3.53, transform=t)
# 定义的变量g,v,d,c均未调用，该行可删除
g, v, d, c = s1.geometry, s1.vertices, s1.dimensions, s1.constraints
# 窗口显示草图，删除不影响
s1.setPrimaryObject(option=SUPERIMPOSE)
# p变量已经定义，并指向同一个对象，该行可删除
p = mdb.models['Model-1'].parts['Part-1']
# 原有几何投影到草图s1上，保留
p.projectReferencesOntoSketch(sketch=s1, filter=COPLANAR_EDGES)
# 绘制圆，保留
s1.CircleByCenterPerimeter(center=(-7.06, 10.59), point1=(-4.4125, 10.59))
# p变量已经定义，并指向同一个对象，该行可删除
p = mdb.models['Model-1'].parts['Part-1']
# 创建面对象，该变量已经定义，且指向的对象不变，可删除
f = p.faces
# 选择刨切的面对象pickedFaces，保留
pickedFaces = f.getSequenceFromMask(mask=('[#1 ]', ), )
# 定义的变量e1, d2均未调用，该行可删除
e1, d2 = p.edges, p.datums
# 按草图s1刨切选择的面对象pickedFaces，保留
p.PartitionFaceBySketch(faces=pickedFaces, sketch=s1)
# 窗口隐藏草图，与显示行对应，删除不影响
s1.unsetPrimaryObject()
# 删除临时草图，保留
del mdb.models['Model-1'].sketches['__profile__']

同时代码中大量出现“mdb.models['Model-1']”，可创建变量model进行代替，简化代码，清理后的代码：

from abaqus import *
from abaqusConstants import *
# 创建model变量
model = mdb.models['Model-1']
# 创建草图，保留
s = model.ConstrainedSketch(name='__profile__', sheetSize=200.0)
# 创建矩形，保留
s.rectangle(point1=(0.0, 0.0), point2=(50.0, 50.0))
# 创建部件，保留
p = model.Part(name='Part-1', dimensionality=TWO_D_PLANAR, type=DEFORMABLE_BODY)
# 部件添加基础面，保留
p.BaseShell(sketch=s)
# 删除临时草图，保留
del model.sketches['__profile__']
# 定义的变量e,d1均未调用，该行可简化
f = p.faces
# 创建草图投影面，默认草图中心会在f[0]的几何中心，需与全局坐标一致，origin可修改为0，0，保留
t = p.MakeSketchTransform(sketchPlane=f[0], sketchPlaneSide=SIDE1, origin=(0, 0, 0))
# 创建草图
s1 = model.ConstrainedSketch(name='__profile__', 
    sheetSize=141.42, gridSpacing=3.53, transform=t)
# 原有几何投影到草图s1上，保留
p.projectReferencesOntoSketch(sketch=s1, filter=COPLANAR_EDGES)
# 绘制圆，保留
s1.CircleByCenterPerimeter(center=(-7.06, 10.59), point1=(-4.4125, 10.59))
# 选择刨切的面对象pickedFaces，保留
pickedFaces = f.getSequenceFromMask(mask=('[#1 ]', ), )
# 按草图s1刨切选择的面对象pickedFaces，保留
p.PartitionFaceBySketch(faces=pickedFaces, sketch=s1)
# 删除临时草图，保留
del model.sketches['__profile__']
# 视图显示部件，可移到最后
session.viewports['Viewport: 1'].setValues(displayedObject=p)

模型中可参数化的变量包含:

• • 部件名称 partName

• • 部件宽度 width

• • 部件高度 height

• • 颗粒中心 center_x,center_y

• • 颗粒半径 radius

因此单颗粒的参数化代码如下：

from abaqus import *
from abaqusConstants import *
partName = "Part-1"     #部件名称
width = 50              #部件宽度
height=50               #部件高度
radius = 5              #颗粒半径
center_x = 10           #颗粒中心坐标x
center_y = 20           #颗粒中心坐标y
model = mdb.models['Model-1']
s = model.ConstrainedSketch(name='__profile__', sheetSize=200.0)
s.rectangle(point1=(0.0, 0.0), point2=(width, height))
p = model.Part(name=partName, dimensionality=TWO_D_PLANAR, type=DEFORMABLE_BODY)
p.BaseShell(sketch=s)
del model.sketches['__profile__']
f = p.faces
t = p.MakeSketchTransform(sketchPlane=f[0], sketchPlaneSide=SIDE1, origin=(0, 0, 0))
s1 = model.ConstrainedSketch(name='__profile__', 
    sheetSize=141.42, gridSpacing=3.53, transform=t)
p.projectReferencesOntoSketch(sketch=s1, filter=COPLANAR_EDGES)
s1.CircleByCenterPerimeter(center=(center_x, center_y), 
                   point1=(center_x, center_y+radius))
pickedFaces = f.getSequenceFromMask(mask=('[#1 ]', ), )
p.PartitionFaceBySketch(faces=pickedFaces, sketch=s1)
del model.sketches['__profile__']
session.viewports['Viewport: 1'].setValues(displayedObject=p)

第三步：循环随机

在上面单颗粒参数化模型基础上，添加循环并随机生成颗粒中心坐标，即可实现随机多颗粒的生成。为了避免颗粒之间相互重叠，需要将新生成的坐标与原有坐标进行距离判断，二者距离需要大于2倍圆半径，具体代码如下：

import random
from abaqus import *
from abaqusConstants import *
partName = "Part-1"     #部件名称
width = 100             #部件宽度
height=50               #部件高度
radius = 5              #颗粒半径
rockNum = 20            #颗粒数量
model = mdb.models['Model-1']
s = model.ConstrainedSketch(name='__profile__', sheetSize=200.0)
s.rectangle(point1=(0.0, 0.0), point2=(width, height))
p = model.Part(name=partName, dimensionality=TWO_D_PLANAR, type=DEFORMABLE_BODY)
p.BaseShell(sketch=s)
del model.sketches['__profile__']
f = p.faces
t = p.MakeSketchTransform(sketchPlane=f[0], sketchPlaneSide=SIDE1, origin=(0, 0, 0))
s1 = model.ConstrainedSketch(name='__profile__', 
    sheetSize=141.42, gridSpacing=3.53, transform=t)
p.projectReferencesOntoSketch(sketch=s1, filter=COPLANAR_EDGES)
tryTimes = 0            #引入尝试次数，避免进入死循环
rockCenterList = []     #用于存储已经投递成功的颗粒中心
while tryTimes<1000 and len(rockCenterList)<rockNum:
    #在box范围内随机生成颗粒中心坐标
    center_x = random.uniform(radius,width-radius)
    center_y = random.uniform(radius,height-radius)
    #循环判断与其它颗粒距离是否小于目标值
    mark = 1    #用于记录是否符合要求
    for x,y in rockCenterList:
        dist = sqrt((x-center_x)**2+(y-center_y)**2)
        # 距离判断，一个距离小于2倍半径，则直接退出该层循环
        if dist<2*radius:
            mark = 0
            break
    if mark == 1:
        tryTimes = 0                                            #投递成功，重新计数
        rockCenterList.append([center_x,center_y])              #更新rockCenterList
        s1.CircleByCenterPerimeter(center=(center_x, center_y), 
                           point1=(center_x, center_y+radius))  #绘制颗粒
    tryTimes += 1                                               #投递失败，尝试次数累加

pickedFaces = f.getSequenceFromMask(mask=('[#1 ]', ), )
p.PartitionFaceBySketch(faces=pickedFaces, sketch=s1)
del model.sketches['__profile__']
session.viewports['Viewport: 1'].setValues(displayedObject=p)

这种颗粒中心点完全随机的方法，对于骨料颗粒含量较少时，是比较适用的，但随着区域被颗粒占据后，尝试的次数会越来越多，就会导致所需的投递时间越长，甚至进入死循环，虽然引入tryTimes进行控制，超过一定尝试次数后终止投递，但无法实现高含量骨料模型的生成。

为了解决上述问题，可以引入区域离散化方法，在每个区域内随机生成一个投递位置，这样让投递区域内相对均匀的布置着可能的颗粒中心，我们只需要对有限的中心位置进行随机选择，并判断他们与已投递颗粒的关系，从而对可能的点进行删减，解决了完全随机过程中的区域重叠问题。同时这种方法还极易普及到非规则几何模型，实现异形构建内投递骨料。POLARIS_MesoConcrete就是采用的该方法，也能生成体素骨料模型。

我们还可以对上述代码进行丰富，半径可以设置为满足一定的分布规律，这里提醒：颗粒尽量从大尺寸到小尺寸的次序进行投递，可以增大颗粒投递成功的概率。

另外也可以增加ITZ层，只需要在刨切的时候，在颗粒中心位置同时生成一个变半径的同心圆即可。

如果需要生成周期性颗粒，首先颗粒的投递区域就不需要边界减去颗粒的半径范围，而是整个区域，其次需要判断颗粒与四条边界的相互关系，如果穿过一条边界，则需要在其相对的边界位置布置一个相同的颗粒。

大家不妨自己试一试，编写自己特色的随机颗粒代码，期待你的分享。

相关文章：
技文|ABAQUS二次开发小工具推荐
插件|POLARIS_PythonTest
插件|POLARIS_MesoConcrete
技文|ABAQUS结果提取大于某值的区域体积
技文|INP关键字跳转、代码高亮、自动补全
技文|Abaqus中提取裂缝数据并用matplotlib库绘图