draft PR -- hourglass control (updated Lagrangian) for elasticity and plasticity

[Shuaihao-Zhang]

Here, I will show the results when dealing with the hourglass control method for ULSPH.

[Shuaihao-Zhang]

This is the method we discussed just now. We only calculate the inverse of shear strain for particle i when the return mapping is called, i.e., when particle i is plastic. Then, we can calculate the scale matrix. If particle iis elastic, the scale matrix is the identity matrix.

https://github.com/Xiangyu-Hu/SPHinXsys/blob/06d36e86df5301e93e3381d6e9ed49e9a845ff88/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp#L545-L553

The result is as follows:

As Prof. Hu mentioned, if the particle is plastic, then the inverse of shear strain exists.

However, there are some other problems, and the simulation crash after several steps. And it looks like tensile instability.

I think this may because the scale matrix is too small at the bottom of the bar, as shown in the following figure. The color indicates the scale matrix. It's an identity matrix when the particle is elastic, and the magnitude is 1.4 in 2D. If the particle is plastic, the magnitude of the scale matrix is smaller than 1.4.

I suppose, if the scale matrix is very small, the effect of this hourglass control method is also very small. So hourglass and tensile instability occurs.

The figure below shows the plastic indicator.

[Xiangyu-Hu]

Are you sure that you have 1/2 in front of (\phi_i + \phi_j)?

[Xiangyu-Hu]

So the parameter \xi has not effect to the solution?

[Xiangyu-Hu]

This is the method we discussed just now. We only calculate the inverse of shear strain for particle i when the return mapping is called, i.e., when particle i is plastic. Then, we can calculate the scale matrix. If particle iis elastic, the scale matrix is the identity matrix.

https://github.com/Xiangyu-Hu/SPHinXsys/blob/06d36e86df5301e93e3381d6e9ed49e9a845ff88/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp#L545-L553

The result is as follows:

As Prof. Hu mentioned, if the particle is plastic, then the inverse of shear strain exists.

However, there are some other problems, and the simulation crash after several steps. And it looks like tensile instability.

I think this may because the scale matrix is too small at the bottom of the bar, as shown in the following figure. The color indicates the scale matrix. It's an identity matrix when the particle is elastic, and the magnitude is 1.4 in 2D. If the particle is plastic, the magnitude of the scale matrix is smaller than 1.4.

I suppose, if the scale matrix is very small, the effect of this hourglass control method is also very small. So hourglass and tensile instability occurs.

The figure below shows the plastic indicator.

Also, in the first part of the equation should it be \tild{\sigma_i} + \tild{\sigma_j}?

[Shuaihao-Zhang]

Are you sure that you have 1/2 in front of (\phi_i + \phi_j)?

Yes, I do put 1/2 in the code. This actually doesn't matter because we have this coefficient in front of this integration.

[Shuaihao-Zhang]

So the parameter \xi has not effect to the solution?

The parameter 𝜉 do have effects. But the simulation still crash even I used a large 𝜉 . I think this is because the difference between the elastic correction and plastic correction is too large, due to the scale matrix 𝝋 for plastic particle is too small.

[Shuaihao-Zhang]

Also, in the first part of the equation should it be \tild{\sigma_i} + \tild{\sigma_j}?

Yes, this should be the shear stress after return mapping. This is a typo here. The code is correct.

[Xiangyu-Hu]

I have checked you code and found, besides quite intervened coding, a serious consistent issue. That is. If I comment line https://github.com/Shuaihao-Zhang/SPHinXsys/blob/06d36e86df5301e93e3381d6e9ed49e9a845ff88/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp#L533 to skip return mapping and change https://github.com/Shuaihao-Zhang/SPHinXsys/blob/06d36e86df5301e93e3381d6e9ed49e9a845ff88/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp#L555 to Mat3d deviatoric_stress = J2plasticity.ReturnMappingShearStress(shear_stress3D[index_i]);

