yade-users team mailing list archive

[Alexander <question269063@xxxxxxxxxxxxxxxxxxxxx>]

Question #269063 on Yade changed:
https://answers.launchpad.net/yade/+question/269063

Alexander gave more information on the question:
Also code reproduces last picture

import __builtin__
from yade import export

###################################################
# define materials and model configuration

E = 1.338e11#2e11 # Young's modulus of model
v = 0.3#0.3 # Poisson's ratio
p = 150e6 # initial stress value (positive - tension, negative - compression)
d = 7850 # density

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=d,
						  sigmaT=3.5e6,
						  epsCrackOnset=1e-4,
						  neverDamage=True,
						  isoPrestress=0,
						  relDuctility=30,
						  label = 'mat'))
	
# 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',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='xyzXYZ' # Fixed all spheres connected to KGHO 
           O.forces.addF(id,(0,-p,0),permanent=True)		   
		
      
###################################################
# define engines

# simulation loop 
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interaction_radius,label='bo1s')]),
 InteractionLoop(
	[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interaction_radius,label='ig2ss')],
    [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 one single step
O.step()

# reset interaction radius to the default value
bo1s.aabbEnlargeFactor=1.0
ig2ss.interactionDetectionFactor=1.0

# run simulation, until static equilibrium will not reached
while unbalancedForce()>1e-2:
   O.run(10,True)

# store positions after simulation (deformed state)
TW.setState(1)  

# compute deformation for each body
TW.computeDeformations()  

# compute stress tensor for each body
stresses = bodyStressTensors() 

# get strain tensor for each body
strains = [Matrix3(*[TW.deformation(b.id,i,j) for i in range(0,3) for j in range(0,3)]) for b in O.bodies]

###################################################
# save data to vtk. file
__builtin__.my_stresses = stresses 
__builtin__.my_strains = strains 
vtk = export.VTKExporter('result')
vtk.exportSpheres(what=[('radius','b.shape.radius'),('displacement','b.state.displ()'),('stress','my_stresses[b.id]'),('strains','my_strains[b.id]')])

with regards, Alexander

-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.