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,
Here is the code again.
Best
V.
******************************************************************
# 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')
print 'End building packing'
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();O.wait()
#*********************************************************************************************************************
# 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
