stl geometry problem

Asked by Panos on 2019-11-04

Hello everyone,
I am trying to simulate the rock cutting process with a disk cutter.
During the cut the disk translates and rotates. I have done that by using
Kinematic Engine.
The problem is with the geometry.
When disk rotates, deforms. The facets forming the disk seem to detach.
I have created the geometry (Solid), with 3 different ways (Blender, 3d builder and Autocad) without solving the problem.
Any suggestions on how can I fix it;
I am wondering it might affect forces exerted on the disk.
I am using Ubuntu 18.04 and Yade 2018.02b.
Thank you
I am sending the following: A file of the disk geometry in .stl format, and my code.

https://1drv.ms/u/s!AiKcYaY3wMDqgV9IE7I2RFUWVTxp

from __future__ import print_function
from __future__ import division
from yade import pack, plot, geom, pack,utils,ymport,export,bodiesHandling
from yade.pack import *
from yade import _packPredicates
import math
from yade import *

################# SIMULATIONS DEFINED HERE
#### Simulation
OUT='LinearCutTest_JCFPM'

#### Simulation Control
saveVTK=10000 # saving output files for paraview

#### Material microproperties
readParamsFromTable(noTableOk=True, # unknownOk=True,
 Sphere_Radius = 1.5e-3,
 intR=1.25, # allows near neighbour interaction (GammaInt)
 DENSITY=2700, # Density of the material (kg/m^3),could be adapted to match material density: dens_DEM=dens_rock*(V_rock/V_particles)=dens_rock*1/(1-poro_DEM) -> packing porosity as to be computed?
 YOUNG=68e9, # Eeq(GPa),elastic modulus
 FRICTION_ANGLE=12, #(Degrees),Contact friction angle
 POISSON=1,
 TENS=22e6, #(Tensile Strenght T (MPa)),Defines the maximum admissible normal force in traction in the matrix (FnMax = tensileStrength * crossSection). (Pa)
 COH=220e6, #(Local Cohesion C (MPa)),Defines the maximum admissible tangential force in shear, for Fn=0, in the matrix (FsMax = cohesion * crossSection). (Pa)

### Cutting_Specimen_dimensions
 X = 54.85e-3,
 Y=54.85e-3,
 Z=180e-3,

)

from yade.params.table import *

#### material definition
Sample = O.materials.append(JCFpmMat(
 young=YOUNG,
 poisson=POISSON,
 frictionAngle=radians(FRICTION_ANGLE),
 cohesion=COH,
 tensileStrength=TENS,
 density = DENSITY,
 label='spheres'
))

#### create the specimen
specimen=pack.randomDensePack(
 pack.inAlignedBox((0,0,0),(X,Y,Z)),
 radius = Sphere_Radius,
 spheresInCell = 1000,
 material=Sample,
 rRelFuzz=0.34,
 memoizeDb = '/tmp/Linear-Cutting_Test.db',
 returnSpherePack = True,
)
specimen.toSimulation()

#### Disk Creation
Disk=ymport.stl('./p(5mm).stl')
O.bodies.append(Disk)
ids = [b.id for b in Disk]

################# ENGINES DEFINED HERE

O.engines=[
 ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb'),Bo1_Facet_Aabb()],verletDist=.05*Sphere_Radius),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom'),Ig2_Facet_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.5,defaultDt=utils.PWaveTimeStep()),
 NewtonIntegrator(damping=0.5,label='newton'),
 PyRunner(iterPeriod=int(100),initRun=True,command='recorder()',label='data'),
     VTKRecorder(iterPeriod=int(saveVTK),initRun=True,fileName=OUT+'-',recorders=['spheres','jcfpm','cracks','facets'],Key=OUT,label='vtk'),
 PyRunner(iterPeriod=int(100), command='Stop()',label='data'),
 CombinedKinematicEngine(ids=ids,comb=[
 TranslationEngine(translationAxis=[0,0,1],velocity=-580),
 RotationEngine(angularVelocity=6700,rotationAxis=[1,0,0],rotateAroundZero=False)],label='Eng')
]

################# RECORDER DEFINED HERE
def recorder():
   global Fz, Fx, Fy
   global dz
   Fz = 0.0
   Fz = abs(sum(O.forces.f(facet.id)[2] for facet in Disk))
   Fx = abs(sum(O.forces.f(facet.id)[0] for facet in Disk))
   Fy = abs(sum(O.forces.f(facet.id)[1] for facet in Disk))
   yade.plot.addData({'i':O.iter,
        'FS':Fx,
        'FN':Fy,
        'FR':Fz,
        'dz': abs(Disk[2].state.pos[2]),
        'tc':interactionLaw.nbTensCracks,
        'sc':interactionLaw.nbShearCracks,
        'te':interactionLaw.totalTensCracksE,
        'se':interactionLaw.totalShearCracksE,
        'unbF':utils.unbalancedForce()})
   plot.saveDataTxt(OUT)

def Stop():
    if Disk[2].state.pos[2] == 0:
        O.pause()

# if you want to plot during simulation
plot.plots={'dz':('FR','FN')}
plot.plot()

#### manage interaction detection factor during the first timestep and then set default interaction range ((cf. A DEM model for soft and hard rock, Scholtes & Donze, JMPS 2013))
O.step();
### initializes the interaction detection factor
SSgeom.interactionDetectionFactor=1.
Saabb.aabbEnlargeFactor=1.

# run simulation
plot.plot()
#O.run()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2019-11-06
Last query:
2019-11-06
Last reply:
2019-11-04
Best Jan Stránský (honzik) said : #1

Hello,

please read [1] and (next time) try to:
- make the question self-contained (see section "Please, no external links!")
- try to make the code MWE

M in MWE stands for minimal. In this case of moving a disk, specimen, material, recorder... are irrelevant.
In this light, you can create a few facets manually, still preserving the core of the problem and breaking the need of external files.
Creating a MWE and localizing the problem, often you can solve the problem yourself.

Solution:
- RotationEngine(...,rotateAroundZero=True) # [2]
- specify zeroPoint of RotationEngine [3]
- (probably) you will need to update zeroPoint each time step according to the translation of the cutter

a MWE:
###
O.bodies.append((
 facet(((0,0,0),(0,-1,0),(0,0,-1))),
 facet(((0,0,0),(0,-1,0),(0,0,+1))),
 facet(((0,0,0),(0,+1,0),(0,0,-1))),
 facet(((0,0,0),(0,+1,0),(0,0,+1))),
))
ids = (0,1,2,3)
O.engines = [
   ForceResetter(),
   CombinedKinematicEngine(ids=ids,comb=[
      TranslationEngine(translationAxis=[0,0,1],velocity=-580,label="trans"),
      RotationEngine(angularVelocity=6700,rotationAxis=[1,0,0],rotateAroundZero=True,zeroPoint=(0,0,0),label="rot")
   ],label='Eng'),
   NewtonIntegrator(),
   PyRunner(iterPeriod=1,command="rot.zeroPoint+=(0,0,trans.velocity*O.dt)")
]
O.dt = 1e-8
###

cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask
[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.RotationEngine.rotateAroundZero
[3] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.RotationEngine.zeroPoint

Panos (panoschristo) said : #2

Thank you Jan Stránský for your response and sorry for replying late to it, i was trying to
make it work.
I did what you point me and it seems that the problem is fixed.
I only changed O.dt=PWaveTimeStep()

I am having a different problem a mismatch between Disk’s position and
the data plotted on a graph. But I wil open a new thread if I don’t figure it.
Thank you.

Panos (panoschristo) said : #3

Thanks Jan Stránský, that solved my question.