Fix Saint-Venant Kirchhoff implementation

[FabienPean-Virtonomy]

StressCauchy is inherited from LinearElasticSolid which is an incorrect formula.

This PR fixes that and restores $\mathbf{F}\mathbf{S}\mathbf{F}^T/J = \boldsymbol\sigma$ invariant

Implementation follows https://help.febio.org/Manuals/FEBioTheory/FEBio_tm_3-0-Subsection-5.2.1.html

[FabienPean-Virtonomy]

The CI output and the issue ought to be discussed in the weekly meeting,

The following tests FAILED:

1205   52 - test_2d_oscillating_beam_cauchy (Failed)
1206   55 - test_2d_plate (Failed)
1207   60 - test_2d_shell (Subprocess aborted)
1208  122 - test_3d_dambreak_elastic_plate_shell (Failed)
1209  156 - sphere_compression.half_sphere (Failed)
1210  157 - test_3d_shell_stability_half_sphere (Failed)
1211  158 - test_3d_slender_beam (Failed)
1212  167 - test_3d_thin_plate (Failed)

[Xiangyu-Hu]

@DongWuTUM Could you also have check? Thanks.

[DongWuTUM]

@FabienPean-Virtonomy Hi Fabien, I totally agree with you. The expression $S = \lambda (tr \mathbf{E}) \mathbf{I} + 2 \mu \mathbf{E}$ cannot be directly converted to $\sigma = \lambda (tr \mathbf{\epsilon}) \mathbf{I} + 2 \mu \mathbf{\epsilon}$ using the transformation $E = F^T \epsilon F$. Both constitutive relations are indeed valid, but the key question is whether we aim to obtain an "invariant" form, as you've shown, or if we acknowledge the difference and maintain both forms for specific applications.

Have you checked how much the results differ when using these two constitutive relations? I think the difference is small.

[WeiyiVirtonomy]

The CI output and the issue ought to be discussed in the weekly meeting,

The following tests FAILED:
1205   52 - test_2d_oscillating_beam_cauchy (Failed)
1206   55 - test_2d_plate (Failed)
1207   60 - test_2d_shell (Subprocess aborted)
1208  122 - test_3d_dambreak_elastic_plate_shell (Failed)
1209  156 - sphere_compression.half_sphere (Failed)
1210  157 - test_3d_shell_stability_half_sphere (Failed)
1211  158 - test_3d_slender_beam (Failed)
1212  167 - test_3d_thin_plate (Failed)

It seems that all the failed tests are using the shell model. Are the failures all caused by the negative B matrix computed from the corrected Almansi strain, which is based on the linear elastic material assumption?

[FabienPean-Virtonomy]

Both constitutive relations are indeed valid, but the key question is whether we aim to obtain an "invariant" form, as you've shown, or if we acknowledge the difference and maintain both forms for specific applications.

Have you checked how much the results differ when using these two constitutive relations? I think the difference is small.

The equality must obviously be preserved, a material model has to behave exactly as per its definition no matter the stress measure one asks for. Personally, I would shudder at a world where $1+1=-3$ 😄

To the best of my knowledge, the equation implemented for the Cauchy stress in LinearElasticSolid which represents the linear stress-strain relation using the Cauchy-Almansi conjugate pair has no name and is not used in literature.

The equivalence $\mathbf{E}\approx\mathbf{e}\approx\boldsymbol\varepsilon$ is acceptable only at infinitesimal strain. The difference is tremendous at finite strain, indeed the Cauchy stress for the StVK model evolves as $\lambda^4_i$, while the current formula evolves as $\lambda^{-2}_i$, where $\lambda_i$ are the principal stretches.

@WeiyiVirtonomy Not all, some tests seem to finish without a NaN

[DongWuTUM]

Both constitutive relations are indeed valid, but the key question is whether we aim to obtain an "invariant" form, as you've shown, or if we acknowledge the difference and maintain both forms for specific applications.
Have you checked how much the results differ when using these two constitutive relations? I think the difference is small.

The equality must obviously be preserved, a material model has to behave exactly as per its definition no matter the stress measure one asks for. Personally, I would shudder at a world where 1 + 1 = − 3 😄

To the best of my knowledge, the equation implemented for the Cauchy stress in LinearElasticSolid which represents the linear stress-strain relation using the Cauchy-Almansi conjugate pair has no name and is not used in literature.

The equivalence E ≈ e ≈ ε is acceptable only at infinitesimal strain. The difference is tremendous at finite strain, indeed the Cauchy stress for the StVK model evolves as λ i 4 , while the current formula evolves as λ i − 2 , where λ i are the principal stretches.

@WeiyiVirtonomy Not all, some tests seem to finish without a NaN

$\sigma = \lambda (tr \mathbf{\epsilon}) \mathbf{I} + 2 \mu \mathbf{\epsilon}$ is commonly used. Check the same website https://help.febio.org/FEBio/FEBio_tm_2_7/FEBio_tm_2-7-Section-5.1.html. Indeed, the difference is small only when there is small deformation because of the linear assumption.

[FabienPean-Virtonomy]

σ = λ ( t r ϵ ) I + 2 μ ϵ is commonly used. Check the same website help.febio.org/FEBio/FEBio_tm_2_7/FEBio_tm_2-7-Section-5.1.html. Indeed, the difference is small only when there is small deformation because of the linear assumption.

The small strain formula is commonly used, as used in StressPK2 function of class LinearElasticSolid, not the Cauchy-Almansi as in StressCauchy. We both agree that the Cauchy-Almansi is a valid implementation for small strain elasticity yes, not at finite strain for any other material model.

SPHinXsys/src/shared/materials/elastic_solid.cpp

Lines 81 to 90 in 00815cb

	Matd LinearElasticSolid::StressPK2(Matd &F, size_t index_i)
	{
	Matd strain = 0.5 * (F.transpose() + F) - Matd::Identity();
	return lambda0_ * strain.trace() * Matd::Identity() + 2.0 * G0_ * strain;
	}
	//=================================================================================================//
	Matd LinearElasticSolid::StressCauchy(Matd &almansi_strain, size_t index_i)
	{
	return lambda0_ * almansi_strain.trace() * Matd::Identity() + 2.0 * G0_ * almansi_strain;
	}