A friend asked me to write a python script to control the FEM software Abaqus. He wanted to simulate a compound material consisting of a material block with small spheres of a diffrent material inserted. It was the first time I ever saw Abaqus, but got the hang of it quite fast.
The resulting script produces a block of material, generates n spheres wich are distributed randomly and not overlapping inside the block. The radius of the spheres is controled by their number. It is choosen in a way that the n spheres make up 20% of the volume of the block.
Sphere Positions and radii are written out into a file. Im looking forward to his simulation results.
from part import *
from material import *
from section import *
from assembly import *
from step import *
from interaction import *
from load import *
from mesh import *
from job import *
from sketch import *
from visualization import *
from connectorBehavior import *
import random
import math
#Euclidean Distance Measure
def euclidean_py(x, y):
# lightly modified from implementation by Thomas Sicheritz-Ponten.
# This works faster than the Numeric implementation on shorter
# vectors.
if len(x) != len(y):
raise ValueError, "vectors must be same length"
sum = 0
for i in range(len(x)):
sum += (x[i]-y[i])**2
return math.sqrt(sum)
#Create Cube
mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=200.0)
mdb.models['Model-1'].sketches['__profile__'].rectangle(point1=(0.0, 0.0), point2=(1.0, 1.0))
mdb.models['Model-1'].Part(dimensionality=THREE_D, name='Part-1', type=DEFORMABLE_BODY)
mdb.models['Model-1'].parts['Part-1'].BaseSolidExtrude(depth=1.0, sketch=mdb.models['Model-1'].sketches['__profile__'])
del mdb.models['Model-1'].sketches['__profile__']
#Instantiate cube
mdb.models['Model-1'].rootAssembly.Instance(dependent=ON, name='tmpCube', part=mdb.models['Model-1'].parts['Part-1'])
#File Output
text_file = open("c:\Sphere_Positions.txt", "w")
#Number of Spheres
n = 10
text_file.write("Number of Spheres: " + str(n) + "\n")
#Radius of one sphere
r = (3*0.20)/(4*math.pi*n)
r = r**(1.0/3.0)
print "Sphere Radius is " + str(r)
text_file.write("Sphere Radius is " + str(r) + "\n")
#Create Sphere
mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=200.0)
mdb.models['Model-1'].sketches['__profile__'].ConstructionLine(point1=(0.0, -100.0), point2=(0.0, 100.0))
mdb.models['Model-1'].sketches['__profile__'].FixedConstraint(entity=mdb.models['Model-1'].sketches['__profile__'].geometry[2])
mdb.models['Model-1'].sketches['__profile__'].ArcByCenterEnds(center=(0.0, 0.0), direction=CLOCKWISE, point1=(0.0, -r), point2=(0.0, r))
mdb.models['Model-1'].sketches['__profile__'].Line(point1=(0.0, -r), point2=(0.0, r))
mdb.models['Model-1'].sketches['__profile__'].VerticalConstraint(entity=mdb.models['Model-1'].sketches['__profile__'].geometry[4])
mdb.models['Model-1'].Part(dimensionality=THREE_D, name='Part-2', type=DEFORMABLE_BODY)
mdb.models['Model-1'].parts['Part-2'].BaseSolidRevolve(angle=360.0, flipRevolveDirection=OFF, sketch=mdb.models['Model-1'].sketches['__profile__'])
del mdb.models['Model-1'].sketches['__profile__']
#List of Spheres
sphereList = []
sphereInstancesList = []
#Create n instances of the sphere
for i in range(1, n+1):
InstanceName = 'Sphere_' + str(i)
print InstanceName
text_file.write(InstanceName)
#Maximum tries to distribute sphere
maxTries = 1000
while len(sphereList) < i:
maxTries -= 1
if maxTries < 1:
print "Maximum Distribution tries exceded. Error! Restart the Script!"
break;
#Make sure Spheres dont cut cube sides
vecPosition = [r+(random.random()*(1.0-r-r)),r+(random.random()*(1.0-r-r)),r+(random.random()*(1.0-r-r))]
for pos in sphereList:
if euclidean_py(pos, vecPosition) < 2*r:
break
else:
sphereList.append(vecPosition)
print vecPosition
text_file.write("\t" + str(vecPosition) + "\n")
#Instantiate Sphere
mdb.models['Model-1'].rootAssembly.Instance(dependent=ON, name=InstanceName, part=mdb.models['Model-1'].parts['Part-2'])
#Translate Instance of Sphere
mdb.models['Model-1'].rootAssembly.translate(instanceList=(InstanceName, ), vector=vecPosition)
sphereInstancesList.append(mdb.models['Model-1'].rootAssembly.instances[InstanceName])
print "Sphere Radius is " + str(r)
#Cut all spheres with cube
mdb.models['Model-1'].rootAssembly.PartFromBooleanCut(cuttingInstances=(
sphereInstancesList ),
instanceToBeCut=mdb.models['Model-1'].rootAssembly.instances['tmpCube'],
name='Part-3')
#Intantiate Cut Object
mdb.models['Model-1'].rootAssembly.Instance(dependent=ON, name='Cube',
part=mdb.models['Model-1'].parts['Part-3'])
#Delete temporary cube
del mdb.models['Model-1'].rootAssembly.instances['tmpCube']
text_file.close()
Filed under: Abaqus, Python
Word Gems 