Non-hourglass formulation for update Lagrange method

[Shuaihao-Zhang]

I have considered the conservation of angular momentum, but there are still some problems. I will describe my problems here.

The shear acceleration can be express as this:

When considering 2D situation, d = 2, it becomes this:

For previous formulation with hourglass (i.e., the oscillating beam case I uploaded previously), the result of oscillating beam is this: (colored by shear acceleration)

With non-hourglass formulation, when the conservation of angular momentum was not considered (the result I showed you last week), the result is like this https://github.com/Xiangyu-Hu/SPHinXsys/blob/e7336713a388d39f8e9016f236a1481b4d19fa0a/tests/user_examples/ExtraSources/extra_sources_share/granular/granular_dynamics/granular_dynamics_inner.cpp#L33-L56 (colored by shear acceleration) We can see that the beam can not move upwards.

With non-hourglass formulation, when the conservation of angular momentum was considered, the result is like this: https://github.com/Xiangyu-Hu/SPHinXsys/blob/e7336713a388d39f8e9016f236a1481b4d19fa0a/tests/user_examples/ExtraSources/extra_sources_share/granular/granular_dynamics/granular_dynamics_inner.cpp#L57-L78 (colored by shear acceleration) The beam move upwards too much.

[Shuaihao-Zhang]

Sorry, I found a bug in my code when calculating shear acceleration considering the conservation of angular momentum. I will check the result and update the result later.

[Xiangyu-Hu]

is there the same bug for the angular momentum non-conservative form?

[Shuaihao-Zhang]

No, for the angular momentum non-conservative form, the code is correct.

For the angular momentum conservative form, I guess I made some mistakes when coding this equation in your paper (2006, POF).

[Xiangyu-Hu]

My may consider to insert rigid body rotation to the formulation and derive shear acceleration.

[Shuaihao-Zhang]

Yeah. I haven't consider the rigid body rotation yet. Calculating the rigid body rotation needs the shear stress. But we skip the calculation of shear stress. I will try to calculate the shear stress according to the previous method and see how's the result.

[Shuaihao-Zhang]

Fix a bug in calculating the shear acceleration when considering the angular momentum conservative form. https://github.com/Xiangyu-Hu/SPHinXsys/blob/de930d1d4d1afe739353b1418a8af68de60c2405/tests/user_examples/ExtraSources/extra_sources_share/granular/granular_dynamics/granular_dynamics_inner.cpp#L60-L77

Although still not correct, but the result seems better.

It seems that the total momentum is not conservative.

Maybe because I didn't consider rigid body rotation, or maybe because the time step. I used a dt for integral acceleration. https://github.com/Xiangyu-Hu/SPHinXsys/blob/de930d1d4d1afe739353b1418a8af68de60c2405/tests/user_examples/ExtraSources/extra_sources_share/granular/granular_dynamics/granular_dynamics_inner.cpp#L76

[Shuaihao-Zhang]

My may consider to insert rigid body rotation to the formulation and derive shear acceleration.

I will try to incorporate the rigid body rotation into the formulation of shear acceleration.

[Xiangyu-Hu]

The stress profile is quite regular. No far from good results.

[Xiangyu-Hu]

I have considered the conservation of angular momentum, but there are still some problems. I will describe my problems here.

The shear acceleration can be express as this:

When considering 2D situation, d = 2, it becomes this:

For previous formulation with hourglass (i.e., the oscillating beam case I uploaded previously), the result of oscillating beam is this: (colored by shear acceleration)

With non-hourglass formulation, when the conservation of angular momentum was not considered (the result I showed you last week), the result is like this

https://github.com/Xiangyu-Hu/SPHinXsys/blob/e7336713a388d39f8e9016f236a1481b4d19fa0a/tests/user_examples/ExtraSources/extra_sources_share/granular/granular_dynamics/granular_dynamics_inner.cpp#L33-L56

(colored by shear acceleration) We can see that the beam can not move upwards. With non-hourglass formulation, when the conservation of angular momentum was considered, the result is like this:

https://github.com/Xiangyu-Hu/SPHinXsys/blob/e7336713a388d39f8e9016f236a1481b4d19fa0a/tests/user_examples/ExtraSources/extra_sources_share/granular/granular_dynamics/granular_dynamics_inner.cpp#L57-L78

