yade-users team mailing list archive

[Nima Goudarzi <question700165@xxxxxxxxxxxxxxxxxxxxx>]

Question #700165 on Yade changed:
https://answers.launchpad.net/yade/+question/700165

    Status: Needs information => Open

Nima Goudarzi gave more information on the question:
Thanks so much, Jan

> please be more specific, there are several approaches / projects on this. ​
My problem falls in the soil-tool interaction category. The original DEM geometry is constructed outside this script and is then imported here. That DEM model is mostly composed of clumps of spheres of different shapes and sizes but there exist some standalone spheres as well.

Soil Geometry 
I am interested to implement a volumetric coupling at the bottom of a 2.5 mm DEM model with their original clumps but am not sure if FEM-DEM coupling works for clumps. As a result, I considered a separate transition zone at the bottom of this 2.5 mm which is purely constructed from spherical particles. In a separate script, the clumps are first deposited on this transition and at the end of the deposition, the DEM geometry which now has two parts (spheres for the transition and clumps for the upper 2.5mm layer) are exported. That is why I have two ymport in my script.

This DEM is then coupled from the bottom to a coarse FEM mesh with constrained BCs (all 6 DOFs) at the bottom and sides.
To constrain the lateral movement of the DEM, I have used a facetBox (wallMask=51) which lacks the top and bottom boundaries. Therefore, the model owns two types of boundaries: nodes constraints of the FEM mesh, and the facetBox for the upper layer DEM.

I hope that I have been able to explain the geometry.

Motions
Upon the generation of the model geometry (soil) and application of boundary conditions, I try to model the interaction between a face wheel over the surface of the DEM terrain:
Using a combinedKinematicsEngine:
1- The facet wheel is first penetrated vertically (in the Y direction) up to a certain depth (1mm). In this phase, the Translation component (of the combinedKinematicsEngine) controls the movement of the wheel. The motion of the wheel is continuously monitored using the zero point of the Rotation Engine component which in the penetration phase has zero angular velocity but is NOT dead. I mean that the difference between the updating coordination of the zero point in direction [1] and the wheel radius, determines when to stop the penetration and start the 2 phase (wheel rotation).
2- Upon the arrival of the wheel to the desired penetration depth, the Translation component (of the combinedKinematicsEngine)  is deactivated and the Rotation Engine component imposes the rotational movement of the wheel using an angular velocity which is calculated from the penetration velocity by means of a SLIDING RATIO.

As you can see, I need functions for enabling/disabling different
motions of the wheel.

THE PREVIOUS SCRIPT HAD SOME ISSUES. THIS ONE IS THE UPDATED ONE:
THE MOST ANNOYING ERROR WHILE RUNNING IS THAT THE WheelPenetration() FUNCTION HAS NOT BEEN DEFINED. Therefore the PyRunner doesn't work. As a consequence, the next functions are not executed as well. 

>what does "other mode" mean?

I meant model not mode. I mean if we add some functions in YADE script
(DEM), do we need to change something in the coupling or FEM scripts?

BTW, functions () are not callable and this is my first and foremost issue. 
Note that a pure DEM model with the same configuration works as a charm. 

Thanks so much
######################################################################
######################################################################
import math
from math import sqrt
from libyade import yade
from yade import *
from yade import pack, plot 
import numpy
from pyquaternion import Quaternion
from yade import pack, export, Vector3
from yade import ymport
import sys
import os
import os.path
import IPython
############################################
###           INPUT PARAMETERS           ###
############################################
## Wheel PROPERTIES
widthWheel=0.00125
Wheel_R=0.0025
Wheel_AdditionalNormalForce=-19.6

##WHELL MOVEMENT
slidingratio=0.3
velocity=0.01 #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

############################################
###   DEFINING 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   #####
###################################
## CREATE WALLS AROUND THE PACKING
radius=8.e-5
XcenterofContainer = 0.0025
YcenterofContainer = 0.0054383384
minY=0.00241406
ZcenterofContainer = 10.0e-3
Container=yade.geom.facetBox((XcenterofContainer,YcenterofContainer,ZcenterofContainer), (XcenterofContainer,(YcenterofContainer-minY),ZcenterofContainer),wallMask=51,material=EOId)
O.bodies.append(Container)
## IMPORTING THE ALREADY COMPACTED BIN
ymport.textClumps("/home/ngoudarzi/Desktop/dem-fem-coupling-master/examples/vol5/DEMPClumps.txt",shift=Vector3(0,0,0),material=SoilId)
O.bodies.append(ymport.text("/home/ngoudarzi/Desktop/dem-fem-coupling-master/examples/vol5/DEMSpheres.txt",shift=Vector3(0,0,0),material=SoilId,color=Vector3(0.6,0.6,0.6)))
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)

XcenterofWheel = widthWheel/2+minX
YcenterofWheel = Wheel_R+maxY+iniWheelDisfromMaxY
ZcenterofWheel = iniDistanceWheelfromBoundary+minZ

Wheel1 = yade.geom.facetCylinder((XcenterofWheel,YcenterofWheel,ZcenterofWheel),radius=Wheel_R,height=widthWheel,orientation=Quaternion((0,1,0),-pi/2),wallMask=7,segmentsNumber=200,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),
  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)),
  PyRunner(iterPeriod=1, command="WheelPenetration()" ,label='checker'), #######THIS IS NOT CALLED
  PyRunner(iterPeriod=1000,command='history()',label='recorder'),
  ]
O.dt = 1e-7
O.step()
# 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
  if (rotEngine.zeroPoint[1]-Wheel_R)<=(maxY-initialPeneterationofWheel):
    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
  if rotEngine.zeroPoint[2]>=maxZ-1.05*Wheel_R:
    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 end of the simulation, with:
  plot.saveDataTxt('./Wheel_OutputData_5by6by20_n057_Relaxed_FreeSurface.txt')
def vtkExport(i):
  from yade import export
  name = '/tmp/vol1_yade'
  export.VTKExporter(name,i).exportSpheres()
  export.VTKExporter(name,i).exportFacets()

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.