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

Shuaihao-Zhang · 2023-11-24T15:58:26Z

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

Shuaihao-Zhang · 2023-11-24T16:38:39Z

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.

SPHinXsys/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp

Lines 545 to 553 in 06d36e8

    
           // calculate scale matrix 
        
           if (plastic_indicator_[index_i]) 
        
           { 
        
               scale_coef_[index_i] = reduceTensor(shear_stress_3D_[index_i]) * reduceTensor(shear_strain_3D_[index_i]).inverse() / (2.0 * G_); 
        
           } 
        
           else 
        
           { 
        
               scale_coef_[index_i] = Matd::Identity(); 
        
           }

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 · 2023-11-24T17:47:38Z

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

Xiangyu-Hu · 2023-11-24T17:48:47Z

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

Xiangyu-Hu · 2023-11-24T17:53:36Z

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.

SPHinXsys/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp

Lines 545 to 553 in 06d36e8

// calculate scale matrix

if (plastic_indicator_[index_i])

{

scale_coef_[index_i] = reduceTensor(shear_stress_3D_[index_i]) * reduceTensor(shear_strain_3D_[index_i]).inverse() / (2.0 * G_);

}

else

{

scale_coef_[index_i] = Matd::Identity();

}

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 · 2023-11-25T02:37:08Z

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 · 2023-11-25T02:45:29Z

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 · 2023-11-25T02:49:39Z

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 · 2023-11-26T16:08:23Z

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.

Xiangyu-Hu · 2023-11-26T16:22:03Z

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 · 2023-11-26T16:41:59Z

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.

SPHinXsys/tests/user_examples/extra_src/shared/general_continuum/general_continuum.cpp

Lines 92 to 114 in 2eac06a

    
           Mat3d J2Plasticity::ConstitutiveRelationShearStress(Mat3d& velocity_gradient, Mat3d& shear_stress) 
        
           { 
        
               Real dim = 3.0; 
        
               Mat3d strain_rate = 0.5 * (velocity_gradient + velocity_gradient.transpose()); 
        
               Mat3d spin_rate = 0.5 * (velocity_gradient - velocity_gradient.transpose()); 
        
               Mat3d deviatoric_strain_rate = strain_rate - (1.0 / dim) * strain_rate.trace() * Mat3d::Identity(); 
        
               //consider the elastic part 
        
               Mat3d shear_stress_rate_elastic = 2.0 * G_ * deviatoric_strain_rate; 
        
               //consider the plastic part 
        
               Mat3d deviatoric_stress_tensor = shear_stress; 
        
               Real stress_tensor_J2 = 0.5 * (deviatoric_stress_tensor.cwiseProduct(deviatoric_stress_tensor.transpose())).sum(); 
        
               Real f = sqrt(2.0 * stress_tensor_J2) - sqrt(2.0 / 3.0) * yield_stress_; 
        
               Real lambda_dot_ = 0; 
        
               Mat3d g = Mat3d::Zero(); 
        
               if (f > TinyReal) 
        
               { 
        
                   Real deviatoric_stress_times_strain_rate = (deviatoric_stress_tensor.cwiseProduct(strain_rate)).sum(); 
        
                   lambda_dot_ = deviatoric_stress_times_strain_rate / sqrt(2.0 * stress_tensor_J2); 
        
                   g = lambda_dot_ * (sqrt(2.0) * G_ * deviatoric_stress_tensor / (sqrt(stress_tensor_J2))); 
        
               } 
        
               return shear_stress_rate_elastic - g; 
        
               //return shear_stress_rate_elastic; 
        
           }

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 · 2023-11-26T16:42:52Z

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 · 2023-11-26T17:28:53Z

@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 · 2023-11-26T17:31:41Z

@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 · 2023-11-26T17:38:13Z

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.

SPHinXsys/tests/user_examples/extra_src/shared/general_continuum/general_continuum.cpp

Lines 92 to 114 in 2eac06a

Mat3d J2Plasticity::ConstitutiveRelationShearStress(Mat3d& velocity_gradient, Mat3d& shear_stress)

