Question #269063 on Yade changed:
https://answers.launchpad.net/yade/+question/269063
Status: Answered => Open
Alexander is still having a problem:
Hi, Jerome
Thank’s a lot for such a detail answer! So I’ll try to be more precise))
So in the first picture (without spheres) I just draw the orthogonal plate
like a real solid body to show what type of problem I have (without any
connection to solving it with DEM).
As I said before boundary conditions are following
1. On face ABCD applied tensile force F parallel to OY axis (red arrow).
it is necessary to say that force F is set for the whole face and not to
specific points.
2. Face KGHO is fixed (Red dots)
So it means that the plate will stretch after force F has been applied to the
side ABCD, because of the opposite side KGHO is fixed.
Objective is to compute distribution of strain, stress and displacement values
on the plate after force has been applied. The final goal of my research is to
solve this problem with 2 methods DEM and FEM and compare results.
In DEM I image it to do like this:
1) Represent plate like a set of spheres
2) Find correct way to set boundary conditions and apply them on spheres
3) Find strain, stress and displacement values for each sphere which
represent the plate and save result data to file. (So now I save each sphere
and it’s strain, stress, displacement attributes to VTK file which can be
analyzed in Paraview)
Below I posted MWE example as you name it)) where permanent force is applied
for spheres connected to boundary ABCD, spheres connected to KGHO is set to be
fixed. (Thank’s for advising regularOrtho() function but still use simple
loops:)
Yade version: 2015-03-17.git-16ecad0
////
###################################################
# define materials and model configuration
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'))
# represent plate like a set of regular monosized set of spheres
# also set boundary conditions via predefined tensile force for spheres on ABCD
and
# fixed spheres on KGHO
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))
spheres.append(O.bodies[id])
if j == 15:
O.forces.addF(id,(0,p,0),permanent=True) # Add force for all spheres
connected to ABCD
if j == 0:
spheres[id].state.blockedDOFs='y' # Fixed all spheres connected to
KGHO
###################################################
# define engines
# 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)
]
###################################################
# 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(100,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')
# 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')
////
This code should give the results shown in the pictures:
displacement - http://s28.postimg.org/ooujic9rh/puc1.png
stress - http://s29.postimg.org/dw4qu3xw7/puc1_stress.png
strain - http://s13.postimg.org/6bbtsjfav/puc1_strain.png
These images shows set of spheres representing the plate (variable “spheres” in
code) where each sphere is drawn with color depending on value of one chosen
attribute (strain, stress or disp)
(as I said above each sphere contains values of 3 attributes strain, stress and
disp)
So I don’t like stress and strain distribution, because for example
spheres connected to ABCD and KGHO have minimal stress values.
And you already said my questions yourself :)
"I used this approach, why are the results so strange?"
“what approach is the best for this simulation?”
Also I wanna say that it’s very important for me to solve this problem:)
Soon I will post a picture with FEM results to compare, because in our
test we consider FEM like main method and verify DEM via FEM.
I’ll be very glad if you can help me, thank’s a lot again
with regards, Alexander
--
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 : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
