星辰技文|一步步教你如何使用35行代碼生成ABAQUS二維隨機顆粒模型

星辰_北極星

2023年7月13日 17:37

瀏覽：3885 評論：6 收藏：32

前面推薦了一些ABAQUS二次開發(fā)小工具，不知道大家是否已經(jīng)安裝使用。

后面以一些小案例帶大家熟悉ABAQUS前后處理相關(guān)的Python庫，以及使用技巧。

星哥開發(fā)的插件大多集中在非均質(zhì)相關(guān)斷裂問題，相信關(guān)注公眾號的很多朋友也都是做這方面，那么我們從最初始的非均質(zhì)幾何模型的案例出發(fā)，來演示一個隨機顆粒模型的代碼編寫的全過程，效果如下所示：

星辰技文|一步步教你如何使用35行代碼生成ABAQUS二維隨機顆粒模型的圖2

在這個案例中，最大的幫手是PythonReader,它能讓初學者能輕松了解GUI界面中的每個操作對應的代碼段是什么，比如，點擊新建文件按鈕，會彈出以下代碼：

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

其中第一行就是創(chuàng)建新數(shù)據(jù)模型的命令，第二行和第三行均為注釋，描述模型的信息，第4行則是視圖設置當前顯示對象為None，即空。大多數(shù)的前處理操作均可以用這種方式進行錄制，只需要了解一些Python基礎知識和接口代碼的風格和結(jié)構(gòu)，小白也能輕松上手。

為實現(xiàn)隨機顆粒模型的代碼編寫，我們分3步進行：

? 首先通過錄制獲得ABAQUS中建模的相關(guān)代碼。
? 然后修改相應代碼，刪除無用代碼，實現(xiàn)代碼的參數(shù)化。
? 最后在代碼中添加循環(huán)和干涉判斷，實現(xiàn)多個顆粒的隨機投遞。

第一步：錄制

界面中操作：

1）新建部件，Modeling Space定義為2D Planar，其它保持默認，進入草圖

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

2）草圖內(nèi)繪制一個矩形，輸入坐標[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）完成草圖，實現(xiàn)顆粒的繪制

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

第二步：參數(shù)化

有一定英語基礎，并了解基本的Python類語法的概念，我們就能輕松讀懂每一行代碼，從而了解哪些代碼可以被簡化或刪除：

from abaqus import *
from abaqusConstants import *
# 創(chuàng)建草圖，保留
s = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=200.0)
# 定義的變量g,v,d,c均未調(diào)用，該行可刪除
g, v, d, c = s.geometry, s.vertices, s.dimensions, s.constraints
# 窗口顯示草圖，刪除不影響
s.setPrimaryObject(option=STANDALONE)
# 創(chuàng)建矩形，保留
s.rectangle(point1=(0.0, 0.0), point2=(50.0, 50.0))
# 創(chuàng)建部件，保留
p = mdb.models['Model-1'].Part(name='Part-1', dimensionality=TWO_D_PLANAR, type=DEFORMABLE_BODY)
# p變量已經(jīng)定義，并指向同一個對象，該行可刪除
p = mdb.models['Model-1'].parts['Part-1']
# 部件添加基礎面，保留
p.BaseShell(sketch=s)
# 窗口隱藏草圖，與顯示行對應，刪除不影響
s.unsetPrimaryObject()
# p變量已經(jīng)定義，并指向同一個對象，該行可刪除
p = mdb.models['Model-1'].parts['Part-1']
# 視圖顯示部件，可移到最后
session.viewports['Viewport: 1'].setValues(displayedObject=p)
# 刪除臨時草圖，保留
del mdb.models['Model-1'].sketches['__profile__']
# p變量已經(jīng)定義，并指向同一個對象，該行可刪除
p = mdb.models['Model-1'].parts['Part-1']
# 定義的變量e,d1均未調(diào)用，該行可簡化
f, e, d1 = p.faces, p.edges, p.datums
# 創(chuàng)建草圖投影面，默認草圖中心會在f[0]的幾何中心，需與全局坐標一致，origin可修改為0，0，保留
t = p.MakeSketchTransform(sketchPlane=f[0], sketchPlaneSide=SIDE1, origin=(
    25.0, 25.0, 0.0))
