ABAQUS PYTHON 用35行代码生成二维随机颗粒模型

@[TOC] ABAQUS PYTHON 用35行代码生成二维随机颗粒模型

前面推荐了一些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库绘图

更新内容:1.增加了设置命令流字体背景和边框颜色的功能;2.增加了可保存整个命令流为文本的功能;3.调整了部分右键菜单的位置和快捷键;4.修正了个别情况下打开命令流显示错误的问题;5.修正了一些情况下垂直滚动条出错的问题;6.修正自动滚屏低部留有空白的问题;7.修正不能实时读取命令流的问题。 操作提示: 1、该程序无须放在ABAQUS的工作目录下,可随意放置。另外程序可以随时打开,无须考虑与ABAQUS CAE的打开次序; 1、程序第一次运时需要指定abaqus.rpy的位置,以后运会自动加载上一次设置; 2、在窗口中拖动右键可以移动窗口位置; 3、把鼠标移动到窗口边缘可以拖动改变窗口大小。 该程序主要是给使用ABAQUS的朋友们学习Python用的,可以作为ABAQUS PDE的辅助工具, 对于ABAQUSPython的关系我就不多说了,在ABAQUS CAE中的每一个菜单或按钮操作都是被解释为Python语句,然后才提交上去。 而这些Python语句被适时地保存在工作目录下的abaqus.rpy文件中,这就给我们提供了一个绝好的Python学习途径:进CAE的操作,然后查看abaqus.rpy文件中的对应的Python语句// 该程序会适时的读取abaqus.rpy文件,以便你把相应的CAE操作对照起来// 如果对该程序有什么好的建议或意见,或需要添加什么样的功能, 或者发现什么bug,可以直接给我联系, Email:ck436#126.com(把#改为@). ========================================================== 如果你根本就不能运本程序,那很有可能你还没有安装.NET Framework 2.0以上的平台. .NET Framework是在Microsoft .NET平台上进程序开发和程序运的基础. 给出解决方法: 你可以通过以下几个网址下载: http://www.onlinedown.net/soft/38669.htm http://www.microsoft.com/downloads/details.aspx?FamilyID=0856eacb-4362-4b0d-8edd-aab15c5e04f5&DisplayLang=zh-cn#QuickInfoContainer 如果安装.NET Framework时提示installer错误,则你需要先安装Windows Installer(一般不会遇到): http://www.onlinedown.net/soft/12668.htm http://www.microsoft.com/downloads/details.aspx?displaylang=zh-cn&FamilyID=889482fc-5f56-4a38-b838-de776fd4138c
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值