Then the scaling factor should be Identity even we are applying https://github.com/Shuaihao-Zhang/SPHinXsys/blob/06d36e86df5301e93e3381d6e9ed49e9a845ff88/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp#L548

However, the simulation still gives crashed results. Please have a check.

[Xiangyu-Hu]

Also. in the TL version, we use full strain to obtain scaling factor. Please have a look. https://github.com/Xiangyu-Hu/SPHinXsys/blob/bbe70dc9411ab91c98b5c08fa3d8848b65112bef/src/shared/particle_dynamics/solid_dynamics/inelastic_dynamics.cpp#L44-L47C49

[Shuaihao-Zhang]

I have checked you code and found, besides quite intervened coding, a serious consistent issue. That is. If I comment line https://github.com/Shuaihao-Zhang/SPHinXsys/blob/06d36e86df5301e93e3381d6e9ed49e9a845ff88/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp#L533 to skip return mapping and change https://github.com/Shuaihao-Zhang/SPHinXsys/blob/06d36e86df5301e93e3381d6e9ed49e9a845ff88/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp#L555 to Mat3d deviatoric_stress = J2plasticity.ReturnMappingShearStress(shear_stress3D[index_i]);

Then the scaling factor should be Identity even we are applying https://github.com/Shuaihao-Zhang/SPHinXsys/blob/06d36e86df5301e93e3381d6e9ed49e9a845ff88/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp#L548

However, the simulation still gives crashed results. Please have a check.

(1) Sorry for the coding. I have deleted some unnecessary code blocks and comments to make it more clear. (2) The plastic behavior is not only controlled by return mapping, but also by the constitutive model. If we want to transfer the material to elasticity, besides skip the return mapping, we also need to modify a line of code in the constitutive model. https://github.com/Xiangyu-Hu/SPHinXsys/blob/2eac06a0ef5745f0c97d47e71832d9114a21e9ef/tests/user_examples/extra_src/shared/general_continuum/general_continuum.cpp#L92-L114 We just change the Line-112 to return shear_stress_rate_elastic;, then we can get this elastic behavior. And the scale matrix is an identity matrix for all particles.

[Shuaihao-Zhang]

Also. in the TL version, we use full strain to obtain scaling factor. Please have a look. https://github.com/Xiangyu-Hu/SPHinXsys/blob/bbe70dc9411ab91c98b5c08fa3d8848b65112bef/src/shared/particle_dynamics/solid_dynamics/inelastic_dynamics.cpp#L44-L47C49

I will check this in detail tomorrow.

[DongWuTUM]

@Shuaihao-Zhang, the returning map you wrote is to return the deviatoric part. The inconsistency between the elastic and plastic behaviors may arise due to the incompressible assumption. Pls, check it. Could you write down the plastic constitutive relation formulation in detail?

[DongWuTUM]

@Xiangyu-Hu @Shuaihao-Zhang Form the derivation, we should use the shear strain to obtain the scaling factor. But we can have a try for full strain.

[Xiangyu-Hu]

Further more, the TL version use the full strain tensor not only the shear tension for the scaling factor. Please see

I have checked you code and found, besides quite intervened coding, a serious consistent issue. That is. If I comment line https://github.com/Shuaihao-Zhang/SPHinXsys/blob/06d36e86df5301e93e3381d6e9ed49e9a845ff88/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp#L533 to skip return mapping and change https://github.com/Shuaihao-Zhang/SPHinXsys/blob/06d36e86df5301e93e3381d6e9ed49e9a845ff88/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp#L555 to Mat3d deviatoric_stress = J2plasticity.ReturnMappingShearStress(shear_stress3D[index_i]); Then the scaling factor should be Identity even we are applying https://github.com/Shuaihao-Zhang/SPHinXsys/blob/06d36e86df5301e93e3381d6e9ed49e9a845ff88/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp#L548 However, the simulation still gives crashed results. Please have a check.

