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.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
