# Velocity gradient

Hi (Bruno),

I see some inconsistencies in Newton Integrator which involve calculation the velocity gradient. I am going to fix them but I would need your approval. The problems are the following ones.

1) line 120. vel is at (t-dt/2) but velGrad is at (t+dt/2)

2) line 178. prevVelGrad is at (t-dt/2) but vel is at (t+dt/2)

I think it would be better to be consistent. As for 1) I would consider prevVelGrad so it is expressed at time (t-dt/2). As for 2) everything should be expressed at time t.

Would you agree to change the code according to the above? Or would you have an explanation why it should not be changed?

Thanks a lot,

Chiara

## Question information

- Language:
- English Edit question

- Status:
- Answered

- For:
- Yade Edit question

- Assignee:
- No assignee Edit question

- Last query:
- 2011-09-02

- Last reply:
- 2011-09-05

I will have to think twice...

Is it related to the problem that Vaclav mentioned for homo=3?

Sorry, you lost me.

line 120 in Newton.cpp, I see:

scene->

line 178:

} YADE_PARALLEL_

I just synced with bzr.

Sorry! I got the wrong numbers. So here they are.

Lines 135:

Vector3r fluctVel=

->state-

and line 193:

if (homoDeform=

It is not related to the problem of homoDeform, that was a mistake in the

code not in the formulation. So what do you think?

Chiara

On 31 August 2011 17:55, Chareyre <email address hidden>wrote:

> Your question #169689 on Yade changed:

> https:/

>

> Chareyre proposed the following answer:

> Sorry, you lost me.

>

> line 120 in Newton.cpp, I see:

> scene->

>

> line 178:

> } YADE_PARALLEL_

>

> I just synced with bzr.

>

> --

> If this answers your question, please go to the following page to let us

> know that it is solved:

> https:/

>

> If you still need help, you can reply to this email or go to the

> following page to enter your feedback:

> https:/

>

> You received this question notification because you asked the question.

>

1) I agree to use prevVelGrad.

2) Indeed it's using v(t+dt/2) (more or less since it's not fully

updated). I think it should be using the value v(t-dt/2). How would it

read at time "t"?

I imagine this cycle where the body moves only due to cell deformation

(no contact forces), with velGrad switching between -0.5 and +0.5 at

each cycle, and assuming dt=1:

* position=1, velgrad=-0.5 : position goes to 0.5 after the first

integration, as velocity is -0.5

* next step, velgrad=+0.5 : velocity is adjusted at line 193:

v=v+0.25=-0.25 (correct, it corresponds to the velocity when the

particle arrives at position 0.5 in the old velocity field); then it's

accelerated the other way by dVelGrad=+2, v=0.25, pos=0.75.

The result after 2 steps looks correct (the position is not back to 1,

this is expected as the final size of the period itself is really 0.75,

not 1 (1*0.5*1.5)).

I commited a fix for 2), tell me if you agree (just moving a few lines

and updating acceleration instead of velocity).

On 1 September 2011 11:30, Chareyre <email address hidden>wrote:

> Your question #169689 on Yade changed:

> https:/

>

> Status: Open => Answered

>

> Chareyre proposed the following answer:

> 1) I agree to use prevVelGrad.

>

Ok, let's do it then. We should simple add the velGrad argument to that

function and then input the prevVelGrad value when it is called in the

integrator.

>

> 2) Indeed it's using v(t+dt/2) (more or less since it's not fully

> updated). I think it should be using the value v(t-dt/2). How would it

> read at time "t"?

>

Well the correction to the velocity due to the homogeneous deformation is

(was, before last changes) made up of two terms.

First term (this is at time t, so it is fine): Vector3r dVel=dVelGrad*state

->pos;

Second term (it should have been at time t but it is/was not): state->vel+=

scene->

The second term could be expressed like

(prevVelGrad+

I shall look at what you propose and I let you know. You also let me know

what you think about it.

I imagine this cycle where the body moves only due to cell deformation

> (no contact forces), with velGrad switching between -0.5 and +0.5 at

> each cycle, and assuming dt=1:

>

> * position=1, velgrad=-0.5 : position goes to 0.5 after the first

> integration, as velocity is -0.5

> * next step, velgrad=+0.5 : velocity is adjusted at line 193:

> v=v+0.25=-0.25 (correct, it corresponds to the velocity when the

> particle arrives at position 0.5 in the old velocity field); then it's

> accelerated the other way by dVelGrad=+2, v=0.25, pos=0.75.

>

> The result after 2 steps looks correct (the position is not back to 1,

> this is expected as the final size of the period itself is really 0.75,

> not 1 (1*0.5*1.5)).

> I commited a fix for 2), tell me if you agree (just moving a few lines

> and updating acceleration instead of velocity).

>

> --

> If this answers your question, please go to the following page to let us

> know that it is solved:

> https:/

>

> If you still need help, you can reply to this email or go to the

> following page to enter your feedback:

> https:/

>

> You received this question notification because you asked the question.

>

For 1), you can also use cell->prevVelGrad in the function and

everywhere in Newton (and update it). This attribute is currently useless.

For 2), I think the two corrections should really be at different times.

The convective term is a correction of velocity accounting for the