(1) Sorry for the coding. I have deleted some unnecessary code blocks and comments to make it more clear. (2) The plastic behavior is not only controlled by return mapping, but also by the constitutive model. If we want to transfer the material to elasticity, besides skip the return mapping, we also need to modify a line of code in the constitutive model.

https://github.com/Xiangyu-Hu/SPHinXsys/blob/2eac06a0ef5745f0c97d47e71832d9114a21e9ef/tests/user_examples/extra_src/shared/general_continuum/general_continuum.cpp#L92-L114

We just change the Line-112 to return shear_stress_rate_elastic;, then we can get this elastic behavior. And the scale matrix is an identity matrix for all particles.

Sorry. If I do as you said the

I have checked you code and found, besides quite intervened coding, a serious consistent issue. That is. If I comment line https://github.com/Shuaihao-Zhang/SPHinXsys/blob/06d36e86df5301e93e3381d6e9ed49e9a845ff88/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp#L533 to skip return mapping and change https://github.com/Shuaihao-Zhang/SPHinXsys/blob/06d36e86df5301e93e3381d6e9ed49e9a845ff88/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp#L555 to Mat3d deviatoric_stress = J2plasticity.ReturnMappingShearStress(shear_stress3D[index_i]); Then the scaling factor should be Identity even we are applying https://github.com/Shuaihao-Zhang/SPHinXsys/blob/06d36e86df5301e93e3381d6e9ed49e9a845ff88/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp#L548 However, the simulation still gives crashed results. Please have a check.

(1) Sorry for the coding. I have deleted some unnecessary code blocks and comments to make it more clear. (2) The plastic behavior is not only controlled by return mapping, but also by the constitutive model. If we want to transfer the material to elasticity, besides skip the return mapping, we also need to modify a line of code in the constitutive model.

https://github.com/Xiangyu-Hu/SPHinXsys/blob/2eac06a0ef5745f0c97d47e71832d9114a21e9ef/tests/user_examples/extra_src/shared/general_continuum/general_continuum.cpp#L92-L114

We just change the Line-112 to return shear_stress_rate_elastic;, then we can get this elastic behavior. And the scale matrix is an identity matrix for all particles.

If I do as you said by: return shear_stress_rate_elastic; Then I do not comments line 533 but move it before, Mat3d deviatoric_stress = shear_stress3D[index_i]; The simulation now should not blowup. but it does. So, there seems still has consistent issue.

[Xiangyu-Hu]

Because it does not blowup if we simply constrain scale factor to identity in any case.

[Shuaihao-Zhang]

If I do as you said by: return shear_stress_rate_elastic; Then I do not comments line 533 but move it before, Mat3d deviatoric_stress = shear_stress3D[index_i]; The simulation now should not blowup. but it does. So, there seems still has consistent issue.

To transfer to elasticity, we should switch to return shear_stress_rate_elastic; and also comment Line 533 shear_stress_3D_[index_i] = J2_plasticity_.ReturnMappingShearStress(shear_stress_3D_[index_i]); to skip return mapping. If we don't, then the return mapping still works, which makes the stress-strain behavior is non-linear. The stress and strain after return mapping is very small, so the scale matrix 𝝋 will be small. Then effect the hourglass control is small, and the simulation crashes because (I guess) tensile instability, as I mention above.

[Shuaihao-Zhang]

@Shuaihao-Zhang, the returning map you wrote is to return the deviatoric part. The inconsistency between the elastic and plastic behaviors may arise due to the incompressible assumption. Pls, check it. Could you write down the plastic constitutive relation formulation in detail?

Yes, this is the constitutive for J2 plasticity. For J2 plasticity, the plastic part only related to the shear stress and the pressure part has nothing to do with the plasticity. This is why the return mapping only applied to the deviatoric part.

[Shuaihao-Zhang]

@Xiangyu-Hu @Shuaihao-Zhang Form the derivation, we should use the shear strain to obtain the scaling factor. But we can have a try for full strain.

