carapace_core/model/materials/
shell_section.rs

1//! `ShellSection`: the section seam for 3D shell elements, the shell counterpart of
2//! `PlaneMaterial`. A value holds only *committed* state; `trial_resultant` is pure and
3//! `commit` returns the new committed value. An element owns one independent copy per
4//! Gauss point.
5//!
6//! The generalized strain and resultant vectors follow OpenSees' `ElasticMembranePlateSection`
7//! so recorded values compare directly: `[eps_x, eps_y, gamma_xy, kappa_x, kappa_y,
8//! 2 kappa_xy, gamma_xz, gamma_yz]` and `[Nx, Ny, Nxy, Mx, My, Mxy, Qx, Qy]` in the element's
9//! local axes. OpenSees' bending resultants carry the opposite sign of the usual
10//! `M = D kappa` (its `kappa` is `+w,xx`), so `D`'s bending block is negative; the element
11//! compensates in its `B` matrix.
12
13use std::fmt;
14
15use nalgebra::{SMatrix, SVector};
16
17/// Generalized shell strain or resultant: membrane (3), bending (3), transverse shear (2).
18pub type ShellVector = SVector<f64, 8>;
19/// The `8x8` section tangent.
20pub type ShellMatrix = SMatrix<f64, 8, 8>;
21
22/// Transverse shear correction factor `5/6`.
23const SHEAR_CORRECTION: f64 = 5.0 / 6.0;
24
25#[derive(Debug, Clone, PartialEq)]
26pub enum ShellSection {
27    /// Isotropic linear elastic membrane, plate bending and transverse shear
28    /// (`ElasticMembranePlateSection`): modulus `e`, Poisson ratio `nu`, thickness `h`,
29    /// mass density `rho`.
30    ElasticMembranePlate { e: f64, nu: f64, h: f64, rho: f64 },
31}
32
33/// Why a `ShellSection`'s constants are unusable.
34#[derive(Debug, Clone, Copy, PartialEq)]
35pub enum ShellSectionError {
36    InvalidModulus,
37    InvalidPoissonRatio,
38    InvalidThickness,
39    InvalidDensity,
40}
41
42impl ShellSectionError {
43    pub fn message(&self) -> &'static str {
44        match self {
45            ShellSectionError::InvalidModulus => "modulus is not finite and positive",
46            ShellSectionError::InvalidPoissonRatio => "Poisson ratio must lie in (-1, 0.5]",
47            ShellSectionError::InvalidThickness => "thickness is not finite and positive",
48            ShellSectionError::InvalidDensity => "density is not finite and non-negative",
49        }
50    }
51}
52
53impl fmt::Display for ShellSectionError {
54    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
55        f.write_str(self.message())
56    }
57}
58
59impl std::error::Error for ShellSectionError {}
60
61impl ShellSection {
62    pub fn elastic_membrane_plate(
63        e: f64,
64        nu: f64,
65        h: f64,
66        rho: f64,
67    ) -> Result<Self, ShellSectionError> {
68        let section = ShellSection::ElasticMembranePlate { e, nu, h, rho };
69        section.validate()?;
70        Ok(section)
71    }
72
73    /// Checks the constants of a value built from the variants directly.
74    pub fn validate(&self) -> Result<(), ShellSectionError> {
75        match *self {
76            ShellSection::ElasticMembranePlate { e, nu, h, rho } => {
77                if !(e.is_finite() && e > 0.0) {
78                    return Err(ShellSectionError::InvalidModulus);
79                }
80                if !(nu.is_finite() && nu > -1.0 && nu <= 0.5) {
81                    return Err(ShellSectionError::InvalidPoissonRatio);
82                }
83                if !(h.is_finite() && h > 0.0) {
84                    return Err(ShellSectionError::InvalidThickness);
85                }
86                if !(rho.is_finite() && rho >= 0.0) {
87                    return Err(ShellSectionError::InvalidDensity);
88                }
89                Ok(())
90            }
91        }
92    }
93
94    /// Whether the response is linear (independent of history).
95    pub fn is_linear(&self) -> bool {
96        true
97    }
98
99    /// Mass per unit area, `rho * h`.
100    pub fn rho_h(&self) -> f64 {
101        match *self {
102            ShellSection::ElasticMembranePlate { h, rho, .. } => rho * h,
103        }
104    }
105
106    /// The section tangent (OpenSees sign convention, see the module doc comment).
107    pub fn tangent(&self) -> ShellMatrix {
108        match *self {
109            ShellSection::ElasticMembranePlate { e, nu, h, .. } => {
110                let m = e / (1.0 - nu * nu) * h;
111                let g = 0.5 * e / (1.0 + nu) * h;
112                let d = e * h * h * h / 12.0 / (1.0 - nu * nu);
113                let gs = g * SHEAR_CORRECTION;
114                let mut t = ShellMatrix::zeros();
115                t[(0, 0)] = m;
116                t[(1, 1)] = m;
117                t[(0, 1)] = nu * m;
118                t[(1, 0)] = nu * m;
119                t[(2, 2)] = g;
120                t[(3, 3)] = -d;
121                t[(4, 4)] = -d;
122                t[(3, 4)] = -nu * d;
123                t[(4, 3)] = -nu * d;
124                t[(5, 5)] = -0.5 * d * (1.0 - nu);
125                t[(6, 6)] = gs;
126                t[(7, 7)] = gs;
127                t
128            }
129        }
130    }
131
132    /// Resultants and tangent at a trial strain.
133    pub fn trial_resultant(&self, strain: &ShellVector) -> (ShellVector, ShellMatrix) {
134        let tangent = self.tangent();
135        (tangent * strain, tangent)
136    }
137
138    /// The drilling penalty stiffness OpenSees uses: the smallest eigenvalue of the
139    /// membrane block of the initial tangent.
140    pub fn drilling_stiffness(&self) -> f64 {
141        match *self {
142            ShellSection::ElasticMembranePlate { e, nu, h, .. } => {
143                let m = e / (1.0 - nu * nu) * h;
144                // Eigenvalues of m * [[1, nu, 0], [nu, 1, 0], [0, 0, (1 - nu) / 2]].
145                (m * (1.0 - nu) / 2.0)
146                    .min(m * (1.0 - nu))
147                    .min(m * (1.0 + nu))
148            }
149        }
150    }
151
152    /// The new committed value at `strain`. Elastic sections carry no history.
153    pub fn commit(&self, _strain: &ShellVector) -> Self {
154        self.clone()
155    }
156}
carapace_core/model/materials/shell_section.rs

carapace_core/model/materials/
shell_section.rs
