carapace_core/model/materials/
steel01.rs

1use super::Material;
2
3#[derive(Debug, Clone, Copy, PartialEq, Eq)]
4pub enum Steel01Loading {
5    None,
6    Loading,
7    Unloading,
8}
9
10impl Material {
11    pub fn steel01(fy: f64, e0: f64, b: f64, a1: f64, a2: f64, a3: f64, a4: f64) -> Self {
12        Material::Steel01 {
13            fy,
14            e0,
15            b,
16            a1,
17            a2,
18            a3,
19            a4,
20            min_strain: 0.0,
21            max_strain: 0.0,
22            shift_p: 1.0,
23            shift_n: 1.0,
24            loading: Steel01Loading::None,
25            strain: 0.0,
26            stress: 0.0,
27            tangent: e0,
28        }
29    }
30}
31
32/// Ported from `Steel01::determineTrialState` (`Steel01.cpp`). `m` stands in
33/// for the `C*` fields; `strain` is the trial (`Tstrain`).
34pub(super) fn evaluate(m: &Material, strain: f64) -> (f64, f64, Material) {
35    let Material::Steel01 {
36        fy,
37        e0,
38        b,
39        a1,
40        a2,
41        a3,
42        a4,
43        min_strain,
44        max_strain,
45        shift_p,
46        shift_n,
47        loading,
48        strain: cstrain,
49        stress: cstress,
50        tangent: ctangent,
51    } = m
52    else {
53        unreachable!()
54    };
55    let (fy, e0, b, a1, a2, a3, a4) = (*fy, *e0, *b, *a1, *a2, *a3, *a4);
56    let (mut min_strain, mut max_strain, mut shift_p, mut shift_n) =
57        (*min_strain, *max_strain, *shift_p, *shift_n);
58    let mut loading = *loading;
59    let (cstrain, cstress) = (*cstrain, *cstress);
60
61    let dstrain = strain - cstrain;
62
63    // `Steel01::setTrialStrain` skips `determineTrialState` below
64    // `DBL_EPSILON`, leaving the trial state at the committed one.
65    if dstrain.abs() <= f64::EPSILON {
66        return (cstress, *ctangent, m.clone());
67    }
68
69    let fy_one_minus_b = fy * (1.0 - b);
70    let esh = b * e0;
71    let epsy = fy / e0;
72
73    let c1 = esh * strain;
74    let c2 = shift_n * fy_one_minus_b;
75    let c3 = shift_p * fy_one_minus_b;
76    let c = cstress + e0 * dstrain;
77
78    // The three-way clamp `max(c1-c2, min(c1+c3, c))` — see the
79    // commented-out simplification in `Steel01.cpp`, which the source
80    // says only exists split into two `if`s as an optimizer
81    // workaround; the single-expression form is used here.
82    let stress = (c1 - c2).max((c1 + c3).min(c));
83    let tangent = if (stress - c).abs() < f64::EPSILON {
84        e0
85    } else {
86        esh
87    };
88
89    if loading == Steel01Loading::None && dstrain != 0.0 {
90        loading = if dstrain > 0.0 {
91            Steel01Loading::Loading
92        } else {
93            Steel01Loading::Unloading
94        };
95    }
96
97    // Transition from loading to unloading.
98    if loading == Steel01Loading::Loading && dstrain < 0.0 {
99        loading = Steel01Loading::Unloading;
100        if cstrain > max_strain {
101            max_strain = cstrain;
102        }
103        shift_n = 1.0 + a1 * ((max_strain - min_strain) / (2.0 * a2 * epsy)).powf(0.8);
104    }
105
106    // Transition from unloading to loading.
107    if loading == Steel01Loading::Unloading && dstrain > 0.0 {
108        loading = Steel01Loading::Loading;
109        if cstrain < min_strain {
110            min_strain = cstrain;
111        }
112        shift_p = 1.0 + a3 * ((max_strain - min_strain) / (2.0 * a4 * epsy)).powf(0.8);
113    }
114
115    (
116        stress,
117        tangent,
118        Material::Steel01 {
119            fy,
120            e0,
121            b,
122            a1,
123            a2,
124            a3,
125            a4,
126            min_strain,
127            max_strain,
128            shift_p,
129            shift_n,
130            loading,
131            strain,
132            stress,
133            tangent,
134        },
135    )
136}
137
138#[cfg(test)]
139mod tests {
140    use super::*;
141
142    #[test]
143    fn steel01_matches_elastic_on_first_loading_within_yield() {
144        let m = Material::steel01(60.0, 29000.0, 0.01, 0.9, 5.0, 0.9, 5.0);
145        let (stress, tangent) = m.trial_stress_tangent(0.001);
146        assert!((stress - 29000.0 * 0.001).abs() < 1e-9, "got {stress}");
147        assert_eq!(tangent, 29000.0);
148    }
149
150    #[test]
151    fn steel01_bilinear_hardening_past_yield() {
152        let m = Material::steel01(60.0, 29000.0, 0.01, 0.9, 5.0, 0.9, 5.0);
153        // epsy = fy/E0 = 60/29000 ≈ 0.002069. Push well past it.
154        let (stress, tangent) = m.trial_stress_tangent(0.02);
155        let esh = 0.01 * 29000.0;
156        // On first loading (shift_p=shift_n=1, min=max=0) the upper
157        // bounding line is c1+c3 = Esh*strain + fy*(1-b).
158        let expected = esh * 0.02 + 60.0 * (1.0 - 0.01);
159        assert!(
160            (stress - expected).abs() < 1e-6,
161            "expected {expected}, got {stress}"
162        );
163        assert!((tangent - esh).abs() < 1e-9);
164    }
165
166    #[test]
167    fn steel01_trial_at_committed_strain_keeps_committed_tangent() {
168        // A displacement-control predictor evaluates the committed state; a
169        // yielded fiber must report its hardening tangent, not e0.
170        let m = Material::steel01(60.0, 29000.0, 0.01, 0.9, 5.0, 0.9, 5.0).commit(0.02);
171        let (stress, tangent) = m.trial_stress_tangent(0.02);
172        assert_eq!(tangent, 0.01 * 29000.0);
173        assert!((stress - (0.01 * 29000.0 * 0.02 + 60.0 * 0.99)).abs() < 1e-9);
174    }
175
176    #[test]
177    fn steel01_isotropic_hardening_shifts_the_bound_on_reversal() {
178        let m = Material::steel01(60.0, 29000.0, 0.01, 0.9, 5.0, 0.9, 5.0);
179        // Load well into the hardening range, commit, then reverse twice —
180        // the first reversal (loading -> unloading) shifts shift_n, the
181        // second (unloading -> loading) shifts shift_p.
182        let m = m.commit(0.02);
183        let m = m.commit(-0.02);
184        let Material::Steel01 { shift_n, .. } = m else {
185            panic!()
186        };
187        assert!(
188            shift_n > 1.0,
189            "shift_n should have grown after the first reversal: {shift_n}"
190        );
191
192        let m = m.commit(0.02);
193        let Material::Steel01 { shift_p, .. } = m else {
194            panic!()
195        };
196        assert!(
197            shift_p > 1.0,
198            "shift_p should have grown after the second reversal: {shift_p}"
199        );
200    }
201}
carapace_core/model/materials/steel01.rs

carapace_core/model/materials/
steel01.rs
