yade-users team mailing list archive

[Nima Goudarzi <question700164@xxxxxxxxxxxxxxxxxxxxx>]

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.