{

Real dim = 3.0;

Mat3d strain_rate = 0.5 * (velocity_gradient + velocity_gradient.transpose());

Mat3d spin_rate = 0.5 * (velocity_gradient - velocity_gradient.transpose());

Mat3d deviatoric_strain_rate = strain_rate - (1.0 / dim) * strain_rate.trace() * Mat3d::Identity();

//consider the elastic part

Mat3d shear_stress_rate_elastic = 2.0 * G_ * deviatoric_strain_rate;

//consider the plastic part

Mat3d deviatoric_stress_tensor = shear_stress;

Real stress_tensor_J2 = 0.5 * (deviatoric_stress_tensor.cwiseProduct(deviatoric_stress_tensor.transpose())).sum();

Real f = sqrt(2.0 * stress_tensor_J2) - sqrt(2.0 / 3.0) * yield_stress_;

Real lambda_dot_ = 0;

Mat3d g = Mat3d::Zero();

if (f > TinyReal)

{

Real deviatoric_stress_times_strain_rate = (deviatoric_stress_tensor.cwiseProduct(strain_rate)).sum();

lambda_dot_ = deviatoric_stress_times_strain_rate / sqrt(2.0 * stress_tensor_J2);

g = lambda_dot_ * (sqrt(2.0) * G_ * deviatoric_stress_tensor / (sqrt(stress_tensor_J2)));

}

return shear_stress_rate_elastic - g;

//return shear_stress_rate_elastic;

}

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.

SPHinXsys/tests/user_examples/extra_src/shared/general_continuum/general_continuum.cpp

Lines 92 to 114 in 2eac06a

Mat3d J2Plasticity::ConstitutiveRelationShearStress(Mat3d& velocity_gradient, Mat3d& shear_stress)

{

Real dim = 3.0;

Mat3d strain_rate = 0.5 * (velocity_gradient + velocity_gradient.transpose());

Mat3d spin_rate = 0.5 * (velocity_gradient - velocity_gradient.transpose());

Mat3d deviatoric_strain_rate = strain_rate - (1.0 / dim) * strain_rate.trace() * Mat3d::Identity();

//consider the elastic part

Mat3d shear_stress_rate_elastic = 2.0 * G_ * deviatoric_strain_rate;

//consider the plastic part

Mat3d deviatoric_stress_tensor = shear_stress;

Real stress_tensor_J2 = 0.5 * (deviatoric_stress_tensor.cwiseProduct(deviatoric_stress_tensor.transpose())).sum();

Real f = sqrt(2.0 * stress_tensor_J2) - sqrt(2.0 / 3.0) * yield_stress_;

Real lambda_dot_ = 0;

Mat3d g = Mat3d::Zero();

if (f > TinyReal)

{

Real deviatoric_stress_times_strain_rate = (deviatoric_stress_tensor.cwiseProduct(strain_rate)).sum();

lambda_dot_ = deviatoric_stress_times_strain_rate / sqrt(2.0 * stress_tensor_J2);

g = lambda_dot_ * (sqrt(2.0) * G_ * deviatoric_stress_tensor / (sqrt(stress_tensor_J2)));

}

return shear_stress_rate_elastic - g;

//return shear_stress_rate_elastic;

}

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_stress_3D_[index_i];
The simulation now should not blowup. but it does. So, there seems still has consistent issue.

Xiangyu-Hu · 2023-11-26T17:42:27Z

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

Shuaihao-Zhang · 2023-11-27T01:36:51Z

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 · 2023-11-27T01:44:39Z

@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 · 2023-11-27T02:47:07Z

@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 · 2023-11-27T06:30:34Z

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.

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 · 2023-11-27T06:34:48Z

@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 · 2023-11-27T06:37:04Z

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

Shuaihao-Zhang · 2023-11-27T07:22:42Z

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 · 2023-11-27T09:09:38Z

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

Xiangyu-Hu · 2023-11-27T09:12:35Z

We better begin with the simplest plastic case first.

Shuaihao-Zhang · 2023-11-27T09:15:51Z

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.

DongWuTUM · 2023-11-27T11:34:42Z

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.

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 · 2023-11-27T12:34:38Z

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.

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 · 2023-11-27T13:40:39Z

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 · 2023-11-27T13:55:24Z

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 · 2023-11-27T16:24:50Z

I have summarized the differences between the TL and UL.

Please check whether it is correct or not.
Check the difference in the scaling matrix when running cases of TL and UL.
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).
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.)
The shear stress-strain invertibility is not important here as the relation between them is only needed to times 2G.

Xiangyu-Hu · 2023-11-27T17:21:08Z

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 · 2023-11-28T08:51:22Z

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

Shuaihao-Zhang · 2023-11-28T12:04:31Z

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

Shuaihao-Zhang · 2023-11-29T03:59:50Z

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 · 2023-11-29T06:47:33Z

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 · 2023-11-29T09:03:28Z

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 · 2023-11-29T09:31:12Z

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 · 2023-11-29T09:43:43Z

