Skip to main content

conspire/constitutive/hybrid/elastic/multiplicative/elastic/
mod.rs

1#[cfg(test)]
2mod test;
3use crate::math::Quantity;
4use crate::math::{Current, Intermediate, Reference};
5
6use crate::{
7    constitutive::{
8        ConstitutiveError,
9        hybrid::ElasticMultiplicative,
10        solid::{
11            Solid,
12            elastic::{Elastic, internal_variables::ElasticIV},
13        },
14    },
15    math::{ContractThirdWithFirst, Rank2, TensorArray, TensorRank4},
16    mechanics::{
17        CauchyStress, CauchyTangentStiffness, CauchyTangentStiffness1, DeformationGradient,
18        DeformationGradient2, FirstPiolaKirchhoffStress, FirstPiolaKirchhoffStress1,
19        FirstPiolaKirchhoffStress2, FirstPiolaKirchhoffTangentStiffness,
20        FirstPiolaKirchhoffTangentStiffness2, SecondPiolaKirchhoffStress,
21    },
22    units::Stress,
23};
24
25impl<C1, C2> Solid for ElasticMultiplicative<C1, C2>
26where
27    C1: Elastic,
28    C2: Elastic,
29{
30    fn bulk_modulus(&self) -> Quantity<Stress> {
31        1.0 / (1.0 / self.0.bulk_modulus() + 1.0 / self.1.bulk_modulus())
32    }
33    fn shear_modulus(&self) -> Quantity<Stress> {
34        1.0 / (1.0 / self.0.shear_modulus() + 1.0 / self.1.shear_modulus())
35    }
36}
37
38impl<C1, C2> ElasticIV<DeformationGradient2> for ElasticMultiplicative<C1, C2>
39where
40    C1: Elastic,
41    C2: Elastic,
42{
43    type Residual = FirstPiolaKirchhoffStress2;
44    type TangentVu = TensorRank4<3, Intermediate, Reference, Current, Reference, Stress>;
45    type TangentUv = TensorRank4<3, Current, Reference, Intermediate, Reference, Stress>;
46    type TangentVv = FirstPiolaKirchhoffTangentStiffness2;
47    /// Calculates and returns the Cauchy stress.
48    ///
49    /// ```math
50    /// \boldsymbol{\sigma} = \frac{1}{J_2}\,\boldsymbol{\sigma}_1
51    /// ```
52    fn cauchy_stress(
53        &self,
54        deformation_gradient: &DeformationGradient,
55        deformation_gradient_2: &DeformationGradient2,
56    ) -> Result<CauchyStress, ConstitutiveError> {
57        let (deformation_gradient_2_inverse, jacobian_2) =
58            deformation_gradient_2.inverse_and_determinant();
59        let deformation_gradient_1 = deformation_gradient * &deformation_gradient_2_inverse;
60        Ok(self.0.cauchy_stress(&deformation_gradient_1.into())? / jacobian_2)
61    }
62    /// Calculates and returns the tangent stiffness associated with the Cauchy stress.
63    ///
64    /// ```math
65    /// \boldsymbol{\mathcal{T}} = \frac{1}{J_2}\,\boldsymbol{\mathcal{T}}_1\cdot\mathbf{F}_2^{-T}
66    /// ```
67    fn cauchy_tangent_stiffness(
68        &self,
69        deformation_gradient: &DeformationGradient,
70        deformation_gradient_2: &DeformationGradient2,
71    ) -> Result<CauchyTangentStiffness, ConstitutiveError> {
72        let (deformation_gradient_2_inverse, jacobian_2) =
73            deformation_gradient_2.inverse_and_determinant();
74        let deformation_gradient_1 = deformation_gradient * &deformation_gradient_2_inverse;
75        Ok(CauchyTangentStiffness1::from(
76            self.0
77                .cauchy_tangent_stiffness(&deformation_gradient_1.into())?,
78        ) * (deformation_gradient_2_inverse.transpose() / jacobian_2))
79    }
80    /// Calculates and returns the first Piola-Kirchhoff stress.
81    ///
82    /// ```math
83    /// \mathbf{P} = \mathbf{P}_1\cdot\mathbf{F}_2^{-T}
84    /// ```
85    fn first_piola_kirchhoff_stress(
86        &self,
87        deformation_gradient: &DeformationGradient,
88        deformation_gradient_2: &DeformationGradient2,
89    ) -> Result<FirstPiolaKirchhoffStress, ConstitutiveError> {
90        Ok(
91            self.cauchy_stress(deformation_gradient, deformation_gradient_2)?
92                * deformation_gradient.inverse_transpose()
93                * deformation_gradient.determinant(),
94        )
95    }
96    /// Calculates and returns the second Piola-Kirchhoff stress.
97    ///
98    /// ```math
99    /// \mathbf{S} = \mathbf{F}_2^{-1}\cdot\mathbf{S}_1\cdot\mathbf{F}_2^{-T}
100    /// ```
101    fn second_piola_kirchhoff_stress(
102        &self,
103        deformation_gradient: &DeformationGradient,
104        deformation_gradient_2: &DeformationGradient2,
105    ) -> Result<SecondPiolaKirchhoffStress, ConstitutiveError> {
106        Ok(deformation_gradient.inverse()
107            * self.first_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_2)?)
108    }
109    fn internal_variables_initial(&self) -> DeformationGradient2 {
110        DeformationGradient2::identity()
111    }
112    /// Calculates and returns the residual associated with the second deformation gradient.
113    ///
114    /// ```math
115    /// \mathbf{R} = \mathbf{P}_2 - \mathbf{M}_1\cdot\mathbf{F}_2^{-T}
116    /// ```
117    fn internal_variables_residual(
118        &self,
119        deformation_gradient: &DeformationGradient,
120        deformation_gradient_2: &DeformationGradient2,
121    ) -> Result<FirstPiolaKirchhoffStress2, ConstitutiveError> {
122        let deformation_gradient_2_inverse = deformation_gradient_2.inverse();
123        let deformation_gradient_1 = deformation_gradient * &deformation_gradient_2_inverse;
124        Ok(FirstPiolaKirchhoffStress2::from(
125            self.1
126                .first_piola_kirchhoff_stress(deformation_gradient_2.into())?,
127        ) - deformation_gradient_1.transpose()
128            * FirstPiolaKirchhoffStress1::from(
129                self.0
130                    .first_piola_kirchhoff_stress(&deformation_gradient_1.into())?,
131            )
132            * deformation_gradient_2_inverse.transpose())
133    }
134    /// Calculates and returns the tangents of the coupled system.
135    ///
136    /// ```math
137    /// \mathcal{C}_{iJkL} = \frac{\partial P_{iJ}}{\partial F_{kL}}
138    /// ```
139    /// ```math
140    /// \frac{\partial R_{IJ}}{\partial F_{kL}} = -F_{IL}^{2-T}P_{kJ} - F_{mI}^1\mathcal{C}_{mJkL}
141    /// ```
142    /// ```math
143    /// \frac{\partial P_{iJ}}{\partial F_{KL}^2} = -P_{iL}F_{KJ}^{2-T} - \mathcal{C}_{iJmL}F_{mK}^1
144    /// ```
145    /// ```math
146    /// \frac{\partial R_{IJ}}{\partial F_{KL}^2} = \mathcal{C}_{IJKL}^2 + F_{IM}^1P_{ML}{F_{KJ}^{2-T}} - \frac{\partial R_{IJ}}{\partial F_{mL}}\,F_{mK}^1
147    /// ```
148    fn tangents(
149        &self,
150        deformation_gradient: &DeformationGradient,
151        deformation_gradient_2: &DeformationGradient2,
152    ) -> Result<
153        (
154            FirstPiolaKirchhoffTangentStiffness,
155            TensorRank4<3, Intermediate, Reference, Current, Reference, Stress>,
156            TensorRank4<3, Current, Reference, Intermediate, Reference, Stress>,
157            FirstPiolaKirchhoffTangentStiffness2,
158        ),
159        ConstitutiveError,
160    > {
161        let deformation_gradient_2_inverse = deformation_gradient_2.inverse();
162        let deformation_gradient_2_inverse_transpose = deformation_gradient_2_inverse.transpose();
163        let deformation_gradient_1 = deformation_gradient * &deformation_gradient_2_inverse;
164        let deformation_gradient_1_transpose = deformation_gradient_1.transpose();
165        let first_piola_kirchhoff_stress =
166            self.first_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_2)?;
167        let tangent_0 = self.first_piola_kirchhoff_tangent_stiffness(
168            deformation_gradient,
169            deformation_gradient_2,
170        )?;
171        let tangent_1 = TensorRank4::dyad_il_kj(
172            &(deformation_gradient_2_inverse_transpose * -1.0),
173            &first_piola_kirchhoff_stress,
174        ) - &deformation_gradient_1_transpose * &tangent_0;
175        let tangent_2 = TensorRank4::dyad_il_jk(
176            &first_piola_kirchhoff_stress,
177            &(&deformation_gradient_2_inverse * -1.0),
178        ) - tangent_0.contract_third_with_first(&deformation_gradient_1);
179        let tangent_3 = FirstPiolaKirchhoffTangentStiffness2::from(
180            self.1
181                .first_piola_kirchhoff_tangent_stiffness(deformation_gradient_2.into())?,
182        ) - tangent_1.contract_third_with_first(&deformation_gradient_1)
183            + TensorRank4::dyad_il_jk(
184                &(deformation_gradient_1_transpose * first_piola_kirchhoff_stress),
185                &deformation_gradient_2_inverse,
186            );
187        Ok((tangent_0, tangent_1, tangent_2, tangent_3))
188    }
189    /// The second deformation gradient is lower triangular to fix rotational freedom.
190    fn internal_variables_fixed(&self) -> &[usize] {
191        &[1, 2, 5]
192    }
193}