Question #688203 on Yade changed:
https://answers.launchpad.net/yade/+question/688203
Status: Needs information => Open
Rioual gave more information on the question:
...Here is the code:
#Start
import time
from yade import qt
from yade import ymport
#Import packing in the cylinder from prepare-packing.py
O.load('init_Final_state_packing.yade')
###########################################
## 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 Warren-Spring
global fctIdsWS
fctIdsWS=O.bodies.append(ymport.stl('WS-FR-1.stl',color=(1,0,0),material=facetMat))
#fctIdscylinder=O.bodies.append(ymport.stl('PotN.stl',color=(1,0,0),material=facetMat))
## Timestep
O.dt=.2*tc
## Engines
O.engines=[
## Resets forces and momenta the act on bodies
ForceResetter(),
## Associates bounding volume to each body.
#BoundDispatcher([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
## Using bounding boxes find possible body collisions.
#InsertionSortCollider(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_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()],
[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,-1,0],velocity=10)
O.engines=O.engines+[TransEngload2]
while h_WS > hmaxSpheres:
print 'test01'
TransEngload2.dead = False
print 'test02'
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
minY = 1e99
maxY = -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
minY = min(minY,min(v[1] for v in vs))
maxY = max(maxY,max(v[1] for v in vs)) ###
# print 'maxY=',maxY,'minY=',minY
h_WS = maxY
###########################################
#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[1] # 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
#*****************************************************************************************************************************
...So this is the complete script, maybe the problem is in relation with
the problem of dealing with an aditional engine in a while loop (?) . I
am not sure of myself on this....
Thanks for your feedback,
best
Vincent,
--
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
