FrictPhys disable in contact

Asked by Pawel on 2019-01-17

Hi everyone,
I have very simple simulation. At the beginning there are two spheres, with material:

O.materials.append(FrictMat(young=young_k, poisson=poisson_k, density=density_k, frictionAngle=atan(0.3),))

and positions:
O.bodies.append([
  utils.sphere(center=(0,0,0),radius=5,fixed=True,material=kula[0])
  ],)
O.bodies.append([
  utils.sphere(center=(0.2*5,0,2*5),radius=promien_k,fixed=False,material=kula[0])
  ],)

In loop:

InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=1.3)
],
  [Ip2_FrictMat_FrictMat_FrictPhys(),
  Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=True)
  ],
  [Law2_ScGeom_FrictPhys_CundallStrack(),
  ]
 ),
   # apply gravity force to particles
   GravityEngine(gravity=(0,0,-9.81)),
   # damping: numerical dissipation of energy
   NewtonIntegrator(damping=0.2)
]

The upper sphere is rotating over the first one - so the simulation is working.

I'm trying to write short program in C with detection of interaction. Two spheres are in interaction but the interaction all of the time is not real - i get 0 in I->isReal()

I can't figure out how to get this to work
FrictPhys* physFP = dynamic_cast<FrictPhys*>(I->phys.get());
tmp = physFP->shearForce;

tmp - is vector

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2019-01-18
Last reply:
2019-01-20
Jan Stránský (honzik) said : #1

Hello Pawel,

I think much more information would be needed to solve your problem.
Please read [1] and provide a complete MWE, ideally both Python and C++.

Is the C (C++ I guess?) program a standalone program using Yade? Is it C++ function to be called from Python? ...?

thanks
Jan

[1] https://yade-dem.org/wiki/Howtoask

Pawel (pawel-p) said : #2

Hi Jan,
thank you for your message

in c++ i have function VppSaveInteractions with code
FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
    if(I)
     {
     const shared_ptr<Body>& b1 = Body::byId(I->getId1(),scene);
     const shared_ptr<Body>& b2 = Body::byId(I->getId2(),scene);

     State* st1 = static_cast<State*>(b1->state.get());
     State* st2 = static_cast<State*>(b2->state.get());

     FrictPhys* physFP = dynamic_cast<FrictPhys*>(I->phys.get());
     cout << "phys "<< physFP;
}
It is only a part of longer code but here i need to gest phys becouse I want to save normal and shear force from connection.
 the pointer physe is 0

here is the python script

# basic simulation showing sphere falling ball gravity,
# bouncing against another sphere representing the support
os.environ['OMP_NUM_THREADS']='8'
import math
import random
# DATA COMPONENTS

promien_k = 0.0002
promien_f = 0.00015/2
young_k = 1e9
poisson_k = 0.25
density_k = 1040

frictAngle = atan(0.6)

O.materials.append(FrictMat(young=young_k, poisson=poisson_k, density=density_k, frictionAngle=atan(0.3), label='kula'))

O.bodies.append([
  utils.sphere(center=(0,0,0),radius=promien_k,fixed=True,material='kula')
  ],)
O.bodies.append([
  utils.sphere(center=(0.2*promien_k,0,2*promien_k),radius=promien_k,fixed=False,material='kula')
  ],)

O.dt=utils.PWaveTimeStep()

#########################################
nazwaPliku = "wyniki/kule"
nazwaPlikuIter = "wyniki/Kule_iter_30z"
############### Parametry zapisu c++ ########
print("DT do zapisu")
step = O.dt
PiterSave = 10 #ile zapisow w pliku
PtimeSr = 3*step
PtimeStart = 0. #start zapisu
PdeltaTime = 0.005 #1000*step #skok miedzy zapisami

TimeStartAll=1000000
DeltaTimeAll=10000
################################
################################
# FUNCTIONAL COMPONENTS

# simulation loop -- see presentation for the explanation
O.engines=[
 VppSaveInteractions(fileName=nazwaPlikuIter, iterSave= PiterSave, timeSr = PtimeSr, timeStart = PtimeStart, deltaTime= PdeltaTime,timeStartAll = TimeStartAll,deltaTimeAll = DeltaTimeAll ),
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1),Bo1_Wall_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=1.3),
  Ig2_Wall_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom6D()],
  [Ip2_FrictMat_FrictMat_FrictPhys(),
  Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=True),
  ],
  [Law2_ScGeom_FrictPhys_CundallStrack(),
  ]
 ),
   # apply gravity force to particles
   GravityEngine(gravity=(0,0,-9.81)),
   # damping: numerical dissipation of energy
   NewtonIntegrator(damping=0.2)
]

# set timestep to a fraction of the critical timestep
# the fraction is very small, so that the simulation is not too fast
# and the motion can be observed

# save the simulation, so that it can be reloaded later, for experimentation
O.saveTmp()
Gl1_Sphere.stripes=True
Gl1_Sphere.quality=1
yade.qt.Controller(), yade.qt.View();

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

Hello Pawel,
the C++ seems ok. I think "if (I)" is not needed, would put "if(I->isReal())" instead, but don't know all the details.

adding this line to O.engines (basically printing O.iter when the interaction is real)
PyRunner(iterPeriod=1,command="if [i for i in O.interactions]: print O.iter"),
I got output:
3506
3507
3508
3509
3510
7839
7840
7841
7842
7843
7844
10835
10836
10837
10838
10839
10840
12857
12858
12859
12860
12861
12862
14274
14275
14276
14277
14278
14279

so the interaction is indeed not real most of the time, probably performing some kind of "micro-bouncing".
You can test the interaction reality and if false, using zero vectors for the forces.

cheers
Jan

Pawel (pawel-p) said : #4

Hi Jan,
Thank you very much for your answer.

Is there any possibility to eliminate that micro-bouncing?

All the best
Pawel

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

I think it is an intrinsic feature of such two sphere interaction..
You can try to place the spheres in just touching position (to prevent the initial bounce), but still I think there will be moments without contact..
cheers
Jan

Can you help with this problem?

Provide an answer of your own, or ask Pawel for more information if necessary.

To post a message you must log in.