Question #688203 on Yade changed:
https://answers.launchpad.net/yade/+question/688203
Status: Needs information => Open
Rioual gave more information on the question:
Hello Jan,
Sorry about all that. I modified my code to make it independant with a faceted
cylinder for penetration and Oedometric-test.py used to make the packing. Here
it is below.
The simulation seems to be blocked for more than several minutes. Thee number
of particles here is very low
(example of packing taken from the tutorial examples oedometric-test.py).
All the best,
Vincent
The code:
*****************************************************************************************************************************
# gravity deposition, continuing with oedometric test after stabilization
# shows also how to run parametric studies with yade-batch
# The components of the batch are:
# 1. table with parameters, one set of parameters per line (ccc.table)
# 2. readParamsFromTable which reads respective line from the parameter file
# 3. the simulation muse be run using yade-batch, not yade
#
# $ yade-batch --job-threads=1 03-oedometric-test.table 03-oedometric-test.py
#
# load parameters from file if run in batch
# default values are used if not run from batch
readParamsFromTable(rMean=.05,rRelFuzz=.3,maxLoad=1e6,minLoad=1e4)
# make rMean, rRelFuzz, maxLoad accessible directly as variables later
from yade.params.table import *
# create box with free top, and ceate loose packing inside the box
from yade import pack, plot
O.bodies.append(geom.facetBox((.5,.5,.5),(.5,.5,.5),wallMask=31))
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,1),rMean=rMean,rRelFuzz=rRelFuzz)
sp.toSimulation()
O.engines=[
ForceResetter(),
# sphere, facet, wall
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
InteractionLoop(
# the loading plate is a wall, we need to handle sphere+sphere,
sphere+facet, sphere+wall
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
NewtonIntegrator(gravity=(0,0,-9.81),damping=0.5),
# the label creates an automatic variable referring to this engine
# we use it below to change its attributes from the functions called
PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),
]
O.dt=.5*PWaveTimeStep()
# the following checkUnbalanced, unloadPlate and stopUnloading functions are
all called by the 'checker'
# (the last engine) one after another; this sequence defines progression of
different stages of the
# simulation, as each of the functions, when the condition is satisfied,
updates 'checker' to call
# the next function when it is run from within the simulation next time
# check whether the gravity deposition has already finished
# if so, add wall on the top of the packing and start the oedometric test
def checkUnbalanced():
# at the very start, unbalanced force can be low as there is only few
contacts, but it does not mean the packing is stable
if O.iter<5000: return
# the rest will be run only if unbalanced is < .1 (stabilized packing)
if unbalancedForce()>.1: return
# add plate at the position on the top of the packing
# the maximum finds the z-coordinate of the top of the topmost particle
O.bodies.append(wall(max([b.state.pos[2]+b.shape.radius for b in
O.bodies if isinstance(b.shape,Sphere)]),axis=2,sense=-1))
global plate # without this line, the plate variable would only
exist inside this function
plate=O.bodies[-1] # the last particles is the plate
# Wall objects are "fixed" by default, i.e. not subject to forces
# prescribing a velocity will therefore make it move at constant
velocity (downwards)
plate.state.vel=(0,0,-.1)
# start plotting the data now, it was not interesting before
O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=200)]
# next time, do not call this function anymore, but the next one
(unloadPlate) instead
checker.command='unloadPlate()'
def unloadPlate():
# if the force on plate exceeds maximum load, start unloading
if abs(O.forces.f(plate.id)[2])>maxLoad:
plate.state.vel*=-1
# next time, do not call this function anymore, but the next
one (stopUnloading) instead
checker.command='stopUnloading()'
def stopUnloading():
if abs(O.forces.f(plate.id)[2])<minLoad:
# O.tags can be used to retrieve unique identifiers of the
simulation
# if running in batch, subsequent simulation would overwrite
each other's output files otherwise
# d (or description) is simulation description (composed of
parameter values)
# while the id is composed of time and process number
# plot.saveDataTxt(O.tags['d.id']+'.txt')
O.pause()
def addPlotData():
if not isinstance(O.bodies[-1].shape,Wall):
plot.addData(); return
Fz=O.forces.f(plate.id)[2]
#
plot.addData(Fz=Fz,w=plate.state.pos[2]-plate.state.refPos[2],unbalanced=unbalancedForce(),i=O.iter)
# besides unbalanced force evolution, also plot the displacement-force diagram
#plot.plots={'i':('unbalanced',),'w':('Fz',)}
#plot.plot()
O.run()
#*********************************************************************************************************************
# Let the object penetrate the packing
#Start penetration
import time
from yade import qt
from yade import ymport
from yade.gridpfacet import *
###########################################
## PhysicalParameters
Density=2400
frictionAngle=radians(35)
tc = 0.001
en = 0.05
et = 0.05
## Import wall's geometry
facetMat=O.materials.append(ViscElMat(frictionAngle=frictionAngle,tc=tc, en=en,
et=et)) # **params sets kn, cn, ks, cs
sphereMat=O.materials.append(ViscElMat(density=Density,frictionAngle=frictionAngle,tc=tc,en=en,et=et))
from yade import ymport
#******************************************
global TransEngload2
#******************************************
# Importation de la roue
global fctIdsWS
fctIdsWS=O.bodies.append( facetCylinder(center=(0.5,0.5,1.2),
radius=0.3, height=0.2, orientation=Quaternion((1, 0, 0), 0),
segmentsNumber=10, wallMask=7, angleRange=None, closeGap=False,
radiusTopInner=-1, radiusBottomInner=-1, material=facetMat) )
## Timestep
O.dt=.2*tc
## Engines
O.engines=[
## Resets forces and momenta the act on bodies
ForceResetter(),
## Using bounding boxes find possible body collisions.
#InsertionSortCollider(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_GridConnection_Aabb()]),
## Interactions
InteractionLoop(
# the loading plate is a facet, we need to handle
sphere+sphere, sphere+facet
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom(),Ig2_Sphere_GridConnection_ScGridCoGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
## Apply gravity
GravityEngine(gravity=[0,0,-9.81]),
## Cundall damping must been disabled!
NewtonIntegrator(damping=0),
## Apply kinematics to wheel
PyRunner(command='kinematics_WS()',realPeriod=1,label='kine'),
]
from yade import qt
qt.View()
def kinematics_WS():
h_WS = calc_h()[0]
hmaxSpheres = calc_h()[1]
print 'test0','h_WS=',h_WS,'hmaxSpheres=',hmaxSpheres
TransEngload2 =
TranslationEngine(ids=fctIdsWS,translationAxis=[0,0,-1],velocity=10)
O.engines=O.engines+[TransEngload2]
while h_WS > hmaxSpheres:
print 'test01'
TransEngload2.dead = False
print 'test02'
######################### Here the simulation seems to be
blocked#################################################################
O.run(1,True)
# O.engines=O.engines+[PyRunner(command='calc_h()')]
print 'test03'
h_WS = calc_h()[0]
hmaxSpheres = calc_h()[1]
print 'test1','h_WS=',h_WS,'hmaxSpheres=',hmaxSpheres
TransEngload2.dead = True
else:
amTOT = sum((O.bodies[facetid].state.angMom)[1] for
facetid in fctIdsWS)
while (amTOT< 1e-10):
O.engines=O.engines+[PyRunner(command='addTorque()',iterPeriod=1)]
amTOT = sum((O.bodies[facetid].state.angMom)[1] for
facetid in fctIdsWS)
else:
# Stop simulation et measurement of the shear torque and cohesion
C_WS = 3* (imposed_T)/(2*Pi*(R0^3-R1^3))
print 'C_WS=',C_WS
print '********End of the calculation for
Cohesion***********' #
O.pause()
#*********************************************************************************************************************************
#This function calculate the height of the boundary h_WS and the maximum height
of the packing hmax_Spheres
def calc_h():
###########################################
#Calculate h_WS
minZ = 1e99
maxZ = -1e99
for facet in fctIdsWS:
# print 'facet=', facet
facet= O.bodies[facet]
vs = [facet.state.pos + facet.state.ori * v for v in
facet.shape.vertices] # vertices in global coord system
# print 'vs=',vs
minZ = min(minZ,min(v[2] for v in vs))
maxZ = max(maxZ,max(v[2] for v in vs)) ###
# print 'maxZ=',maxZ,'minZ=',minZ
h_WS = maxZ
###########################################
#Calculate hmax_Spheres
idHMax=0 # on definit une variable pour
identifier quel corps a la hauteur la plus haute
hMax=0.0 # initialisation de la hauteur
a zero
for i in O.bodies: # on parcours tout les corps du
systeme
h=i.state.pos[2] # on extrait la position selon
l'axe y
if (type(i.shape).__name__ == 'Sphere' and h>hMax): #si le
solide est une sphere et sa position est plus haut que hmax
hMax=h # on enregistre sa hauteur
idHMax=i.id # on garde en memoire de quel
corps il s'agit
# if (idHMax == 0): return rMax # valeur de retour par default,
rMAX
# else: return (hMax+O.bodies[idHMax].shape.radius)
hmaxSpheres = hMax+O.bodies[idHMax].shape.radius
return h_WS,hmaxSpheres
#*****************************************************************************************************************************
#*****************************************************************************************************************************
# Add incremental torque on the WS boundary
def addTorque():
increment_T= 0.001
for Idfacets in fctIdsWS:
O.forces.addT(Idfacets,(0,increment_T,0))
imposed_T = O.forces.permT(Idfacets)
return imposed_T
#*****************************************************************************************************************************
--
You received this question notification because your team yade-users 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