# 創(chuàng)建草圖
s1 = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', 
    sheetSize=141.42, gridSpacing=3.53, transform=t)
# 定義的變量g,v,d,c均未調(diào)用，該行可刪除
g, v, d, c = s1.geometry, s1.vertices, s1.dimensions, s1.constraints
# 窗口顯示草圖，刪除不影響
s1.setPrimaryObject(option=SUPERIMPOSE)
# p變量已經(jīng)定義，并指向同一個對象，該行可刪除
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變量已經(jīng)定義，并指向同一個對象，該行可刪除
p = mdb.models['Model-1'].parts['Part-1']
# 創(chuàng)建面對象，該變量已經(jīng)定義，且指向的對象不變，可刪除
f = p.faces
# 選擇刨切的面對象pickedFaces，保留
pickedFaces = f.getSequenceFromMask(mask=('[#1 ]', ), )
# 定義的變量e1, d2均未調(diào)用，該行可刪除
e1, d2 = p.edges, p.datums
# 按草圖s1刨切選擇的面對象pickedFaces，保留
p.PartitionFaceBySketch(faces=pickedFaces, sketch=s1)
# 窗口隱藏草圖，與顯示行對應，刪除不影響
s1.unsetPrimaryObject()
# 刪除臨時草圖，保留
del mdb.models['Model-1'].sketches['__profile__']

同時代碼中大量出現(xiàn)“mdb.models['Model-1']”，可創(chuàng)建變量model進行代替，簡化代碼，清理后的代碼：

from abaqus import *
from abaqusConstants import *
# 創(chuàng)建model變量
model = mdb.models['Model-1']
# 創(chuàng)建草圖，保留
s = model.ConstrainedSketch(name='__profile__', sheetSize=200.0)
# 創(chuàng)建矩形，保留
s.rectangle(point1=(0.0, 0.0), point2=(50.0, 50.0))
# 創(chuàng)建部件，保留
p = model.Part(name='Part-1', dimensionality=TWO_D_PLANAR, type=DEFORMABLE_BODY)
# 部件添加基礎面，保留
p.BaseShell(sketch=s)
# 刪除臨時草圖，保留
del model.sketches['__profile__']
# 定義的變量e,d1均未調(diào)用，該行可簡化
f = p.faces
# 創(chuàng)建草圖投影面，默認草圖中心會在f[0]的幾何中心，需與全局坐標一致，origin可修改為0，0，保留
t = p.MakeSketchTransform(sketchPlane=f[0], sketchPlaneSide=SIDE1, origin=(0, 0, 0))
# 創(chuàng)建草圖
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)

模型中可參數(shù)化的變量包含:

? 部件名稱 partName
? 部件寬度 width
? 部件高度 height
? 顆粒中心 center_x,center_y
? 顆粒半徑 radius

因此單顆粒的參數(shù)化代碼如下：

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)

第三步：循環(huán)隨機

在上面單顆粒參數(shù)化模型基礎上，添加循環(huán)并隨機生成顆粒中心坐標，即可實現(xiàn)隨機多顆粒的生成。為了避免顆粒之間相互重疊，需要將新生成的坐標與原有坐標進行距離判斷，二者距離需要大于2倍圓半徑，具體代碼如下：