@Xiangyu-Hu @DongWuTUM I have check the full strain, and there is a problem. For shear strain, we can get the shear strain form shear stress by dividing by 2G. While for total strain, we can not simply get the total strain from total stress.

[Xiangyu-Hu]

If I do as you said by: return shear_stress_rate_elastic; Then I do not comments line 533 but move it before, Mat3d deviatoric_stress = shear_stress3D[index_i]; The simulation now should not blowup. but it does. So, there seems still has consistent issue.

To transfer to elasticity, we should switch to return shear_stress_rate_elastic; and also comment Line 533 shear_stress_3D_[index_i] = J2_plasticity_.ReturnMappingShearStress(shear_stress_3D_[index_i]); to skip return mapping. If we don't, then the return mapping still works, which makes the stress-strain behavior is non-linear. The stress and strain after return mapping is very small, so the scale matrix 𝝋 will be small. Then effect the hourglass control is small, and the simulation crashes because (I guess) tensile instability, as I mention above.

I think that I am wrong here because the shear stress became plastic at next time step then the scale factor is not identity anymore. However, the scaling in TL version should have check.

[Xiangyu-Hu]

@Xiangyu-Hu @Shuaihao-Zhang Form the derivation, we should use the shear strain to obtain the scaling factor. But we can have a try for full strain.

@Xiangyu-Hu @DongWuTUM I have check the full strain, and there is a problem. For shear strain, we can get the shear strain form shear stress by dividing by 2G. While for total strain, we can not simply get the total strain from total stress.

You can have check the TL version. For the scaling factor, the method add the isotropic part back to the right hand cauchy strain.

[Xiangyu-Hu]

It is clear that if we included the isotropic part the variation of the scaling factor will be smaller.

[Shuaihao-Zhang]

It is clear that if we included the isotropic part the variation of the scaling factor will be smaller.

Sure, I will check this.

[Xiangyu-Hu]

Another thing is that the transformation between 2 and 3 d may introduce consistent issue.

[Xiangyu-Hu]

We better begin with the simplest plastic case first.

[Shuaihao-Zhang]

Another thing is that the transformation between 2 and 3 d may introduce consistent issue.

I think this is ok. We can use the case zsh_2d_oscillating_beam_plastic to test. The J2 plasticity material is used in this case. If we skip the return mapping and also switch to return shear_stress_rate_elastic; in the constitutive relative. We can get exactly the same result as the case test_2d_oscillating_beam_ul, which use the elastic material.

[Shuaihao-Zhang]

We better begin with the simplest plastic case first.

Sure, I will try to find some other cases, maybe one dimensional stretching or compression of metals.

[Xiangyu-Hu]

We better begin with the simplest plastic case first.

Sure, I will try to find some other cases, maybe one dimensional stretching or compression of metals.

No. As our issue is stability, we need just simply the material properties for testing.

[Shuaihao-Zhang]

We better begin with the simplest plastic case first.

Sure, I will try to find some other cases, maybe one dimensional stretching or compression of metals.

No. As our issue is stability, we need just simply the material properties for testing.

I guess J2 plasticity is one of the simplest plastic models. I will check.

[DongWuTUM]

If I do as you said by: return shear_stress_rate_elastic; Then I do not comments line 533 but move it before, Mat3d deviatoric_stress = shear_stress3D[index_i]; The simulation now should not blowup. but it does. So, there seems still has consistent issue.

To transfer to elasticity, we should switch to return shear_stress_rate_elastic; and also comment Line 533 shear_stress_3D_[index_i] = J2_plasticity_.ReturnMappingShearStress(shear_stress_3D_[index_i]); to skip return mapping. If we don't, then the return mapping still works, which makes the stress-strain behavior is non-linear. The stress and strain after return mapping is very small, so the scale matrix 𝝋 will be small. Then effect the hourglass control is small, and the simulation crashes because (I guess) tensile instability, as I mention above.

I think that I am wrong here because the shear stress became plastic at next time step then the scale factor is not identity anymore. However, the scaling in TL version should have check.

