GlobalStiffnessTimeStepper with Potential Blocks

Asked by Vasileios Angelidakis on 2019-09-16


I want to use the GlobalStiffnessTimeStepper in conjunction with the Potential Blocks code.
Currently, this engine assumes that the stiffnesses "phys->kn" and "phys->ks" are in [N/m] units [1].
In my case, in the KnKsPBLaw.cpp the "phys->kn" and "phys->ks" are expressed in [N/m^3], while when we multiply them with the contact area, we get the stiffness parameters: "phys->Knormal_area" and "phys->Kshear_area", which are in [N/m] units [2, 3].

I would like your opinion whether I could modify the GlobalStiffnessTimeStepper to consider the "phys->Knormal_area" instead of the "phys->Kn" (and respectively for the shear), when the particle shape is "PotentialBlock" instead of "Sphere". Do you think this would slow down the engine or find any theoretical drawbacks?



Question information

English Edit question
Yade Edit question
No assignee Edit question
Solved by:
Bruno Chareyre
Last query:
Last reply:
Robert Caulk (rcaulk) said : #1

> Do you think this would slow down the engine or find any theoretical drawbacks?

If it is just a single logical operator checking the for a potentialblock, the engine won't slow down. From a units standpoint, what you suggest makes more sense.

I guess I just have two questions:
  - what is currently happening when you GlobalSitffnessTimeStepper with potential blocks? Is it over estimating the time step, resulting in instability? If not, the change you suggest will simply decrease the estimated time step which will slow down yours (and maybe others) simulations.

 - are other people depending on the code being the way it is now? (I guess no, but you know better than me)


Janek Kozicki (cosurgi) said : #2

How about this:

1. reorganize a little bit GlobalStiffnessTimeStepper in following manner:
1.1. separate the code into another function which takes as an argument the [N/m] parameters, then returns the value.
1.2. The original code in GlobalStiffnessTimeStepper only calls this function with previously used [N/m] parameters

2. Then call the same function with your [N/m^3] values multiplied by contact area.

Now the question is if it is better to derive from class GlobalStiffnessTimeStepper, and call that function in parent class. Or maybe if it’s better to handle all cases already in GlobalStiffnessTimeStepper, distinguishing between the N/m and N/m^3.


Thank you both for your advice!


> Is it over estimating the time step, resulting in instability?
- Yes, currently the GlobalStiffnessTimeStepper can lead to instability, by overestimating the timestep. In particular, if we consider the cases where the contactArea > 1.0 (mm^2 or cm^2 or m^2), the timestep we currently get can be orders of magnitude larger than what it should be, if we considered the actual stiffness of the system.

> are other people depending on the code being the way it is now?
- I do not think anyone else currently uses the GlobalStiffnessTimeStepper with the Potential Blocks, so this shouldn't create any backward compatibility problems. If someone used it, we would have heard about this problem by now :)

Jan, thanks for the step-by step advice

> separate the code into another function
This is a good point, in order to not slow down the engine for those who use it with spheres.
I only have one reservation: If we separate the code into another function, we will be able to use the GlobalStiffnessTimeStepper efficiently only for models with one contact physics type. Now, in my case, the Potential Blocks only interact with themselves, so this is not a problem for me. But, if we wanted to use the GlobalStiffnessTimeStepper in a model with: Spheres+Polyhedra (for which the stiffnesses "phys->kn" and "phys->ks" are also expressed in different units than [N/m]), we would need to deal with both shapes in the same timestepper. I think an integrated solution would be to have one function, in which we iterate for each interaction, and depending on the type of the physics involved, calculate the stiffness of the interaction.


We got a similar situation some time ago with Hertzian models. Their tangent stiffness (kn=dfn/dun) depends on the normal force itself, and it's actually not needed when calculating the contact force.

The solution was as follows, by far the simpler:
- reserve attributes "kn" and "ks" to reflect the derivatives dfn/dun and dft/dut (can be updated by the law functor) to be used by GSTS, and compute them just in order to determine stability condition.
- give other names to any other stiffness used when actually computing the contact forces.

I would suggest the same approach. This way GSTS will work as is with every contact model, provided that kn, ks always have the same meaning.

If would start to introduce conditionals in GSTS to specialize it, OTOH, it will soon be a big mess.


In other words:
1/ rename the current kn, ks [force/volume] differently (automatic search replace), say volKn and volKs
2/ introduce new attributes named kn, ks, and assign relevant values to them. They'll be picked up by GSTS.

Hi Bruno,

Many thanks for your reply.
This indeed sounds simpler. The idea of having the same units for kn, ks in every contact model sounds more straightforward too.
I will update the KnKsPBLaw.cpp script to do that.


Thanks Bruno Chareyre, that solved my question.

Janek Kozicki (cosurgi) said : #8

Great solution, Bruno! :)