import random
from abaqus import *
from abaqusConstants import *
partName = "Part-1"     #部件名稱
width = 100             #部件寬度
height=50               #部件高度
radius = 5              #顆粒半徑
rockNum = 20            #顆粒數(shù)量
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            #引入嘗試次數(shù)，避免進入死循環(huán)
rockCenterList = []     #用于存儲已經(jīng)投遞成功的顆粒中心
while tryTimes<1000 and len(rockCenterList)<rockNum:
    #在box范圍內(nèi)隨機生成顆粒中心坐標
    center_x = random.uniform(radius,width-radius)
    center_y = random.uniform(radius,height-radius)
    #循環(huán)判斷與其它顆粒距離是否小于目標值
    mark = 1    #用于記錄是否符合要求
    for x,y in rockCenterList:
        dist = sqrt((x-center_x)**2+(y-center_y)**2)
        # 距離判斷，一個距離小于2倍半徑，則直接退出該層循環(huán)
        if dist<2*radius:
            mark = 0
            break
    if mark == 1:
        tryTimes = 0                                            #投遞成功，重新計數(shù)
        rockCenterList.append([center_x,center_y])              #更新rockCenterList
        s1.CircleByCenterPerimeter(center=(center_x, center_y), 
                           point1=(center_x, center_y+radius))  #繪制顆粒
    tryTimes += 1                                               #投遞失敗，嘗試次數(shù)累加

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

這種顆粒中心點完全隨機的方法，對于骨料顆粒含量較少時，是比較適用的，但隨著區(qū)域被顆粒占據(jù)后，嘗試的次數(shù)會越來越多，就會導致所需的投遞時間越長，甚至進入死循環(huán)，雖然引入tryTimes進行控制，超過一定嘗試次數(shù)后終止投遞，但無法實現(xiàn)高含量骨料模型的生成。

為了解決上述問題，可以引入區(qū)域離散化方法，在每個區(qū)域內(nèi)隨機生成一個投遞位置，這樣讓投遞區(qū)域內(nèi)相對均勻的布置著可能的顆粒中心，我們只需要對有限的中心位置進行隨機選擇，并判斷他們與已投遞顆粒的關(guān)系，從而對可能的點進行刪減，解決了完全隨機過程中的區(qū)域重疊問題。同時這種方法還極易普及到非規(guī)則幾何模型，實現(xiàn)異形構(gòu)建內(nèi)投遞骨料。POLARIS_MesoConcrete就是采用的該方法，也能生成體素骨料模型。

星辰技文|一步步教你如何使用35行代碼生成ABAQUS二維隨機顆粒模型的圖3

星辰技文|一步步教你如何使用35行代碼生成ABAQUS二維隨機顆粒模型的圖4

我們還可以對上述代碼進行豐富，半徑可以設置為滿足一定的分布規(guī)律，這里提醒：顆粒盡量從大尺寸到小尺寸的次序進行投遞，可以增大顆粒投遞成功的概率。

另外也可以增加ITZ層，只需要在刨切的時候，在顆粒中心位置同時生成一個變半徑的同心圓即可。

星辰技文|一步步教你如何使用35行代碼生成ABAQUS二維隨機顆粒模型的圖5

如果需要生成周期性顆粒，首先顆粒的投遞區(qū)域就不需要邊界減去顆粒的半徑范圍，而是整個區(qū)域，其次需要判斷顆粒與四條邊界的相互關(guān)系，如果穿過一條邊界，則需要在其相對的邊界位置布置一個相同的顆粒。

星辰技文|一步步教你如何使用35行代碼生成ABAQUS二維隨機顆粒模型的圖6

大家不妨自己試一試，編寫自己特色的隨機顆粒代碼，期待你的分享。

相關(guān)文章：
技文|ABAQUS二次開發(fā)小工具推薦
 插件|POLARIS_PythonTest
插件|POLARIS_MesoConcrete
技文|ABAQUS結(jié)果提取大于某值的區(qū)域體積
 技文|INP關(guān)鍵字跳轉(zhuǎn)、代碼高亮、自動補全
 技文|Abaqus中提取裂縫數(shù)據(jù)并用matplotlib庫繪圖