yade-users team mailing list archive

[Alexander <question269063@xxxxxxxxxxxxxxxxxxxxx>]

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.