@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 · 2023-11-29T11:11:31Z

Xiangyu-Hu · 2023-11-29T11:25:55Z

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

Shuaihao-Zhang · 2023-11-29T17:02:14Z

@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:

SPHinXsys/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp

Lines 551 to 561 in 4352844

    
           /*shear strain and shear strain rate befor and after return mapping*/ 
        
           shear_strain_[index_i] = reduceTensor(shear_strain_3D_[index_i]); 
        
           shear_strain_return_map_[index_i] = reduceTensor(shear_stress_3D_[index_i]) / (2.0 * G_); 
        
           //calculate the shear_strain_rate_return_map_ based on the strain increment between time step n and n-1 
        
           Real hydrostatic_strain_rate_i = strain_rate_3D_[index_i].trace() / 3.0; 
        
           shear_strain_rate_return_map_[index_i] = (shear_strain_return_map_[index_i] - shear_strain_pre_return_map_[index_i]) / dt; 
        
           //calculate total strain rate 
        
           Matd total_strain_rate_i = reduceTensor(strain_rate_3D_[index_i]); 
        
           Matd total_strain_rate_i_return_map = shear_strain_rate_return_map_[index_i] + hydrostatic_strain_rate_i * Matd::Identity(); 
        
           //shear strain in the time step n-1 
        
           shear_strain_pre_return_map_[index_i] = shear_strain_return_map_[index_i];

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 · 2023-11-29T17:16:42Z

@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 · 2023-11-29T23:57:53Z

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 · 2023-11-30T07:22:18Z

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.

SPHinXsys/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp

Line 547 in 7dacc3d

    
           scale_coef_[index_i] = reduceTensor(shear_stress_3D_[index_i]) * (velocity_gradient_[index_i] + Eps * Matd::Identity()).inverse();

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 reduceTensor and 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 · 2023-11-30T08:34:20Z

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.

SPHinXsys/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp

Line 547 in 7dacc3d

scale_coef_[index_i] = reduceTensor(shear_stress_3D_[index_i]) * (velocity_gradient_[index_i] + Eps * Matd::Identity()).inverse();

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 reduceTensor and 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 · 2023-11-30T09:16:23Z

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.

SPHinXsys/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp

Line 547 in 7dacc3d

scale_coef_[index_i] = reduceTensor(shear_stress_3D_[index_i]) * (velocity_gradient_[index_i] + Eps * Matd::Identity()).inverse();

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 reduceTensor and 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.

It works for elastic (both elastic shear stress rate and no returnning map)
3D case as shown in the test branch https://github.com/Xiangyu-Hu/SPHinXsys/tree/hourglass-control-ul_xy
Note that the correction is numerical force term, there is no requirement for the consistency mentioned by you. More important is they approaching zero when resolution is high and it is small value compared to the main stress induced force.
I made small change on the scaling factor, however i do not known which is better.

DongWuTUM · 2023-11-30T09:40:42Z

@Shuaihao-Zhang It is the plane strain (not plane stress) problem for 2D case. The stress in third coordinate is not needed for calculation. (I remember the stress in third coordinate is needed for soil. You can ignore this stress now.)

Shuaihao-Zhang · 2023-12-01T01:19:08Z

It works for elastic (both elastic shear stress rate and no returnning map) 3D case as shown in the test branch https://github.com/Xiangyu-Hu/SPHinXsys/tree/hourglass-control-ul_xy Note that the correction is numerical force term, there is no requirement for the consistency mentioned by you. More important is they approaching zero when resolution is high and it is small value compared to the main stress induced force. I made small change on the scaling factor, however i do not known which is better.

Thanks. I checked you branch. The 3D Taylor bar seems elastic if we use elastic shear stress rate and no return mapping.
But the scale matrix if super large, and theoretically, it should be an identity matrix for elasticity.

Also for 2D Taylor bar, the scale matrix is very large although is didn't crash.

For 2D oscillating beam, if I use the plastic constitutive relation (with elastic shear stress rate and no return mapping), the result should be exactly the same with the elastic beam (for previous method, we can total reproduce this). But it crashes and I think this is related to the unreasonable scale matrix.

Shuaihao-Zhang · 2023-12-01T01:19:45Z

@Shuaihao-Zhang It is the plane strain (not plane stress) problem for 2D case. The stress in third coordinate is not needed for calculation. (I remember the stress in third coordinate is needed for soil. You can ignore this stress now.)

Thanks! Then I will just use 2D matrix for 2D cases.

Xiangyu-Hu · 2023-12-01T23:14:13Z

