Plotted variables mismatch

Asked by Panos on 2019-12-16

Hello everyone,
I am having trouble plotting the forces, exerted on a facet object (Disk), due to its interaction with synthetic rock (DEM particles). The forces are plotted against distance (dz), with distance being the position of the disk, at every time step during the cutting process.
The thing is that while simulation runs, at some point plotted forces, seem to bounce back, not agreeing with the current disk position (are plotted in the opposite direction).
This only happens, if the facet object rotates and translates at the same time. By setting its angular velocity to zero, it seems to work fine.
What could I do to fix this problem?
Could you take a look in my code? Thanks a lot.

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'

#### Material microproperties
readParamsFromTable(noTableOk=True,
 Sphere_Radius = 2e-3,
 intR=1.1,
 DENSITY=2600,
 YOUNG=30e9,
 FRICTION_ANGLE=10,
 POISSON=2.5,
 TENS=10e6,
 COH=100e6,
 Z=165.5e-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.02735,-0.1097,0.025),(0.02735,-0.055,Z)),
 radius = Sphere_Radius,
 spheresInCell = 1000,
 material=Sample,
 rRelFuzz=0.34,
 memoizeDb = '/tmp/Test.db',
 returnSpherePack = True,
)
specimen.toSimulation()

#### Disk Creation
Disk=geom.facetCylinder(center=(0,0,0),radius=0.06,height=0.01,segmentsNumber=50,wallMask=7,closeGap=True,orientation=Quaternion((0,1,0), pi/2))
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(1),initRun=True,command='recorder()',label='data'),
CombinedKinematicEngine(ids=ids,comb=[
 TranslationEngine(translationAxis=[0,0,1],velocity=+500,label="trans"),
 RotationEngine(angularVelocity=60000,rotationAxis=[1,0,0],rotateAroundZero=True,zeroPoint=(0.0,0.0,0.0),label="rot")],label='Eng'),
PyRunner(iterPeriod=1,command="rot.zeroPoint+=(0,0,trans.velocity*O.dt)"),
 PyRunner(iterPeriod=int(1), command='Stop()',label='data')
]
O.dt=.5e-4*PWaveTimeStep()

#################
def recorder():
   global Fz, Fx, Fy
   global dz
   Fz = 0.0
   dz = 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,
        'FN':Fy,
        'FR':Fz,
        'dz':(Disk[2].state.pos[2])})
   plot.saveDataTxt(OUT)

def Stop():
    if Disk[2].state.pos[2] >= 190e-3:
        O.pause()

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

O.step();
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

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-12-17
Last query:
2019-12-17
Last reply:
2019-12-16
Jan Stránský (honzik) said : #1

Hello,

> What could I do to fix this problem?

is it really a problem? From the description it seems that it is just the result.. if it is related to rotation of the disk, probably it is related to friction (maybe something like on one side it is present and on the other side not).
I suggest to investigate more the situation when the "problem" occur (investigate inidividual facets, if there is some excessive penetration depth etc.)
cheers
Jan

Panos (panoschristo) said : #2

Hello Dr. Jan Stránský
Thank you for your response,

You commented if that’s a problem: I believe it is. During experimental simulations that does not occur. It is not the result i was hoping for. That’s why I believe there is something wrong. Maybe in the way I record forces?
Further more:
>if it is related to rotation of the disk, probably it is related to friction
(Friction between Facets and the discrete elements?)
But would not friction also exist in the case, which the Disk only translates?
>if there is some excessive penetration depth
I need this to be the penetration depth. It is a variable, that I have to change and compare the results. But the "problem" also occurs in less excessive penetration depths.

Thanks for your comments

Jan Stránský (honzik) said : #3

Sorry,
I have run your simulation (have not done it before) and now I see what you have meant. The problem is not forces (as I understood from the OP), but recorded "dz".
Disk[2].state.pos is "center" of a facet, which is rotating around disk center and naturally "bouncing back" in one half of the rotation.

do
'dz': trans.velocity*O.dt # disk center, same as setting zeroPoint
instead

Also consider to do this change in the stop condition

cheers
Jan

Panos (panoschristo) said : #4

Dr. Stránský,

I did as you suggested, and replaced:
 'dz':(Disk[2].state.pos[2]) with 'dz': trans.velocity*O.dt
As well as in the stopping condition.
Am afraid it does not seem to work. The plot starts from a non zero point(=trans.velocity*O.dt) and then forces plot only in the y – axis of the graph, rather than, along the direction of the cut.

Best Jan Stránský (honzik) said : #5

sorry, I just copy-pasted it from the PyRunner, should have been
'dz': trans.velocity*O.time
cheers
Jan

Panos (panoschristo) said : #6

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