SPHinXsys
SPHinXsys copied to clipboard
Xiangyu-Hu
draft PR -- hourglass control (updated Lagrangian) for elasticity and plasticity
Here, I will show the results when dealing with the hourglass control method for ULSPH.
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 i is 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.
This is the method we discussed just now. We only calculate the inverse of shear strain for
particle iwhen the return mapping is called, i.e., whenparticle iis plastic. Then, we can calculate the scale matrix. Ifparticle 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
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}?
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.
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.
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.
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 = J2_plasticity_.ReturnMappingShearStress(shear_stress_3D_[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.
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 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 = J2_plasticity_.ReturnMappingShearStress(shear_stress_3D_[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.
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.
@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?
@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.
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 = J2_plasticity_.ReturnMappingShearStress(shear_stress_3D_[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 = J2_plasticity_.ReturnMappingShearStress(shear_stress_3D_[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.
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_stress_3D_[index_i]; The simulation now should not blowup. but it does. So, there seems still has consistent issue.
Because it does not blowup if we simply constrain scale factor to identity in any case.
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_stress_3D_[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, 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.
@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.
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_stress_3D_[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 533shear_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 @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.
It is clear that if we included the isotropic part the variation of the scaling factor will be smaller.
It is clear that if we included the isotropic part the variation of the scaling factor will be smaller.
Sure, I will check this.
Another thing is that the transformation between 2 and 3 d may introduce consistent issue.
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.
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.
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.
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.
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_stress_3D_[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 533shear_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.
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.