Question #269063 on Yade changed:
https://answers.launchpad.net/yade/+question/269063
Description changed to:
Hello, everybody.
I’m new in Yade and I need help to simulate following situation.
A) Given: metallic rectangular plate with sides parallel or
perpendicular to coordinate axes (the plate may consist 2 differ materials)
such as in linked picture.
B) Boundary conditions:
1. On face ABCD applied tensile force F parallel
to OY axis
2. Face KGHO is fixed.
C) Objective is: to compute distribution of strain, stress and
displacement values on the plate.
As I understand in DEM it can be done through computing these values for each
sphere in the plate.
I tried to use:
- UniaxialStrainer but it doesn’t allow to apply force (or stress) just
constant strain.
- TriaxialStressController but it works only for compression not for
tension (or may be I do something wrong)
- Also i set permanent force for spheres on the boundary ABCD (but I
think it’s not correct because it’s necessary to set force like pressure for
every point on the surface and size of spheres may be different, or may be I’m
wrong again)
- In addition, I used constant velocity for boundary spheres but it’s
just for testing and doesn’t deserve any attention.
I compute stress and strain with
TesselationWrapper().computeDeformations() and bodyStressTensors()
functions, displacement with displ(), then save values manually to VTK
file. I didn't use VTKExporter() because at a later stage I'll need
another file extension to save the result data.
Now I represent the plate as a set of uniform arrows of spheres with
equal radius connected one by one:
spheres=[]
for i in range(0, 16):
for j in range(0, 16):
for k in range(0, 2):
id =
O.bodies.append(sphere([i+0.5,j+0.5,k+0.5],material='mat_mod',radius=r))
Below is a code of executing file “example.py”, where I tried to test 4
variants described above to set boundary conditions (parameter “mode”).
-----------------------------------------------------------------------------------------------------------------------------------------------
###################################################
# define materials and model configuration
mode = 3 # 1 - constant velocity applied to boundary spheres, 2 -
constant force ..., 3 - uniaxial tension, 4 - triaxial test
E = 1161e9 # Young's modulus of model
v = 0.33 # Poisson's ratio
p = -150e6 # initial stress value (positive - tension, negative - compression)
e = 0.02 # loading rate (strain rate) (positive - tension, negative -
compression)
r = 0.5 # spheres radius
# Enlarge interaction radius between spheres using "interaction_radius"
parameter (for example in uniax.py this value is 1.5)
interaction_radius = 1.5
# define plate material, create "dense" packing by setting friction to zero
initially
O.materials.append(CpmMat(young=E,
frictionAngle=0,
poisson=v,
density=4800,
sigmaT=3.5e6,
epsCrackOnset=1e-4,
isoPrestress=0,
relDuctility=30,
label = 'mat_mod'))
# define constraints material
O.materials.append(CpmMat(young=E,
frictionAngle=0,
poisson=v,
density=4800,
sigmaT=3.5e6,
epsCrackOnset=1e-4,
isoPrestress=0,
relDuctility=30,
label = 'mat_con'))
# define walls material
O.materials.append(CpmMat(young=E,
frictionAngle=0,
poisson=v,
density=0,
sigmaT=3.5e6,
epsCrackOnset=1e-4,
isoPrestress=0,
relDuctility=30,
label = 'mat_wal'))
# create walls around the packing if necessary (must be used before spheres
have been added!)
if mode==4:
mn,mx=Vector3(0,0,0),Vector3(16,16,2) # corners of the initial packing
walls=aabbWalls([mn,mx],thickness=0,material='mat_wal')
wallIds=O.bodies.append(walls)
# pack box with spheres
spheres=[]
for i in range(0, 16):
for j in range(0, 16):
for k in range(0, 2):
# for spheres which position satisfies condition below set constraint
material
if (i == 4 or i == 5 or i == 10 or i == 11) and (j == 4 or j == 5 or j
== 10 or j == 11):
id =
O.bodies.append(sphere([i+0.5,j+0.5,k+0.5],material='mat_mod',radius=r))
spheres.append(O.bodies[id])
else:
id =
O.bodies.append(sphere([i+0.5,j+0.5,k+0.5],material='mat_mod',radius=r))
spheres.append(O.bodies[id])
###################################################
# define engines
# mode = 1
if mode == 1:
# define function to keep permanent velocity of boundary spheres during
simulation..
vel_list=[]
vel_list2=[]
def updateStrain():
for b in vel_list:
b.state.vel=(0,e,0)
for b in vel_list2: # DELETE
b.state.vel=(0,-e,0) # DELETE
# ..and also set start velocity for them
for b in spheres:
if b.state.pos[1] == 0.5:
#b.state.blockedDOFs='y'
b.state.vel=(0,-e,0) # DELETE
vel_list2.append(b) # DELETE
if b.state.pos[1] == 15.5:
b.state.vel=(0,e,0)
vel_list.append(b)
# simulation loop
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interaction_radius)]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interaction_radius)],
[Ip2_CpmMat_CpmMat_CpmPhys()],
[Law2_ScGeom_CpmPhys_Cpm()]
),
PyRunner(iterPeriod=1,command='updateStrain()'),
NewtonIntegrator(damping=0.4)
]
# mode = 2
if mode == 2:
# set initial force for boundary spheres
for b in spheres:
if b.state.pos[1] == 15.5:
O.forces.addF(b.id,(0,p,0),permanent=True)
if b.state.pos[1] == 0.5:
b.state.blockedDOFs='y'
# simulation loop
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interaction_radius)]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interaction_radius)],
[Ip2_CpmMat_CpmMat_CpmPhys()],
[Law2_ScGeom_CpmPhys_Cpm()]
),
CpmStateUpdater(realPeriod=1),
NewtonIntegrator(damping=0.4)
]
# mode = 3
if mode == 3:
# get basic settings for UniaxialStrainer
bb=utils.uniaxialTestFeatures()
negIds,posIds,axle,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
# define UniaxialStrainer main parameters
uniax=UniaxialStrainer(
strainRate=e,
axis=axle,
asymmetry=0,
posIds=posIds,
negIds=negIds,
crossSectionArea=crossSectionArea,
blockDisplacements=False,
blockRotations=False,
setSpeeds=True,
stressUpdateInterval=1
)
# simulation loop
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interaction_radius)]),#,sweepLength=.05*r),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interaction_radius)],
[Ip2_CpmMat_CpmMat_CpmPhys()],
[Law2_ScGeom_CpmPhys_Cpm()]
),
NewtonIntegrator(damping=0.4),
CpmStateUpdater(realPeriod=1),
uniax
]
# mode = 4
if mode == 4:
# define TriaxialStressController main parameters
triax=TriaxialStressController(
maxMultiplier=1.+2e4/E, # spheres growing factor (fast growth)
finalMaxMultiplier=1.+2e3/E, # spheres growing factor (slow growth)
thickness = 0,
stressMask = 7, # (0 - stress mode, 7 - strain mode)
internalCompaction=False, # If true the confining pressure is
generated by growing particles !!!
goal1 = 0,
goal2 = p, # positive is tension, negative is compression
goal3 = 0,
wall_bottom_activated=True, # this wall is like upper boundary along
OY axis
wall_top_activated=True,
wall_left_activated=True,
wall_right_activated=True,
wall_back_activated=True,
wall_front_activated=True
)
# simulation loop
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interaction_radius),Bo1_Box_Aabb()]),#,sweepLength=.05*r),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interaction_radius),Ig2_Box_Sphere_ScGeom()],
[Ip2_CpmMat_CpmMat_CpmPhys()],
[Law2_ScGeom_CpmPhys_Cpm()]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
triax,
NewtonIntegrator(damping=0.2)
]
###################################################
# start simulation and compute strain and stress
# try to run script with qt graphical interface
try:
yade.qt.Controller(), yade.qt.View()
except:
print 'Qt graphical interface is not avaliable'
# set the integration timestep to be 1/2 of the "critical" timestep
O.dt=.5*utils.PWaveTimeStep()
# compute strain via tesselation wrapper.
TW=TesselationWrapper()
# store current positions before simulation
TW.setState(0)
# run simulation
O.run(10000,1)
# store positions after simulation (deformed state)
TW.setState(1)
# compute deformation for each body
TW.computeDeformations()
# compute stress tensor for each body
stresses = bodyStressTensors()
###################################################
# save data to vtk. file
n = len(spheres) # get number of particles
f = open('result.vtk','w')
# header
f.write('# vtk DataFile Version 3.0.\n')
f.write('comment\n')
f.write('ASCII\n\n')
# center position
string = str(n)
f.write('DATASET POLYDATA\n')
f.write('POINTS ')
f.write(string )
f.write(' double\n')
for b in spheres:
for i in range(0,3):
value = b.state.pos[i]
string = str(value)
f.write(string)
f.write(' ')
f.write('\n')
f.write('\n')
# radius
string = str(n)
f.write('POINT_DATA ')
f.write(string)
f.write('\n')
f.write('SCALARS radius double 1\n')
f.write('LOOKUP_TABLE default\n')
for b in spheres:
value = b.shape.radius
string = str(value)
f.write(string)
f.write('\n')
f.write('\n')
# velocity
f.write('VECTORS velocity double\n')
for b in spheres:
for i in range(0,3):
value = b.state.vel[i]
string = str(value)
f.write(string)
f.write(' ')
f.write('\n')
f.write('\n')
# strain
f.write('TENSORS strain_tensor float\n')
for b in spheres:
for i in range(0,3):
for j in range(0,3):
value = TW.deformation(b.id,i,j)
string = str(value)
f.write(string)
f.write(' ')
f.write('\n')
f.write('\n')
f.write('\n')
# stress
f.write('TENSORS stress_tensor float\n')
for b in spheres:
for i in range(0,3):
for j in range(0,3):
value = stresses[b.id][i,j]
string = str(value)
f.write(string)
f.write(' ')
f.write('\n')
f.write('\n')
f.write('\n')
# displacement
f.write('VECTORS displacement double\n')
for b in spheres:
for i in range(0,3):
value = b.state.displ()[i]
string = str(value)
f.write(string)
f.write(' ')
f.write('\n')
f.write('\n')
# material - young modulus
f.write('SCALARS young double 1\n')
f.write('LOOKUP_TABLE default\n')
for b in spheres:
value = b.material.young
string = str(value)
f.write(string)
f.write('\n')
f.write('\n')
-----------------------------------------------------------------------------------------------------------------------------------------------
I’ll be glad any help
with regards, Alexander
Picture of the plate: http://s18.postimg.org/omoiobfcp/pic.jpg
--
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