@Shuaihao-Zhang @DongWuTUM
Finally, matrix inverse regularization is done. Now oscillation beam and Taylor bar both work.
shear_stress_try

Shuaihao-Zhang · 2023-12-02T04:52:47Z

@Shuaihao-Zhang @DongWuTUM Finally, matrix inverse regularization is done. Now oscillation beam and Taylor bar both work. shear_stress_try

Thanks! Now the elastic deformation can be totally reproduced.

draft PR -- hourglass control for elasticity and plasticity

06d36e8

Shuaihao-Zhang requested review from Xiangyu-Hu and DongWuTUM November 24, 2023 16:39

Shuaihao-Zhang added 2 commits November 27, 2023 00:27

remove unnecessary code blocks and comments

9dc4811

remove unnecessary code blocks and comments

2eac06a

new formulation with strain rate

4352844

change to plane stress condition

58782c3

Xiangyu-Hu closed this Nov 30, 2023

Xiangyu-Hu reopened this Nov 30, 2023

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

Are you sure you want to change the base?

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

Conversation

Shuaihao-Zhang commented Nov 24, 2023 • edited

Shuaihao-Zhang commented Nov 24, 2023

Xiangyu-Hu commented Nov 24, 2023

Xiangyu-Hu commented Nov 24, 2023

Xiangyu-Hu commented Nov 24, 2023

Shuaihao-Zhang commented Nov 25, 2023

Shuaihao-Zhang commented Nov 25, 2023

Shuaihao-Zhang commented Nov 25, 2023

Xiangyu-Hu commented Nov 26, 2023 • edited

Xiangyu-Hu commented Nov 26, 2023

Shuaihao-Zhang commented Nov 26, 2023

Shuaihao-Zhang commented Nov 26, 2023

DongWuTUM commented Nov 26, 2023

DongWuTUM commented Nov 26, 2023 • edited

Xiangyu-Hu commented Nov 26, 2023

Xiangyu-Hu commented Nov 26, 2023 • edited

Shuaihao-Zhang commented Nov 27, 2023

Shuaihao-Zhang commented Nov 27, 2023

Shuaihao-Zhang commented Nov 27, 2023

Xiangyu-Hu commented Nov 27, 2023

Xiangyu-Hu commented Nov 27, 2023

Xiangyu-Hu commented Nov 27, 2023

Shuaihao-Zhang commented Nov 27, 2023

Xiangyu-Hu commented Nov 27, 2023

Xiangyu-Hu commented Nov 27, 2023

Shuaihao-Zhang commented Nov 27, 2023

DongWuTUM commented Nov 27, 2023

Shuaihao-Zhang commented Nov 27, 2023

DongWuTUM commented Nov 27, 2023

Shuaihao-Zhang commented Nov 27, 2023

DongWuTUM commented Nov 27, 2023

I have summarized the differences between the TL and UL.

Xiangyu-Hu commented Nov 27, 2023 • edited

Xiangyu-Hu commented Nov 28, 2023

Shuaihao-Zhang commented Nov 28, 2023

Shuaihao-Zhang commented Nov 29, 2023

Xiangyu-Hu commented Nov 29, 2023 • edited

Shuaihao-Zhang commented Nov 29, 2023

Xiangyu-Hu commented Nov 29, 2023 • edited

DongWuTUM commented Nov 29, 2023

Xiangyu-Hu commented Nov 29, 2023

Xiangyu-Hu commented Nov 29, 2023

Shuaihao-Zhang commented Nov 29, 2023

Shuaihao-Zhang commented Nov 29, 2023

Xiangyu-Hu commented Nov 29, 2023

Shuaihao-Zhang commented Nov 30, 2023

Xiangyu-Hu commented Nov 30, 2023

Xiangyu-Hu commented Nov 30, 2023

DongWuTUM commented Nov 30, 2023

Shuaihao-Zhang commented Dec 1, 2023

Shuaihao-Zhang commented Dec 1, 2023

Xiangyu-Hu commented Dec 1, 2023

Shuaihao-Zhang commented Dec 2, 2023

Shuaihao-Zhang commented Nov 24, 2023 •

edited

Xiangyu-Hu commented Nov 26, 2023 •

edited

DongWuTUM commented Nov 26, 2023 •

edited

Xiangyu-Hu commented Nov 26, 2023 •

edited

Xiangyu-Hu commented Nov 27, 2023 •

edited

Xiangyu-Hu commented Nov 29, 2023 •

edited

Xiangyu-Hu commented Nov 29, 2023 •

edited