yade-users team mailing list archive

[Nima Goudarzi <question700165@xxxxxxxxxxxxxxxxxxxxx>]

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

Hi,

I'm trying to model a soil-wheel interaction using FEM-DEM. The upper part of my model is DEM clumps which are then are coupled (volume coupling) from the bottom to a FEM mesh. 

Here I used a spherical transition zone and avoided clumps coupling but am also wondering if clump coupling is also supported.

The motions of the facet wheel consist of:
1- Penetration in Y direction up to a certain depth (in Y direction)
2- Rotation in the longitudinal direction using a sliding ratio (in Z direction)

If I model the same problem with pure DEM particles, I define a CombinedKinematicsEngine and play with its components in different functions declared after the O.engines setup. 
Unfortunately, when I use FEM-DEM coupling, the functions in YADE model are not recognized and this excessively limits my options for modeling different motions.
Noteworthy that  I need to use PyRunner to update the zero point of the wheel to disable/enable my engine's components.
The basic question is this:

Do functions literally work in FEM-DEM coupling? 
If yes do I need to adjust something scripts other than the YADE mode?

Here is my scrip with non-functional functions
The geometry is missing but I only need to know if something is wrong with my functions declaration. 
######################################################################
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

######################################################
###   DEFINING VARIABLES AND MATERIAL PARAMETERS   ###
######################################################
deposFricDegree = 28.5 # INITIAL CONTACT FRICTION FOR SAMPLE PREPARATION
normalDamp=0.4 # NUMERICAL DAMPING
shearDamp=0.4
youngSoil=0.7e8# CONTACT STIFFNESS FOR SOIL
youngEO=210e9 # CONTACT STIFFNESS FOR EXTERNAL OBJECTS
poissonSoil=0.3 # POISSION'S RATIO FOR SOIL
poissionEO=0.25 # POISSION'S RATIO FOR EXTERNAL OBJECTS
densSoil=2650 # DENSITY FOR SOIL
densEO=500 # DENSITY FOR EXTERNAL OBJECTS
numDamp=0.4
IniDistanceBladefromBoundary = 0.001
HeightofBlade=0.003
WidthofBlade=0.002
ThicknessofBlade=0.0002
InitialPeneterationofBlade = 0.0005
InitialDistanceofBladefromTopSoil = 0
HorizentalvelocityofBlade = 10
PenetrationvelocityofBlade = 10

SoilId=FrictMat(young=youngSoil,poisson=poissonSoil,frictionAngle=math.radians(deposFricDegree),density=densSoil)
O.materials.append(SoilId)

EOId=FrictMat(young=youngEO,poisson=poissionEO,frictionAngle=math.radians(0),density=densEO)
O.materials.append(EOId)
#O.bodies.append([utils.sphere(c,r) for c,r in sp])
#for b in O.bodies:
#  b.state.blockedDOFs = 'zXY'
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-radius),ZcenterofContainer),wallMask=51,material=EOId)
O.bodies.append(Container)
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)
## CREATE WALLS AROUND THE PACKING
# radius=8e-5
# mn,mx=Vector3(minX,minY,minZ),Vector3(maxX, maxY, maxZ)
# walls=yade.utils.aabbWalls([mn,mx],thickness=0,material=EOId)
# wallIds=O.bodies.append(walls)
# O.bodies.erase(517222)
# O.bodies.erase(517223)
XcenterofBlade = 0.0025
YcenterofBlade = maxY+(HeightofBlade/2)
ZcenterofBlade = minZ+IniDistanceBladefromBoundary
Blade=yade.geom.facetBox((XcenterofBlade,YcenterofBlade,ZcenterofBlade), (WidthofBlade/2,HeightofBlade/2,ThicknessofBlade/2),material=EOId)
O.bodies.append(Blade)

idsr = [w.id for w in Blade]
facets = [b for b in O.bodies if isinstance(b.shape,Facet)] # list of facets in simulation
############################
###   DEFINING ENGINES   ###
############################
gravity=(0,-9.81,0)
O.engines = [
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Facet_Aabb()]),
  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')]),
  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)),
  PyRunner(iterPeriod=1, command="WheelPenetration()" ,label='checker'),
  NewtonIntegrator(gravity=gravity, damping=.1),
]
O.dt = 1e-7
O.step()
transEngine, rotEngine = combEngine.comb[0], combEngine.comb[1]

def WheelPenetration():
  transEngine.velocity = 7*velocity
  rotEngine.zeroPoint += Vector3(0,-1,0)*transEngine.velocity*O.dt
  print("rotEngine.zeroPoint[1]-((HeightofBlade/2)):", (rotEngine.zeroPoint[1]-((HeightofBlade/2))),"Iteration()",O.iter)
  if (rotEngine.zeroPoint[1]-Wheel_R)<=(maxY-initialPeneterationofWheel):
    print("Wheel penetrated and is starting to rotate")
    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)
  if rotEngine.zeroPoint[2]>=maxZ-1.05*Wheel_R:
    O.pause()

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.