(colored by shear acceleration) The beam move upwards too much.

I feel that the derivation from line 3 to the last line is not clear.

[Shuaihao-Zhang]

This is the detailed derivation process. You could help me check this when you are available.

[Shuaihao-Zhang]

Maybe I mixed the row matrix and column matrix. The third item should be a column matrix, and it looks like this. Not sure about this. There are no matrix dimension error when I use the previous row matrix.

[Shuaihao-Zhang]

But anyway,

This means there is only one item is the shear acceleration.

[Xiangyu-Hu]

Maybe I mixed the row matrix and column matrix. The third item should be a column matrix, and it looks like this. Not sure about this. There are no matrix dimension error when I use the previous row matrix.

Therefore, this term can not be canceled with the last divergence term in the original derivation.

[Xiangyu-Hu]

Another possibility for the increasing magnitude may be due the inconsistency of the displacement and the shear stress. One may consider to use incremental (dr_ij e_ij)e_ij other than use (v_ij e_ij) dt * e_ij, because dr_i is actually updated by two sub-steps with updating of r_i.

[Shuaihao-Zhang]

Another possibility for the increasing magnitude may be due the inconsistency of the displacement and the shear stress. One may consider to use incremental (dr_ij e_ij)e_ij other than use (v_ij e_ij) dt * e_ij, because dr_i is actually updated by two sub-steps with updating of r_i.

Thanks. I will try this. I have verified that the result is not relevant to the limiter in the Riemann solver. Because I actually used the DissipativeRiemannSolver, there was no limiter already.

[Shuaihao-Zhang]

Now the stain tensor can be calculated correctly. Previous, when we calculate the velocity gradient, we didn't exclude the rotation. I modified the form to: Then, the result seems good. The same thing for the rotation tensor.

Not very sure if this modification is physically correct.

[Shuaihao-Zhang]

I found two problems when using dual-criterion timestep for solids.

When simulating the oscillating beam: (1) The result cannot converge when using dual-criterion timestep, as shown in the below: I guess there are some errors when using dual-criterion timestep. So I tested the convergence for single timestep, the result is ok (as shown below). I think the result for dual-timestep can be accepted, if the the result of dual-criterion timestep can be consistent with that of single timestep when increasing the resolution. However, it doesn't, which leads to the second problem.

(2)The results of dual-criterion timestep are not consistent with those of the single timestep, even increasing the resolution. Yesterday, I showed a consistent result between dual-timestep and single-timestep. As shown below: But the results are wrong, because I used a wrong anticipated maximum velocity. I just used the coefficient 0.05, but it's actually 2.8. So the results maybe not consistent, they just happen to appear consistent.

Anyway, when I used a correct anticipated maximum velocity to re-simulate, the results are bad (as shown below).

I thought about the reason. The advection timestep is 10-20 times as much as the acoustic timestep. I guess maybe the advection timestep should be reduced.

[Shuaihao-Zhang]

Another thing is about using B correction when calculating acceleration. I used this formulation and found the result is not totally correctly. The amplitude is unexpected large and the period is also large compare with previous result. As shown below, the colored one is the result using B correction for calculating the acceleration, and the grey one in the result from previous method. Although the amplitude and period are not correct, the mechanical energy is still conservative. I think there also should be a coefficient when calculating the the second-order derivative for velocity: I found that, with a coefficient greater than 1, the amplitude and period can be reduced.

[DongWuTUM]

Another thing is about using B correction when calculating acceleration. I used this formulation and found the result is not totally correctly. The amplitude is unexpected large and the period is also large compare with previous result. As shown below, the colored one is the result using B correction for calculating the acceleration, and the grey one in the result from previous method. Although the amplitude and period are not correct, the mechanical energy is still conservative. I think there also should be a coefficient when calculating the the second-order derivative for velocity: I found that, with a coefficient greater than 1, the amplitude and period can be reduced.

Please check if correction matrix B is placed in right place, and If the transpose of B is used.

[Shuaihao-Zhang]

Please check if correction matrix B is placed in right place, and If the transpose of B is used.

I think the B_ is applied at the correct place. The results are the same when I use the transpose of B_, and the amplitude is still too large.