displacement that occured between t-dt and t. This displacement is

exactly v(t-dt/2)*dt.

The term with dVelGrad reflects an "impulse" (or acceleration) occuring

at time t, yes, but it actually defines the velocity between t and t+dt

(what we call v(t+dt/2).

In other words, we correct the old, then we define the new.

Computing an average velGrad for time t would not make much sense since

velocities are never defined for integer steps.

If you apply this idea in the cycle I described, if I'm not wrong, it

would result in:

* first step: v=v-0.25*0=0 (convective); v=dVelGrad*pos=-0.5 (like before)

* second step: v=v+0*meanVel=-0.5 (no correction); v=v+dVelGrad*pos=0

As you see the particle is stuck at second step. The cell is deforming

but the particle doesn't move.

I'm not really sure I'm applying your suggestion correctly because in

fact this equation looks recursive to me:

v+=(prevVelGrad

It means:

v(t+dt/

How can we use v(t+dt/2) to define v(t+dt/2)?

Same remark for a(t), which is unknown before we apply this correction.

On 1 September 2011 17:00, Chareyre <email address hidden>wrote:

> Your question #169689 on Yade changed:

> https:/

>

> Status: Open => Answered

>

> Chareyre proposed the following answer:

> For 1), you can also use cell->prevVelGrad in the function and

> everywhere in Newton (and update it). This attribute is currently useless.

>

> For 2), I think the two corrections should really be at different times.

> The convective term is a correction of velocity accounting for the

> displacement that occured between t-dt and t. This displacement is

> exactly v(t-dt/2)*dt.

> The term with dVelGrad reflects an "impulse" (or acceleration) occuring

> at time t, yes, but it actually defines the velocity between t and t+dt

> (what we call v(t+dt/2).

> In other words, we correct the old, then we define the new.

>

> Computing an average velGrad for time t would not make much sense since

> velocities are never defined for integer steps.

>

In fact I was considering to compute the two corrections to the velocity

both at time t. As you say, we use dVelGrad to update velocity at mid-step.

For me it was correct to consider _both_ corrections at time t, being the

value we update at time (t+dt/2). Anyway there is the problem below where I

would use v(t+dt/2) to update v(t+dt/2). In fact the latter is not really

v(t+dt/2) since this is not yet updated.

I can only partially see what you are doing by splitting these two

contributions, one acting on the acceleration and the other one on the

velocity. Reading the chapter by Radjai in the Discrete Element book about

PBCs, I can see that they apply ONE single corrective term to the

acceleration (please see eq. 7.26 at pag. 190). However there is a

difference between their formulation and ours since they distinguish the

fluctuating velocity of the fluid particle to the velocity field applied

whereas we do not do that if not when we compute kinetic energy.

Can we somehow be consistent with their formulation which seems good to me?

Chiara

> If you apply this idea in the cycle I described, if I'm not wrong, it

> would result in:

>

> * first step: v=v-0.25*0=0 (convective); v=dVelGrad*pos=-0.5 (like before)

> * second step: v=v+0*meanVel=-0.5 (no correction); v=v+dVelGrad*pos=0

>

> As you see the particle is stuck at second step. The cell is deforming

> but the particle doesn't move.

> I'm not really sure I'm applying your suggestion correctly because in

> fact this equation looks recursive to me:

>

> v+=(prevVelGrad

>

>

> It means:

>

> v(t+dt/

>

>

> How can we use v(t+dt/2) to define v(t+dt/2)?

> Same remark for a(t), which is unknown before we apply this correction.

>

> --

> If this answers your question, please go to the following page to let us

> know that it is solved:

> https:/

>

> If you still need help, you can reply to this email or go to the

> following page to enter your feedback:

> https:/

>

> You received this question notification because you asked the question.

>

For me eq. 7.26 is very close to:

Vector3r dVel=dVelGrad*

state->vel+=dVel;

dVel is an impulse at time t right?

Now the other correction we have is either hidden in their equations,

due to the difference between fluctuational and total velocity, or due a

mistake in a derivation on one side.

I wouldn't bet my life on that but when I checked I found a mistake in

their derivation. I'll try and find it again and will let you know.

> Now the other correction we have is either hidden in their equations,

> due to the difference between fluctuational and total velocity, or due a

> mistake in a derivation on one side.

The first is true. In Radjai, the contact laws must compute the absolute

velocity = fluctuation+

In Yade we assign the absolute velocity in Newton, so that in contact

laws we don't care summing two terms.

In turn we have to add velG*pos and remove it at the next time-step.

You can see that by rewriting the operations differently. We have now:

v+=VG(t-

v+=(VG(

Due the fact that v(t-1/2)

v+= -VG(t-1/2)*pos(t-1) + VG(t+1/2)*pos(t)

i.e. we remove the affine velocity of the previous step from the

absolute velocity (it leaves the fluctuation only) and we add the affine

velocity of the next step.

We are exactly equivalent I think.

When you say Radjai applies only one correction, it is not true. He has

only one correction in Newton, but he has another one in contact laws.

We have two in Newton but none in contact laws.

Since one body's velocity will be used for many contacts, it seems more

optimal to compute it once in Newton. I don't think we have to change

anything.

## Can you help with this problem?

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