New question #700164 on Yade:
https://answers.launchpad.net/yade/+question/700164
Hi There,
I'm trying to model a soil-wheel interaction using MPI. Originally, my DEM
model consisted of a larger of highly irregular-shaped clump particles but
since YADE MPI has not been designed for clumps, I decided to run the same
problem with spherical particles from those clumps. I have set up the model and
it started working but no collision occurs between my external moving object
which is a facet wheel and the underlying spheres. This is a two-phase problem
where the wheel first penetrates into the soil and then starts rotating in the
longitudinal direction where the position of the zero point is continually
updated.
Here is my code missing the DEM geometry (which is imported).
As most MPI exampmples in YADE deal with gravity deposition only, I'm wondering
if it is possible to model soil-tool interactions like the one I have here.
1- Am I missing something regarding collision establishment?
2- Do I need any adjustments for the rank of my facet wheel to enable collision
and getting forces on the wheel?
Thanks so much
from yade import pack, export
from yade import pack, plot
from yade import ymport
import sys
import os
import os.path
import numpy
import math
import os
from yade import mpy as mp
from yade import timing
NSTEPS=100 #turn it >0 to see time iterations, else only initilization
numThreads = 9 # number of threads to be spawned, (in interactive mode).
############################################
### INPUT PARAMETERS ###
############################################
## Wheel PROPERTIES
rmean=5.56539971721455E-05
widthWheel=0.00125
Wheel_R=0.0025
Wheel_AdditionalNormalForce=-19.6
##WHELL MOVEMENT
slidingratio=0.3
velocity=0.1 #m/s THAT SHOULD BE 0.01 m/s
angularVelocity=velocity/((1-slidingratio)*Wheel_R) #RAD/S
gravityAcc=-9.81
deposFricDegree = 28.5 # INITIAL CONTACT FRICTION FOR SAMPLE PREPARATION
normalDamp=0.7 # NUMERICAL DAMPING
shearDamp=0.7
youngSoil=0.7e8# CONTACT STIFFNESS FOR SOIL
youngContainer=210e9 # CONTACT STIFFNESS FOR CONTAINER
poissonSoil=0.3 # POISSION'S RATIO FOR SOIL
poissionContainer=0.25 # POISSION'S RATIO FOR CONTAINER
densSoil=2650 # DENSITY FOR SOIL
densContainer=7850 # DENSITY FOR CONTAINER
numDamp=0.4
binSize=0.4*Wheel_R
activeDistance=0.01
iniDistanceWheelfromBoundary=0.004
differenceWCenterActiveEnd=activeDistance - iniDistanceWheelfromBoundary
initialPeneterationofWheel=0.001
iniWheelDisfromMaxY=0.0001
############################################
### DDEFINING VARIABLES AND MATERIALS ###
############################################
SoilId=FrictMat(young=youngSoil,poisson=poissonSoil,frictionAngle=radians(deposFricDegree),density=densSoil)
O.materials.append(SoilId)
ContainerId=FrictMat(young=youngContainer,poisson=poissionContainer,frictionAngle=radians(0),density=densContainer)
O.materials.append(ContainerId)
###################################
##### CREATING GEOMETRIES #####
###################################
## IMPORTING THE ALREADY COMPACTED BIN
O.bodies.append(ymport.text("/home/ngoudarzi/Desktop/SWI MPI/Half
Width/5by6by20_n057_Relaxed_FreeSurface_Spheres.txt",shift=Vector3(0,0,0),material=SoilId))
minX=min([b.state.pos[0]-b.shape.radius for b in O.bodies if
isinstance(b.shape,Sphere)])
maxX=max([b.state.pos[0]+b.shape.radius for b in O.bodies if
isinstance(b.shape,Sphere)])
minY=min([b.state.pos[1]-b.shape.radius for b in O.bodies if
isinstance(b.shape,Sphere)])
maxY=max([b.state.pos[1]+b.shape.radius for b in O.bodies if
isinstance(b.shape,Sphere)])
minZ=min([b.state.pos[2]-b.shape.radius for b in O.bodies if
isinstance(b.shape,Sphere)])
maxZ=max([b.state.pos[2]+b.shape.radius for b in O.bodies if
isinstance(b.shape,Sphere)])
print
("minX:",minX,"maxX:",maxX,"minY:",minY,"maxY:",maxY,"minZ:",minZ,"maxZ:",maxZ)
mn,mx=Vector3(minX,minY,minZ),Vector3(maxX,0.01,maxZ)
## CREATE WALLS AROUND THE PACKING
walls=aabbWalls([mn,mx],thickness=0,material=ContainerId)
wallIds=O.bodies.append(walls)
for b in O.bodies:
if b.state.pos[1]>0.009:
O.bodies.erase(b.id)
widthWheel=maxX-minX
XcenterofWheel = widthWheel/2+minX
YcenterofWheel = Wheel_R+maxY+iniWheelDisfromMaxY
ZcenterofWheel = iniDistanceWheelfromBoundary+minZ
Wheel1 =
geom.facetCylinder((XcenterofWheel,YcenterofWheel,ZcenterofWheel),radius=Wheel_R,height=widthWheel,orientation=Quaternion((0,1,0),-pi/2),wallMask=7,segmentsNumber=1000,material=ContainerId,wire=True,color=Vector3(255,255,255))
O.bodies.append(Wheel1)
idsr = [w.id for w in Wheel1]
facets = [b for b in O.bodies if isinstance(b.shape,Facet)] # list of facets in
simulation
for b in facets:
O.forces.setPermF(b.id,(0,Wheel_AdditionalNormalForce/800,0))
############################
### DEFiniNG ENGINES ###
############################
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Facet_Aabb()],
label="collider"),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_MindlinPhys(betan=normalDamp,betas=shearDamp,label='ContactModel')],
[Law2_ScGeom_MindlinPhys_Mindlin(label='Mindlin')]
),
## WE WILL USE THE GLOBAL STIFFNESS OF EACH BODY TO DETERMINE AN OPTIMAL
TIMESTEP (SEE
HTTPS://YADE-DEM.ORG/W/IMAGES/1/1B/CHAREYRE&VILLARD2005_LICENSED.PDF)
#GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.25,parallelMode=True,
label = 'timeStepper'),
CombinedKinematicEngine(ids=idsr,label='combEngine') +
TranslationEngine(translationAxis=(0,-1,0),velocity=velocity) +\
RotationEngine(rotationAxis=(1,0,0), angularVelocity=0, rotateAroundZero=True,
zeroPoint=(XcenterofWheel,YcenterofWheel,ZcenterofWheel)),
NewtonIntegrator(damping=numDamp,gravity=(0,gravityAcc,0)),
VTKRecorder(fileName='/home/ngoudarzi/Desktop/SWI MPI/Half
Width/VTKWheel/Output', recorders=['all'], parallelMode=True,iterPeriod=500),
#use .pvtu to open spheres, .pvtp for ints, and .vtu for boxes.
PyRunner(iterPeriod=1, command="WheelPenetration()" ,label='checker'),
PyRunner(iterPeriod=100,command='history()',label='recorder'),
PyRunner(iterPeriod=1000,command='eraseOffs()'),
]
O.dt = 1e-7
collider.verletDist = 0.05*rmean
# get TranslationEngine and RotationEngine from CombinedKinematicEngine
transEngine, rotEngine = combEngine.comb[0], combEngine.comb[1]
def WheelPenetration():
transEngine.velocity = 7*velocity
rotEngine.zeroPoint += Vector3(0,-1,0)*transEngine.velocity*O.dt
timing.reset()
t=time.time()
print("rotEngine.zeroPoint[1]-Wheel_R:
",(rotEngine.zeroPoint[1]-Wheel_R),"Iteration",O.iter,"Ellapsed Time",O.time)
if (rotEngine.zeroPoint[1]-Wheel_R)<=(maxY-initialPeneterationofWheel):
print("PENETRATION PERFORMED, ROLLING IS BEING STARTED")
checker.command='WheelRolling()'
def WheelRolling():
transEngine.translationAxis=(0,0,1)
transEngine.velocity = velocity
rotEngine.angularVelocity = angularVelocity
rotEngine.zeroPoint += Vector3(0,0,1)*velocity*O.dt
print("rotEngine.zeroPoint[2]:
",rotEngine.zeroPoint[2],"Iteration",O.iter,"Ellapsed Time",O.time)
if rotEngine.zeroPoint[2]>=maxZ-1.05*Wheel_R:
print("ROLLING PERFORMED")
O.pause()
def history():
global Fx,Fy,Fz,Dx,Dy,Dz
Fx=0
Fy=0
Fz=0
for b in facets:
Fx+=O.forces.f(b.id,sync=True)[0]
Fy+=O.forces.f(b.id,sync=True)[1]
Fz+=O.forces.f(b.id,sync=True)[2]
Dx=rotEngine.zeroPoint[0]
Dy=rotEngine.zeroPoint[1]
Dz=rotEngine.zeroPoint[2]
yade.plot.addData({'i':O.iter,'Fx':Fx,'Fy':Fy,'Fz':Fz,'Dx':Dx,'Dy':Dy,'Dz':Dz,})
## In that case we can still save the data to a text file at the the end of
the simulation, with:
plot.saveDataTxt('/home/ngoudarzi/Desktop/SWI MPI/Half
Width/Wheel_OutputData_5by6by20_n057_Relaxed_FreeSurface.txt')
## declare what is to plot. "None" is for separating y and y2 axis
plot.plots={'i':('Fx','Fy','Fz','Dx','Dy','Dz')}
plot.plot()
######### RUN ##########
# customize mpy
mp.ERASE_REMOTE_MASTER = True #keep remote bodies in master?
mp.DOMAIN_DECOMPOSITION= True #automatic splitting/domain decomposition
#mp.mpirun(NSTEPS) #passive mode run
mp.MERGE_W_INTERACTIONS = False
mp.mpirun(NSTEPS,numThreads,withMerge=True) # interactive run, numThreads is
the number of workers to be initialized, see below for withMerge explanation.
mp.mergeScene() #merge scene after run.
if mp.rank == 0: O.save('/home/ngoudarzi/Desktop/SWI MPI/Half
Width/mergedScene.yade')
#demonstrate getting stuff from workers, here we get kinetic energy from worker
subdomains,notice that the master (mp.rank = 0), uses the sendCommand to tell
workers to compute kineticEnergy.
if mp.rank==0:
print("kinetic energy from workers:
"+str(mp.sendCommand([1,2],"kineticEnergy()",True)))
--
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 : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