I think it is Okay that the scale factor is identity because we have checked whether it is the plastic stage or not just before calculating the scale factor. And it is just related to this time step, not the next step.

[DongWuTUM]

We better begin with the simplest plastic case first.

Sure, I will try to find some other cases, maybe one dimensional stretching or compression of metals.

No. As our issue is stability, we need just simply the material properties for testing.

I guess J2 plasticity is one of the simplest plastic models. I will check.

Yes, it is already the simplest.

[DongWuTUM]

If I do as you said by: return shear_stress_rate_elastic; Then I do not comments line 533 but move it before, Mat3d deviatoric_stress = shear_stress3D[index_i]; The simulation now should not blowup. but it does. So, there seems still has consistent issue.

To transfer to elasticity, we should switch to return shear_stress_rate_elastic; and also comment Line 533 shear_stress_3D_[index_i] = J2_plasticity_.ReturnMappingShearStress(shear_stress_3D_[index_i]); to skip return mapping. If we don't, then the return mapping still works, which makes the stress-strain behavior is non-linear. The stress and strain after return mapping is very small, so the scale matrix 𝝋 will be small. Then effect the hourglass control is small, and the simulation crashes because (I guess) tensile instability, as I mention above.

I am still confused about why "the stress and strain after return mapping are very small, so the scale matrix 𝝋 will be small". If it is perfect plasticity, the shear strain after the return map should be the same as the elastic strain. Does it mean the elastic strain is very small? Please check this.

[Shuaihao-Zhang]

If I do as you said by: return shear_stress_rate_elastic; Then I do not comments line 533 but move it before, Mat3d deviatoric_stress = shear_stress3D[index_i]; The simulation now should not blowup. but it does. So, there seems still has consistent issue.

To transfer to elasticity, we should switch to return shear_stress_rate_elastic; and also comment Line 533 shear_stress_3D_[index_i] = J2_plasticity_.ReturnMappingShearStress(shear_stress_3D_[index_i]); to skip return mapping. If we don't, then the return mapping still works, which makes the stress-strain behavior is non-linear. The stress and strain after return mapping is very small, so the scale matrix 𝝋 will be small. Then effect the hourglass control is small, and the simulation crashes because (I guess) tensile instability, as I mention above.

I am still confused about why "the stress and strain after return mapping are very small, so the scale matrix 𝝋 will be small". If it is perfect plasticity, the shear strain after the return map should be the same as the elastic strain. Does it mean the elastic strain is very small? Please check this.

Yeah, this may be the issue here. Here, we calculate the shear strain after return mapping based on the shear stress after return mapping. Initially, we use the way because we thought that, for elasticity, the scale matrix 𝝋 will be the identity matrix; while for plasticity, there is a reduction for the shear stress (due to return mapping and plastic constitutive relation), then the shear strain after return mapping will also decrease. Then we can get a smaller scale matrix 𝝋 for plasticity, i.e., smaller hourglass control for plasticity, which is what we want.

But, as you mentioned just now, the return mapping only works for shear stress, and theoretically, the shear strain should not change. I guess the stress-strain relationship is not reversible. We can get the stress from strain, but we can not get strain from stress, because one stress state may corresponds to multiple strains.

[DongWuTUM]

The shear stress is changing, but the J2 remains the same. I am not sure. Please check it. I think the stress-strain relationship is reversible as they are tensors.

[Shuaihao-Zhang]

The shear stress is changing, but the J2 remains the same. I am not sure. Please check it. I think the stress-strain relationship is reversible as they are tensors.

The second invariant J2 also changes after return mapping. If you check the return mapping process in I1-J2 (first Invariant and second Invariant ) space, it's clear. The following figure is the returning mapping in I1-J2 space for DP constitutive model [1].

The process B-C is the return mapping for deviatoric stress (I1 keeps unchanged). We can see the J2 (y-axis) also changed. And for J2 plasticity model, the return mapping process is similar with the DP model.

