Material API and Shell kinematics

[FabienPean-Virtonomy]

The issue I want to raise is twofold.

First, the balance equation for shell implicitly relies on material parameters to set some part of the Almansi strain measure

SPHinXsys/src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.cpp

Line 169 in 9b947a1

current_local_almansi_strain = getCorrectedAlmansiStrain(current_local_almansi_strain, nu_);

This is a problem for nonlinear or anisotropic materials where the notion does not exist per say (or on a per-step basis). The formulation is obtained from a linear elastic assumption according to the manuscript on arxiv. Besides the plane stress condition is enforced after. The strain tensor is only one of the measure that may be used for a material formulation.

We often end up to the point where we need $\mathbf{B}$ which clutter the material with

SPH::Mat3d B = (I-2*e).inverse();

Differentiating for Cauchy stress and 2nd PK stress creates an asymmetry in the interface which is confusing during development, in simpler cases one should be able to compute on as function of the other.

The second point is about the frame into which those values are expressed. The Almansi strain is given in current shell-local coordinates. However anisotropic materials rely on one or more preferred directions, often expressed in reference coordinates. With shell working in shell-local coordinates, it means those directions must be transformed at some point of the stress calculation pipeline.

SPHinXsys/src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.cpp

Lines 161 to 181 in 9b947a1

	Matd F_gaussian_point = F_[index_i] + gaussian_point_[i] * F_bending_[index_i] * thickness_[index_i] * 0.5;
	Matd dF_gaussian_point_dt = dF_dt_[index_i] + gaussian_point_[i] * dF_bending_dt_[index_i] * thickness_[index_i] * 0.5;
	Matd inverse_F_gaussian_point = F_gaussian_point.inverse();
	Matd current_local_almansi_strain = current_transformation_matrix * transformation_matrix0_[index_i].transpose() * 0.5 *
	(Matd::Identity() - inverse_F_gaussian_point.transpose() * inverse_F_gaussian_point) *
	transformation_matrix0_[index_i] * current_transformation_matrix.transpose();
	/** correct Almansi strain tensor according to plane stress problem. */
	current_local_almansi_strain = getCorrectedAlmansiStrain(current_local_almansi_strain, nu_);
	/** correct out-plane numerical damping. */
	numerical_damping_scaling_matrix_(Dimensions - 1, Dimensions - 1) = thickness_[i] < smoothing_length_ ? thickness_[i] : smoothing_length_;
	Matd cauchy_stress = elastic_solid_.StressCauchy(current_local_almansi_strain, index_i) +
	current_transformation_matrix * transformation_matrix0_[index_i].transpose() * F_gaussian_point *
	elastic_solid_.NumericalDampingRightCauchy(F_gaussian_point, dF_gaussian_point_dt, numerical_damping_scaling_matrix_, index_i) *
	F_gaussian_point.transpose() * transformation_matrix0_[index_i] * current_transformation_matrix.transpose() / F_gaussian_point.determinant();
	/** Impose modeling assumptions. */
	cauchy_stress.col(Dimensions - 1) *= shear_correction_factor_;
	cauchy_stress.row(Dimensions - 1) *= shear_correction_factor_;
	cauchy_stress(Dimensions - 1, Dimensions - 1) = 0.0;

Currently SPHinXsys dodges calculation of Cauchy stress for complex materials (see elastic_solid.cpp)

Standard approach and my recommendation is to provide the deformation gradient in global coordinates as input to the material (and return stress in global coordinates). It allows a uniform material formulation independent of the kinematics (beam,shell,volume) and would allow to reuse existing implementation of material (because then $\sigma=\mathbf{FSF}^T/J$ is true)

virtual SPH::Mat3d StressPK2(SPH::Mat3d const& F, size_t i) = 0;
virtual SPH::Mat3d StressCauchy(SPH::Mat3d const& F, size_t i) {return F * StressPK2(F,i) * F.transpose() / F.determinant();};

[DongWuTUM]

Hi Fabien,

We can simply use the linear elastic assumption for plane stress conditions. I think it will not introduce much error.

We can also use other constitutive relations. After that, the calculated stress can be transformed into Cauchy stress. The deformation gradient tensor F is defined in initial local coordinates, which is not relevant to the applied constitutive relations.
I mean if we input the local F into constitutive relations, then we get local stress. We have to get the local F, as you can see from the code or paper, the last column of F is related to the local normal direction.

We can have further discussion.

[FabienPean-Virtonomy]

Following-up the discussion at the conference, here is what I could find on shell theory and reusing 3D material formulation that I believe relevant:

• Simo, J.C., Rifai, M.S., Fox, D.D., 1990. On a stress resultant geometrically exact shell model. Part IV: Variable thickness shells with through-the-thickness stretching. Comput. Methods Appl. Mech. Eng. 81, 91–126. https://doi.org/10.1016/0045-7825(90)90143-A
• Betsch, P., Gruttmann, F., Stein, E., 1996. A 4-node finite shell element for the implementation of general hyperelastic 3D-elasticity at finite strains. Comput. Methods Appl. Mech. Eng. 130, 57–79. https://doi.org/10.1016/0045-7825(95)00920-5
• Kiendl, J., Hsu, M.C., Wu, M.C.H., Reali, A., 2015. Isogeometric Kirchhoff-Love shell formulations for general hyperelastic materials. Comput. Methods Appl. Mech. Eng. 291, 280–303. https://doi.org/10.1016/j.cma.2015.03.010
• Klinkel, S., Govindjee, S., 2002. Using finite strain 3D-material models in beam and shell elements. Eng. Comput. 19, 902–921. https://doi.org/10.1108/02644400210450341
• FEBio Theory Manual https://help.febio.org/docs/FEBioTheory-4-5/TM45-Section-4.2.html

[DongWuTUM]

Following-up the discussion at the conference, here is what I could find on shell theory and reusing 3D material formulation that I believe relevant:

• Simo, J.C., Rifai, M.S., Fox, D.D., 1990. On a stress resultant geometrically exact shell model. Part IV: Variable thickness shells with through-the-thickness stretching. Comput. Methods Appl. Mech. Eng. 81, 91–126. https://doi.org/10.1016/0045-7825(90)90143-A
• Betsch, P., Gruttmann, F., Stein, E., 1996. A 4-node finite shell element for the implementation of general hyperelastic 3D-elasticity at finite strains. Comput. Methods Appl. Mech. Eng. 130, 57–79. https://doi.org/10.1016/0045-7825(95)00920-5
• Kiendl, J., Hsu, M.C., Wu, M.C.H., Reali, A., 2015. Isogeometric Kirchhoff-Love shell formulations for general hyperelastic materials. Comput. Methods Appl. Mech. Eng. 291, 280–303. https://doi.org/10.1016/j.cma.2015.03.010
• Klinkel, S., Govindjee, S., 2002. Using finite strain 3D-material models in beam and shell elements. Eng. Comput. 19, 902–921. https://doi.org/10.1108/02644400210450341

Got it! Thanks.