conspire/constitutive/solid/hyperelastic/eight_chain/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::{
5    constitutive::{
6        ConstitutiveError,
7        solid::{FIVE_THIRDS, Solid, TWO_THIRDS, elastic::Elastic, hyperelastic::Hyperelastic},
8    },
9    math::{IDENTITY, Rank2},
10    mechanics::{CauchyStress, CauchyTangentStiffness, Deformation, DeformationGradient, Scalar},
11    physics::molecular::single_chain::Thermodynamics as SingleChainThermodynamics,
12};
13use std::{
14    any::type_name,
15    fmt::{self, Debug, Formatter},
16};
17
18#[doc = include_str!("doc.md")]
19#[derive(Clone)]
20pub struct EightChain<T>
21where
22    T: SingleChainThermodynamics,
23{
24    /// The bulk modulus $`\kappa`$.
25    pub bulk_modulus: Scalar,
26    /// The shear modulus $`\mu`$.
27    pub shear_modulus: Scalar,
28    /// The single-chain model.
29    pub single_chain_model: T,
30}
31
32impl<T> Debug for EightChain<T>
33where
34    T: SingleChainThermodynamics,
35{
36    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
37        write!(
38            f,
39            "EightChain {{ bulk_modulus: {}, shear_modulus: {}, single_chain_model: {} }}",
40            self.bulk_modulus,
41            self.shear_modulus,
42            type_name::<T>()
43                .rsplit("::")
44                .next()
45                .unwrap()
46                .split("<")
47                .next()
48                .unwrap(),
49        )
50    }
51}
52
53impl<T> EightChain<T>
54where
55    T: SingleChainThermodynamics,
56{
57    /// Returns the nondimensional force in a single chain.
58    fn nondimensional_force(
59        &self,
60        nondimensional_extension: Scalar,
61    ) -> Result<Scalar, ConstitutiveError> {
62        match SingleChainThermodynamics::nondimensional_force(
63            &self.single_chain_model,
64            nondimensional_extension,
65        ) {
66            Ok(nondimensional_force) => Ok(nondimensional_force),
67            Err(error) => Err(ConstitutiveError::Upstream(
68                format!("{error}"),
69                format!("{self:?}"),
70            )),
71        }
72    }
73    /// Returns the nondimensional stiffness in a single chain.
74    fn nondimensional_stiffness(
75        &self,
76        nondimensional_extension: Scalar,
77    ) -> Result<Scalar, ConstitutiveError> {
78        match SingleChainThermodynamics::nondimensional_stiffness(
79            &self.single_chain_model,
80            nondimensional_extension,
81        ) {
82            Ok(nondimensional_stiffness) => Ok(nondimensional_stiffness),
83            Err(error) => Err(ConstitutiveError::Upstream(
84                format!("{error}"),
85                format!("{self:?}"),
86            )),
87        }
88    }
89    /// Returns the number of links in a single chain.
90    pub fn number_of_links(&self) -> Scalar {
91        self.single_chain_model.number_of_links() as Scalar
92    }
93    /// Returns the single-chain model.
94    pub fn single_chain_model(&self) -> &T {
95        &self.single_chain_model
96    }
97}
98
99impl<T> Solid for EightChain<T>
100where
101    T: SingleChainThermodynamics,
102{
103    fn bulk_modulus(&self) -> Scalar {
104        self.bulk_modulus
105    }
106    fn shear_modulus(&self) -> Scalar {
107        self.shear_modulus
108    }
109}
110
111impl<T> Elastic for EightChain<T>
112where
113    T: SingleChainThermodynamics,
114{
115    #[doc = include_str!("cauchy_stress.md")]
116    fn cauchy_stress(
117        &self,
118        deformation_gradient: &DeformationGradient,
119    ) -> Result<CauchyStress, ConstitutiveError> {
120        let jacobian = self.jacobian(deformation_gradient)?;
121        let (
122            deviatoric_isochoric_left_cauchy_green_deformation,
123            isochoric_left_cauchy_green_deformation_trace,
124        ) = (deformation_gradient.left_cauchy_green() / jacobian.powf(TWO_THIRDS))
125            .deviatoric_and_trace();
126        let gamma =
127            (isochoric_left_cauchy_green_deformation_trace / 3.0 / self.number_of_links()).sqrt();
128        let gamma_0 = (1.0 / self.number_of_links()).sqrt();
129        Ok(deviatoric_isochoric_left_cauchy_green_deformation
130            * (self.shear_modulus() * self.nondimensional_force(gamma)?
131                / self.nondimensional_force(gamma_0)?
132                * gamma_0
133                / gamma
134                / jacobian)
135            + IDENTITY * self.bulk_modulus() * 0.5 * (jacobian - 1.0 / jacobian))
136    }
137    #[doc = include_str!("cauchy_tangent_stiffness.md")]
138    fn cauchy_tangent_stiffness(
139        &self,
140        deformation_gradient: &DeformationGradient,
141    ) -> Result<CauchyTangentStiffness, ConstitutiveError> {
142        let jacobian = self.jacobian(deformation_gradient)?;
143        let inverse_transpose_deformation_gradient = deformation_gradient.inverse_transpose();
144        let left_cauchy_green_deformation = deformation_gradient.left_cauchy_green();
145        let deviatoric_left_cauchy_green_deformation = left_cauchy_green_deformation.deviatoric();
146        let (
147            deviatoric_isochoric_left_cauchy_green_deformation,
148            isochoric_left_cauchy_green_deformation_trace,
149        ) = (left_cauchy_green_deformation / jacobian.powf(TWO_THIRDS)).deviatoric_and_trace();
150        let gamma =
151            (isochoric_left_cauchy_green_deformation_trace / 3.0 / self.number_of_links()).sqrt();
152        let gamma_0 = (1.0 / self.number_of_links()).sqrt();
153        let eta = self.nondimensional_force(gamma)?;
154        let scaled_shear_modulus =
155            gamma_0 / self.nondimensional_force(gamma_0)? * self.shear_modulus() * eta
156                / gamma
157                / jacobian.powf(FIVE_THIRDS);
158        let scaled_deviatoric_isochoric_left_cauchy_green_deformation =
159            deviatoric_left_cauchy_green_deformation * scaled_shear_modulus;
160        let term = CauchyTangentStiffness::dyad_ij_kl(
161            &scaled_deviatoric_isochoric_left_cauchy_green_deformation,
162            &(deviatoric_isochoric_left_cauchy_green_deformation
163                * &inverse_transpose_deformation_gradient
164                * ((self.nondimensional_stiffness(gamma)? / eta - 1.0 / gamma)
165                    / 3.0
166                    / self.number_of_links()
167                    / gamma)),
168        );
169        Ok(
170            (CauchyTangentStiffness::dyad_ik_jl(&IDENTITY, deformation_gradient)
171                + CauchyTangentStiffness::dyad_il_jk(deformation_gradient, &IDENTITY)
172                - CauchyTangentStiffness::dyad_ij_kl(&IDENTITY, deformation_gradient)
173                    * (TWO_THIRDS))
174                * scaled_shear_modulus
175                + CauchyTangentStiffness::dyad_ij_kl(
176                    &(IDENTITY * (0.5 * self.bulk_modulus() * (jacobian + 1.0 / jacobian))
177                        - scaled_deviatoric_isochoric_left_cauchy_green_deformation
178                            * (FIVE_THIRDS)),
179                    &inverse_transpose_deformation_gradient,
180                )
181                + term,
182        )
183    }
184}
185
186impl<T> Hyperelastic for EightChain<T>
187where
188    T: SingleChainThermodynamics,
189{
190    #[doc = include_str!("helmholtz_free_energy_density.md")]
191    fn helmholtz_free_energy_density(
192        &self,
193        deformation_gradient: &DeformationGradient,
194    ) -> Result<Scalar, ConstitutiveError> {
195        let jacobian = self.jacobian(deformation_gradient)?;
196        let isochoric_left_cauchy_green_deformation =
197            deformation_gradient.left_cauchy_green() / jacobian.powf(TWO_THIRDS);
198        let gamma =
199            (isochoric_left_cauchy_green_deformation.trace() / 3.0 / self.number_of_links()).sqrt();
200        let eta = self.nondimensional_force(gamma)?;
201        let gamma_0 = (1.0 / self.number_of_links()).sqrt();
202        let eta_0 = self.nondimensional_force(gamma_0)?;
203        //
204        // If end up re-using so much of ArrudaBoyce should make helper function to share.
205        //
206        Ok(3.0 * gamma_0 / eta_0
207            * self.shear_modulus()
208            * self.number_of_links()
209            * (gamma * eta - gamma_0 * eta_0 - (eta_0 * eta.sinh() / (eta * eta_0.sinh())).ln())
210            + 0.5 * self.bulk_modulus() * (0.5 * (jacobian.powi(2) - 1.0) - jacobian.ln()))
211    }
212}