I think the stress-strain relationship is not reversible. We can get the stress from strain as one strain only corresponds to one stress, but we can not get strain from stress because on stress corresponds to several strains, as shown in the figure below.

[1]Bui, H. H., Fukagawa, R., Sako, K., Ohno, S., 2008. Lagrangian meshfree particles method (SPH) for large deformation and failure flows of geomaterial using elastic–plastic soil constitutive model. International journal for numerical and analytical methods in geomechanics 32(12), 1537-1570.

[DongWuTUM]

I have summarized the differences between the TL and UL.

1. Please check whether it is correct or not.
2. Check the difference in the scaling matrix when running cases of TL and UL.
3. Try: a. include the isotropic term before inverse operation (Please show me how you do it); b. Ɛ^(-1) = (Ɛ+I)^(-1) - I (Here, we can try Ɛ is shear strain or full strain); c. scaling matrix = (Ɛ_corrected+I)(Ɛ+I)^(-1).
4. The TL code is in branch "try_plastic_decomposed". (Note that the factor zeta in the code has not been adjusted. I will upload the newest code tomorrow.)
5. The shear stress-strain invertibility is not important here as the relation between them is only needed to times 2G.

[Xiangyu-Hu]

I have the feeling that the inconsistency is coming from the \phi term and \hat{v}_ij term. Since we are using full strain tensor for the second part of \hat{V}_ij, the \phi term should also be the same. The derivation is e_ij = \epsilon_ij^-1 \epsilon_ij e_ij -> \epsilon_ij^-1 (v_ij - \epsilon_ij e_ij) = \epsilon_ij^-1 * \hat{v}_ij.

[Xiangyu-Hu]

The suggested formulation on the correction term. Note that we need compute the effective strain rate within return mapping.

[Shuaihao-Zhang]

@Xiangyu-Hu @DongWuTUM Thanks a lot. I will check these in detail.

[Shuaihao-Zhang]

The suggested formulation on the correction term. Note that we need compute the effective strain rate within return mapping.

The new form is different from the previous formulation even for elastic case. I tested the new formulation and it doesn't work for elastic oscillating beam. It seems the correction term is too large. I do use B_correction here

[Xiangyu-Hu]

The suggested formulation on the correction term. Note that we need compute the effective strain rate within return mapping.

The new form is different from the previous formulation even for elastic case. I tested the new formulation and it doesn't work for elastic oscillating beam. It seems the correction term is too large. I do use B_correction here

It seems you misunderstood the new formulation, it should recover the elastic correction exactly before plasticity shows up. Because in that case, the two strain rates, yield identity.

[Shuaihao-Zhang]

It seems you misunderstood the new formulation, it should recover the elastic correction exactly before plasticity shows up. Because in that case, the two strain rates, yield identity.

I see. I guess there is a typo in you new formulation. The new scale matrix 𝝋 based on the strain rate in ok. But for the \hat{v}_ij, we should still use the old one. The new \hat{v}_ij used the strain rate. The strain rate incudes the transpose of velocity gradient, which we don't have in the old formulation.

Then, I tested the following formulation: The formulation can reproduce the elastic deformation, But for plasticity, it still crashes like this.

[Xiangyu-Hu]

It seems you misunderstood the new formulation, it should recover the elastic correction exactly before plasticity shows up. Because in that case, the two strain rates, yield identity.

I see. I guess there is a typo in you new formulation. The new scale matrix 𝝋 based on the strain rate in ok. But for the \hat{v}_ij, we should still use the old one. The new \hat{v}_ij used the strain rate. The strain rate incudes the transpose of velocity gradient, which we don't have in the old formulation.

Then, I tested the following formulation: The formulation can reproduce the elastic deformation, But for plasticity, it still crashes like this.

In the second case, how you obtain the effective strain rate? Also they should be strain rate not shear rate.

[DongWuTUM]

