本文已參與「新人創作禮」活動,一起開啟掘金創作之路.

前言

According to the problems existing in the previous article（random fiber generation、Section assignment、切割）organize learning

代碼

Guide package and parameters

from abaqus import *
from abaqusConstants import *
from caeModules import *
from driverUtils import executeOnCaeStartup
from math import *
import random as rd
# Radius of the fiberrDepth with base length
r = 2
w = 50
depth = 5
# List of recorded fiber center coordinatesxy
xy = [[rd.random()*width, rd.random()*width]]
# Minimum distance between two fibers,Prevent the grid from being too small if the distance is too close
l_min = 2*r+r/4
# Target Fiber Volume Fraction vs Current Fiber Volume Fraction
# The target volume fraction should not be too large,理論上<0.5或<0.6基本不會出錯
vf = 0.3
vf_curr = 0
# 材料參數（Modulus poisson's ratio）
E_m = 50
v_m = 0.3
E_f = 273
v_f = 0.2
復制代碼

Create a matrix and give cross section

s = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=200.0)
s.rectangle(point1=(0.0, 0.0), point2=(w, w))
p = mdb.models['Model-1'].Part(name='matrix', dimensionality=THREE_D, type=DEFORMABLE_BODY)
p.BaseSolidExtrude(sketch=s, depth=depth)
del mdb.models['Model-1'].sketches['__profile__']
mdb.models['Model-1'].Material(name='Material-matrix')
mdb.models['Model-1'].materials['Material-matrix'].Elastic(table=((E_m, v_m),))
mdb.models['Model-1'].HomogeneousSolidSection(name='Section-matrix', material='Material-matrix', thickness=None)
c = p.cells
cells = c.findAt(((w, w, 0),))
region = regionToolset.Region(cells=cells)
p.SectionAssignment(region=region, sectionName='Section-matrix', offset=0.0, offsetType=MIDDLE_SURFACE, offsetField='', thicknessAssignment=FROM_SECTION)
復制代碼

創建纖維

According to your different needs,Fiber coordinate generation code can be modified,Wrote a judging condition here is not very simple,待優化

def judge(cir):
for curr_cir in xy:
if (cir[0]-curr_cir[0])**2+(cir[1]-curr_cir[1])**2 < l_min**2:
return False
else:
continue
return True
def add(cir):
if judge(cir):
# If the fiber and the four sides do not overlap,則直接添加
if cir[0] > r+r/2 and w-cir[0] > r+r/2 and cir[1] > r+r/2 and w-cir[1] > r+r/2:
xy.append(cir)
else:
# Fibers and quadrilateral overlap
# 其中距離4no corner too close
if cir[0]**2+cir[1]**2 > 2*r**2 and cir[0]**2+(cir[1]-w)**2 > 2*r**2 and (cir[0]-w)**2+cir[1]**2 > 2*r**2 and (cir[0]-w)**2+(cir[1]-w)**2 > 2*r**2:
# Determine whether the fibers on the four sides need to be periodic,並添加
if cir[0] < r:
if r-cir[0] > 0.8:
if judge([cir[0]+w, cir[1]]):
xy.append(cir)
xy.append([cir[0]+w, cir[1]])
if w-cir[0] < r:
if r-(w-cir[0]) > 0.8:
if judge([cir[0]-w, cir[1]]):
xy.append(cir)
xy.append([cir[0]-w, cir[1]])
if cir[1] < r:
if r-cir[1] > 0.8:
if judge([cir[0], cir[1]+w]):
xy.append(cir)
xy.append([cir[0], cir[1]+w])
if w-cir[1] < r:
if r-(w-cir[1]) > 0.8:
if judge([cir[0], cir[1]-w]):
xy.append(cir)
xy.append([cir[0], cir[1]-w])
while vf_curr < vf:
cir = [rd.random()*w, rd.random()*w]
add(cir)
vf_curr = len(xy)*math.pi*r**2/w**2
print(vf_curr)
# The center coordinates of the fibers created above are all inxy列表裡了
# 在s1circle in sketch,extrude into solid
s1 = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=200.0)
for cir in xy:
s1.CircleByCenterPerimeter(center=(cir[0], cir[1]), point1=(cir[0]+r, cir[1]))
p = mdb.models['Model-1'].Part(name='fiber', dimensionality=THREE_D,
type=DEFORMABLE_BODY)
p.BaseSolidExtrude(sketch=s1, depth=depth)
del mdb.models['Model-1'].sketches['__profile__']
復制代碼

