Skip to main content

conspire/constitutive/solid/elastic/seth_hill_lagrangian/
mod.rs

1#[cfg(test)]
2mod test;
3
4use super::hencky::Hencky;
5use crate::{
6    constitutive::{
7        ConstitutiveError,
8        solid::{Solid, TWO_THIRDS, elastic::Elastic},
9    },
10    math::{
11        ContractThirdFourthWithFirstSecond, IDENTITY_00, Quantity, Rank2, Spectrum, TensorRank4,
12    },
13    mechanics::{
14        Deformation, DeformationGradient, Scalar, SecondPiolaKirchhoffStress,
15        SecondPiolaKirchhoffTangentStiffness,
16    },
17    units::Stress,
18};
19
20#[doc = include_str!("doc.md")]
21#[derive(Clone, Debug)]
22pub struct SethHillLagrangian {
23    /// The bulk modulus $`\kappa`$.
24    pub bulk_modulus: Quantity<Stress>,
25    /// The shear modulus $`\mu`$.
26    pub shear_modulus: Quantity<Stress>,
27    /// The exponent $`m`$.
28    pub exponent: Scalar,
29}
30
31impl SethHillLagrangian {
32    /// Returns the exponent.
33    pub fn exponent(&self) -> Scalar {
34        self.exponent
35    }
36    fn hencky(&self) -> Hencky {
37        Hencky {
38            bulk_modulus: self.bulk_modulus,
39            shear_modulus: self.shear_modulus,
40        }
41    }
42}
43
44impl Solid for SethHillLagrangian {
45    fn bulk_modulus(&self) -> Quantity<Stress> {
46        self.bulk_modulus
47    }
48    fn shear_modulus(&self) -> Quantity<Stress> {
49        self.shear_modulus
50    }
51}
52
53impl Elastic for SethHillLagrangian {
54    #[doc = include_str!("second_piola_kirchhoff_stress.md")]
55    fn second_piola_kirchhoff_stress(
56        &self,
57        deformation_gradient: &DeformationGradient,
58    ) -> Result<SecondPiolaKirchhoffStress, ConstitutiveError> {
59        if self.exponent() == 0.0 {
60            return self
61                .hencky()
62                .second_piola_kirchhoff_stress(deformation_gradient);
63        }
64        let _jacobian = self.jacobian(deformation_gradient)?;
65        let right_cauchy_green = deformation_gradient.right_cauchy_green();
66        let (deviatoric_strain, strain_trace) = ((right_cauchy_green
67            .powm(0.5 * self.exponent())
68            .map_err(|error| ConstitutiveError::upstream(error, self))?
69            - IDENTITY_00)
70            / self.exponent())
71        .deviatoric_and_trace();
72        Ok(deviatoric_strain * (2.0 * self.shear_modulus())
73            + IDENTITY_00 * (self.bulk_modulus() * strain_trace))
74    }
75    #[doc = include_str!("second_piola_kirchhoff_tangent_stiffness.md")]
76    fn second_piola_kirchhoff_tangent_stiffness(
77        &self,
78        deformation_gradient: &DeformationGradient,
79    ) -> Result<SecondPiolaKirchhoffTangentStiffness, ConstitutiveError> {
80        if self.exponent() == 0.0 {
81            return self
82                .hencky()
83                .second_piola_kirchhoff_tangent_stiffness(deformation_gradient);
84        }
85        let _jacobian = self.jacobian(deformation_gradient)?;
86        let right_cauchy_green = deformation_gradient.right_cauchy_green();
87        let half_exponent = 0.5 * self.exponent();
88        let spectrum = Spectrum::new(&right_cauchy_green)
89            .map_err(|error| ConstitutiveError::upstream(error, self))?;
90        let deformation_gradient_transpose = deformation_gradient.transpose();
91        let scaled_deformation_gradient_transpose =
92            &deformation_gradient_transpose * (2.0 * self.shear_modulus() / self.exponent());
93        let trace_term = deformation_gradient
94            * spectrum
95                .powm(half_exponent - 1.0)
96                .map_err(|error| ConstitutiveError::upstream(error, self))?;
97        Ok(spectrum
98            .dpowm(half_exponent)
99            .map_err(|error| ConstitutiveError::upstream(error, self))?
100            .contract_third_fourth_with_first_second(
101                &(TensorRank4::dyad_il_jk(&IDENTITY_00, &scaled_deformation_gradient_transpose)
102                    + TensorRank4::dyad_ik_jl(
103                        &scaled_deformation_gradient_transpose,
104                        &IDENTITY_00,
105                    )),
106            )
107            + TensorRank4::dyad_ij_kl(
108                &(IDENTITY_00 * (self.bulk_modulus() - TWO_THIRDS * self.shear_modulus())),
109                &trace_term,
110            ))
111    }
112}