# RungeKuttaCashKarp54Integrator with gridConnections is causing errors

Hello,

I want to use the RungeKuttaCashK

The minimum working code is given below. Please help me to fix the issue

Kind regards,

Rohit John

#######

# encoding: utf-8

"An example showing bending beams."

from yade.gridpfacet import *

#### Parameter ####

L=10. # length of the beam

n=12 # number of nodes used to generate the beam

r=L/50. # radius of the beam element

#### Engines ####

integrator=

ForceResetter(),

GeneralIntegra

InteractionLoop(

[Ig2_

[Ip2_

[Law2_

),

GravityEngine(

])

#Tolerances can be set for the optimum accuracy

integrator.

integrator.

O.engines=

]

#### Create materials and set different properties ####

O.materials.

#### Vertices ####

vertices1=[]

for i in range(0,n):

vertices1.append( [i*L/n,0,0] )

#### Create cylinder connections ####

nodesIds=[]

cylIds=[]

cylinderConnect

#### Set boundary conditions ####

for i in range(0,1):

O.bodies[

#### Set a time step ####

O.dt=1e-05

#### Allows to reload the simulation ####

O.saveTmp()

## Question information

- Language:
- English Edit question

- Status:
- Solved

- For:
- Yade Edit question

- Assignee:
- No assignee Edit question

- Solved by:
- Janek Kozicki

- Solved:
- 2021-01-28

- Last query:
- 2021-01-28

- Last reply:
- 2021-01-28

Janek Kozicki (cosurgi) said : | #1 |

See: https:/

> I added viscous damping to the source code of Law2_ScGeom6D_

Quick note:

We can't really try since we don't have your specific damping model. You could push to a gitlab branch at some point if you want someone to be able to try the script with relevant code version.

Maybe you did that already, but before trying another integrator (or higher number precision ;) ), I would go through this checkList:

- what happens with super-small O.dt? (defining what is "super-small" needs to think about mass-viscosity-

- are the directions of forces, moment, velocities and spin consistent, with just 2-3 nodes after 2 iterations?

Cheers

Bruno

Rohit John (rohitkjohn) said : | #3 |

Thanks Janek Kozicki, that solved my question.

Rohit John (rohitkjohn) said : | #4 |

Dear Bruno,

> what happens with super-small O.dt? (defining what is "super-small" needs to think about mass-viscosity-

- It depends on the damping fraction and the material stiffness. I wrote the code such that the damping coefficient is a fraction of the stiffness. (Same fraction for normal, shear, twist and bending stiffness in Law2_ScGeom6D_

> are the directions of forces, moment, velocities and spin consistent, with just 2-3 nodes after 2 iterations?

I checked the values and when the damping and dt were badly set, the forces and velocities grow exponentially, but when I reduce the dt or damping they "calm down".

From my understanding the integrator in YADE is similar to forward Euler which is unstable when there is damping. So I though changing the integrator might help. Janek Kozicki solution works but the RungeKuttaCashK

If possible could you check could the code found here: https:/

Kind regards,

Rohit K. John

You are past-the-end of my checklist. :)

RungeKuttaCashK

Last one maybe: in a working case (small damping), did you find consistency _at the structural scale_ in terms of damping the oscillations? It could be that there is a mismatch between what you input at the interactions scale and what you expect in terms of damping at the structural scale?

> is similar to forward Euler

Yes and no. It is a 2nd order time centered scheme from Newton's point of view. Problem is velocities are not known on step, they are know at mid-steps. So this wonderful scheme for mass-spring systems becomes a poor 1st order Euler when there are viscous terms. If you take mid-steps as your frame of reference and you think to what happens to velocities wrt. visous forces it is plain explicit Euler, you are right.

There are ways to evaluate on-step velocites (something similar to "velocity Verlet" [1]) which we don't use. It should improve your case, it needs some changes in the integrator and it needs two calculations of forces per time iteration.

> the RungeKuttaCashK

Because there are even more than two calculations of forces per time iteration, I guess?

Do you mean that despite larger timesteps it is still slower?

Bruno

> It could be that there is a mismatch between what you input at the interactions scale and what you expect in terms of damping at the structural scale

I would say the bending viscosity of the interactions should be a multiple of ~N/L, do we agree on that?

> https:/

Problem here is we can't see you changes.

From the top of my head, you need to do something like:

0. get sources with "git clone"

1. edit sources

2. "git commit -a -m 'that's my change'"

3. git push to gitlab/github whatever

This way we can see what was done specifically in step 1

Bruno

Rohit John (rohitkjohn) said : | #7 |

Dear Bruno,

> Problem here is we can't see you changes.

- I see. I tried pushing the code into the YADE gitlab and I got this error

remote: You are not allowed to push code to this project.

So I pushed the source into my account. You can view the changes here.

https:/

Kind regards,

Rohit K. John

Rohit John (rohitkjohn) said : | #8 |

Dear Bruno,

>Last one maybe: in a working case (small damping), did you find consistency _at the structural scale_ in terms of damping the oscillations? It could be that there is a mismatch between what you input at the interactions scale and what you expect in terms of damping at the structural scale?

- I cannot say this with certainity because I did not find any literature with data that I can use to verify me results. But in the structural case I did find that the motion is damped. For example, in one simulation I had attached at weight to the tip of a cantilever (a forceEngine at the particle at the end) and let it reach equilibrium after deformation. Then the weight was released (forceEngine dead = True), the oscillations were similar to what we find in a damped oscillator.

The mismatch that you mentions is a bit unclear. Could you please elaborate or maybe give an example?

> I would say the bending viscosity of the interactions should be a multiple of ~N/L, do we agree on that?

- Could you please define N and L. Am I wrong to assume N is the number of gridNodes and L is the length of the cantilever? My first thought is that damping should be a material property so it should be constant. However, it could also depend on the shape of the cantilever, like the area moment of inertia perhaps. It is difficult to find good literature on damping in beams. Or I might looking in the wrong place or the wrong thing.

Kind regards,

Rohit K. John

Hi,

> the oscillations were similar to what we find in a damped oscillator

Cool. :)

> Am I wrong to assume N is the number of gridNodes and L is the length of the cantilever?

You are right.

Sorry, I wrote that comment right after reading another thread on a linked topic, where N and L were clearly defined. Overall, the thing is: interaction properties - if they are supposed to reflect a beam-like thing - necessarily depend on beam geometry and discretisation. It can't be conveyed by material properties alone. The "mismatch" I was mentioning was about that: material properties are not "per-unit length of a beam" properties and so someone can be off by orders of magnitudes if based on wrong expectations.

I couldn't check the github push in depth, yet, but thanks for clarity.

Let you report here if you have more insights on the issue.

Bruno

Rohit John (rohitkjohn) said : | #10 |

Thanks for your advice Bruno, I shall put an update in this thread when I find something useful.

Kind regards,

Rohit K. John