Give the fiber a cross section

mdb.models['Model-1'].Material(name='Material-fiber')
mdb.models['Model-1'].materials['Material-fiber'].Elastic(table=((E_f, v_f),))
mdb.models['Model-1'].HomogeneousSolidSection(name='Section-fiber', material='Material-fiber', thickness=None)
p = mdb.models['Model-1'].parts['fiber']
c = p.cells
for cir in xy:
cells = c.findAt(((cir[0], cir[1], 0.0),))
region = regionToolset.Region(cells=cells)
p.SectionAssignment(region=region, sectionName='Section-fiber', offset=0.0, offsetType=MIDDLE_SURFACE, offsetField='', thicknessAssignment=FROM_SECTION)
復制代碼

裝配

a = mdb.models['Model-1'].rootAssembly
a.DatumCsysByDefault(CARTESIAN)
p = mdb.models['Model-1'].parts['fiber']
a.Instance(name='fiber-1', part=p, dependent=OFF)
p = mdb.models['Model-1'].parts['matrix']
a.Instance(name='matrix-1', part=p, dependent=OFF)
a.InstanceFromBooleanMerge(name='UDcomp', instances=(a.instances['fiber-1'], a.instances['matrix-1'], ), keepIntersections=ON, originalInstances=SUPPRESS, domain=GEOMETRY)
復制代碼

切割

該文章Various methods of selecting geometric elements are introduced,這裡使用findAt()函數

p = mdb.models['Model-1'].parts['UDcomp']
f, e = p.faces, p.edges
# 通過findAtFind canvas position
face = f.findAt(((rd.random()*w, rd.random()*w, depth),))
edge = e.findAt(((w, rd.random()*w, depth),))
# face和edge為GeomSequence object,There is only one element in itface object、edge object
# 可通過face[0].getCentroid()獲取面、barycentric coordinates of edges,Find face to observe、Whether the edge meets the requirements
t = p.MakeSketchTransform(sketchPlane=face[0], sketchUpEdge=edge[0],
sketchPlaneSide=SIDE1, sketchOrientation=RIGHT, origin=(w/2, w/2,
depth))
s = mdb.models['Model-1'].ConstrainedSketch(name='__profile__',
sheetSize=200, transform=t)
p.projectReferencesOntoSketch(sketch=s, filter=COPLANAR_EDGES)
s.rectangle(point1=(-w/2, -w/2), point2=(w/2, w/2))
s.rectangle(point1=(-w/2-2*r, -w/2-2*r), point2=(w/2+2*r, w/2+2*r))
p.CutExtrude(sketchPlane=face[0], sketchUpEdge=edge[0], sketchPlaneSide=SIDE1,
sketchOrientation=RIGHT, sketch=s, flipExtrudeDirection=OFF)
del mdb.models['Model-1'].sketches['__profile__']
復制代碼

mesh

p = mdb.models['Model-1'].parts['UDcomp']
p.seedPart(size=1.0, deviationFactor=0.1, minSizeFactor=0.1)
c = p.cells
pickedRegions = c.getSequenceFromMask(mask=('[#ffffffff #3fffffff ]', ), )
p.setMeshControls(regions=pickedRegions, elemShape=WEDGE)
p.generateMesh()
復制代碼

EasyPBC

## 依然是,Need to remove redundant instances,否則報錯
a = mdb.models['Model-1'].rootAssembly
a.deleteFeatures(('fiber-1', 'matrix-1', ))
import sys
sys.path.insert(8, r'e:/CAE/abaqus_plugins/EasyPBC V.1.4')
import easypbc
easypbc.feasypbc(part='Model-1', inst='UDcomp-1', meshsens=1E-07, CPU=8,
E11=True, E22=True, E33=False, G12=True, G13=False, G23=False,
onlyPBC=False, CTE=False, intemp=0, fntemp=100)
復制代碼

效果

輸出的結果