@Shuaihao-Zhang I have upload the newest houglass control code to my branch https://github.com/DongWuTUM/SPHinXsys/tree/plastic_hourglass_control. I don't have time to test in Linux, but it already works in Windows.

[Xiangyu-Hu]

[Xiangyu-Hu]

@Shuaihao-Zhang Could you update the code with this formulation? Thanks.

[Shuaihao-Zhang]

@Shuaihao-Zhang Could you update the code with this formulation? Thanks.

Sure, I have add the new formulation and the simulation still crashes quick after the bar touches the plate. You can check the code at: https://github.com/Xiangyu-Hu/SPHinXsys/blob/4352844bf976493eee43c7824be28fc2d2a8f7ce/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp#L551-L561 I used the difference of the shear strain after return mapping between time t and t-1 to calculate the effective shear strain rate. Then I also added the hydrostatic strain rate to get the effective total strain rate. I also tried with only the shear strain rate, and it also crashed.

[Shuaihao-Zhang]

@Shuaihao-Zhang I have upload the newest houglass control code to my branch https://github.com/DongWuTUM/SPHinXsys/tree/plastic_hourglass_control. I don't have time to test in Linux, but it already works in Windows.

Thanks! I will check it in detail.

[Xiangyu-Hu]

Please have a check on a simple formulation at the branch https://github.com/Xiangyu-Hu/SPHinXsys/tree/hourglass-control-ul_xy numerical stability seems OK.

[Shuaihao-Zhang]

Please have a check on a simple formulation at the branch https://github.com/Xiangyu-Hu/SPHinXsys/tree/hourglass-control-ul_xy numerical stability seems OK.

Thanks. Although the Taylor bar didn't crash with the new formulation, the scale coefficient is not an identity matrix in elastic situations. https://github.com/Xiangyu-Hu/SPHinXsys/blob/7dacc3d137e6328f1d23ad03e8060f9b65242b03/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp#L547 If I skip return mapping and use return shear_stress_rate_elastic;, ideally, the bar should bounce back like elastic material, but it didn't and the scale matrix is 0 everywhere.

Another thing is that I changed my branch to plane stress condition based on your branch for the J2 plasticity. Previously, I use the plane strain condition, which means I also need to use 3D stress tensor in 2D cases. So I use the reduceTensorand increaseTensor to transfer from 2D and 3D. But this just affect 2D cases.

I checked Dong's code, and he used plane stress assumption for TLSPH, which means that the stress is 2D for 2D cases and 3D for 3D cases. Since his result is good, I suppose the plans stress is also ok for the J2 plasticity (for soils, it's not ok). Then I just changed to plane stress condition based on your new branch. I deleted the transform between 2D and 3D tensors to simplify the code and reduce the source of error.

[Xiangyu-Hu]

Please have a check on a simple formulation at the branch https://github.com/Xiangyu-Hu/SPHinXsys/tree/hourglass-control-ul_xy numerical stability seems OK.

Thanks. Although the Taylor bar didn't crash with the new formulation, the scale coefficient is not an identity matrix in elastic situations.

https://github.com/Xiangyu-Hu/SPHinXsys/blob/7dacc3d137e6328f1d23ad03e8060f9b65242b03/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp#L547

If I skip return mapping and use return shear_stress_rate_elastic;, ideally, the bar should bounce back like elastic material, but it didn't and the scale matrix is 0 everywhere. Another thing is that I changed my branch to plane stress condition based on your branch for the J2 plasticity. Previously, I use the plane strain condition, which means I also need to use 3D stress tensor in 2D cases. So I use the reduceTensorand increaseTensor to transfer from 2D and 3D. But this just affect 2D cases.

I checked Dong's code, and he used plane stress assumption for TLSPH, which means that the stress is 2D for 2D cases and 3D for 3D cases. Since his result is good, I suppose the plans stress is also ok for the J2 plasticity (for soils, it's not ok). Then I just changed to plane stress condition based on your new branch. I deleted the transform between 2D and 3D tensors to simplify the code and reduce the source of error.