carapace_core/model/materials/
concrete02.rs

1use super::Material;
2
3impl Material {
4    /// `fc`/`epsc0`/`fcu`/`epscu` are normalized to be negative
5    /// (compression), matching OpenSees' `Concrete02` constructor. `rat` is
6    /// the ratio between the unloading slope at `epscu` and the initial
7    /// slope; `ft`/`ets` are the tensile strength and tension-softening
8    /// modulus.
9    pub fn concrete02(
10        fc: f64,
11        epsc0: f64,
12        fcu: f64,
13        epscu: f64,
14        rat: f64,
15        ft: f64,
16        ets: f64,
17    ) -> Self {
18        let fc = -fc.abs();
19        let epsc0 = -epsc0.abs();
20        let fcu = -fcu.abs();
21        let epscu = -epscu.abs();
22        Material::Concrete02 {
23            fc,
24            epsc0,
25            fcu,
26            epscu,
27            rat,
28            ft,
29            ets,
30            min_strain: 0.0,
31            tension_strain: 0.0,
32            strain: 0.0,
33            stress: 0.0,
34            tangent: 2.0 * fc / epsc0,
35        }
36    }
37}
38
39/// Ported from `Concrete02::setTrialStrain` (`Concrete02.cpp`).
40pub(super) fn evaluate(m: &Material, strain: f64) -> (f64, f64, Material) {
41    let Material::Concrete02 {
42        fc,
43        epsc0,
44        fcu,
45        epscu,
46        rat,
47        ft,
48        ets,
49        min_strain,
50        tension_strain,
51        strain: cstrain,
52        stress: cstress,
53        tangent: ctangent,
54    } = m
55    else {
56        unreachable!()
57    };
58    let (fc, epsc0, fcu, epscu, rat, ft, ets) = (*fc, *epsc0, *fcu, *epscu, *rat, *ft, *ets);
59    let (mut min_strain, mut tension_strain) = (*min_strain, *tension_strain);
60    let (cstrain, cstress) = (*cstrain, *cstress);
61
62    let ec0 = 2.0 * fc / epsc0;
63
64    // `Concrete02::setTrialStrain` returns early below `DBL_EPSILON`.
65    if (strain - cstrain).abs() < f64::EPSILON {
66        return (cstress, *ctangent, m.clone());
67    }
68
69    let (stress, tangent) = if strain < min_strain {
70        // Further into compression than ever committed: monotonic
71        // compression envelope, and a new low-water mark.
72        let (stress, tangent) = concrete02_compr_envlp(strain, fc, epsc0, fcu, epscu);
73        min_strain = strain;
74        (stress, tangent)
75    } else {
76        // Unloading/reloading branch: point R (`reload_strain`,
77        // `reload_stress`) is where the reloading slope, extended from
78        // the previous minimum, would cross the *initial* tangent line
79        // — the classic EERC-report construction (`Concrete02.cpp`'s
80        // comments cite Eq. 2.31-2.36).
81        let reload_strain = (fcu - rat * ec0 * epscu) / (ec0 * (1.0 - rat));
82        let reload_stress = ec0 * reload_strain;
83        let (min_stress, _) = concrete02_compr_envlp(min_strain, fc, epsc0, fcu, epscu);
84        let reload_slope = (min_stress - reload_stress) / (min_strain - reload_strain);
85        // Strain where the reloading line (through the previous
86        // minimum) crosses zero stress.
87        let zero_stress_strain = min_strain - min_stress / reload_slope;
88
89        if strain <= zero_stress_strain {
90            let stress_lower_bound = min_stress + reload_slope * (strain - min_strain);
91            let stress_upper_bound = reload_slope * 0.5 * (strain - zero_stress_strain);
92            let mut stress = cstress + ec0 * (strain - cstrain);
93            let mut tangent = ec0;
94            if stress <= stress_lower_bound {
95                stress = stress_lower_bound;
96                tangent = reload_slope;
97            }
98            if stress >= stress_upper_bound {
99                stress = stress_upper_bound;
100                tangent = 0.5 * reload_slope;
101            }
102            (stress, tangent)
103        } else {
104            // Tension side: `tension_peak_strain` is the strain at the
105            // peak of the tensile stress-strain relation reached so
106            // far, shifted by `zero_stress_strain` (Eq. 2.42-2.43).
107            let tension_peak_strain = zero_stress_strain + tension_strain;
108            if strain <= tension_peak_strain {
109                let (peak_stress, _) = concrete02_tens_envlp(tension_strain, ft, ec0, ets);
110                let tangent = if tension_strain != 0.0 {
111                    peak_stress / tension_strain
112                } else {
113                    ec0
114                };
115                (tangent * (strain - zero_stress_strain), tangent)
116            } else {
117                tension_strain = strain - zero_stress_strain;
118                concrete02_tens_envlp(tension_strain, ft, ec0, ets)
119            }
120        }
121    };
122
123    (
124        stress,
125        tangent,
126        Material::Concrete02 {
127            fc,
128            epsc0,
129            fcu,
130            epscu,
131            rat,
132            ft,
133            ets,
134            min_strain,
135            tension_strain,
136            strain,
137            stress,
138            tangent,
139        },
140    )
141}
142
143/// `Concrete02::Compr_Envlp`: same Kent-Scott-Park compression envelope as
144/// `Concrete01`'s (see `concrete01::concrete01_envelope`), but note the flat
145/// crushed branch's tangent floor here is `1e-10` (matching `Concrete02.cpp`),
146/// not the exact `0.0` `concrete01_envelope` uses for the same branch — a
147/// genuine difference between the two ports, not a simplification to
148/// reconcile away.
149fn concrete02_compr_envlp(strain: f64, fc: f64, epsc0: f64, fcu: f64, epscu: f64) -> (f64, f64) {
150    let ec0 = 2.0 * fc / epsc0;
151    let ratio = strain / epsc0;
152    if strain >= epsc0 {
153        (fc * ratio * (2.0 - ratio), ec0 * (1.0 - ratio))
154    } else if strain > epscu {
155        let tangent = (fcu - fc) / (epscu - epsc0);
156        (tangent * (strain - epsc0) + fc, tangent)
157    } else {
158        (fcu, 1.0e-10)
159    }
160}
161
162/// `Concrete02::Tens_Envlp`: linear elastic up to `ft`, then linear
163/// softening at slope `-ets` down to zero, then a flat zero branch (with
164/// the same `1e-10` tangent floor as the compression envelope above, for
165/// the same singular-tangent-avoidance reason).
166fn concrete02_tens_envlp(strain: f64, ft: f64, ec0: f64, ets: f64) -> (f64, f64) {
167    let eps0 = ft / ec0;
168    let epsu = ft * (1.0 / ets + 1.0 / ec0);
169    if strain <= eps0 {
170        (strain * ec0, ec0)
171    } else if strain <= epsu {
172        (ft - ets * (strain - eps0), -ets)
173    } else {
174        (0.0, 1.0e-10)
175    }
176}
177
178#[cfg(test)]
179mod tests {
180    use super::*;
181
182    #[test]
183    fn concrete02_matches_compression_envelope_on_first_loading() {
184        let m = Material::concrete02(4.0, 0.002, 3.0, 0.006, 0.1, 0.4, 50.0);
185        let strain = -0.001;
186        let (stress, tangent) = m.trial_stress_tangent(strain);
187        let eta = strain / -0.002;
188        let expected_stress = -4.0 * (2.0 * eta - eta * eta);
189        assert!(
190            (stress - expected_stress).abs() < 1e-9,
191            "expected {expected_stress}, got {stress}"
192        );
193        assert!((tangent - 2000.0).abs() < 1e-9, "got {tangent}");
194    }
195
196    #[test]
197    fn concrete02_is_linear_elastic_in_tension_below_ft() {
198        let m = Material::concrete02(4.0, 0.002, 3.0, 0.006, 0.1, 0.4, 50.0);
199        let ec0 = 2.0 * 4.0 / 0.002;
200        // ft/Ec0 = 0.4/4000 = 0.0001, so 0.00005 is within the linear
201        // range — unlike Concrete01, Concrete02 carries real tension.
202        let (stress, tangent) = m.trial_stress_tangent(0.00005);
203        assert!((stress - ec0 * 0.00005).abs() < 1e-9, "got {stress}");
204        assert!((tangent - ec0).abs() < 1e-9, "got {tangent}");
205    }
206
207    #[test]
208    fn concrete02_softens_in_tension_past_ft() {
209        let m = Material::concrete02(4.0, 0.002, 3.0, 0.006, 0.1, 0.4, 50.0);
210        // Past eps0=ft/Ec0=0.0001: linear softening at slope -Ets from ft.
211        let (stress, tangent) = m.trial_stress_tangent(0.0005);
212        let expected_stress = 0.4 - 50.0 * (0.0005 - 0.0001);
213        assert!(
214            (stress - expected_stress).abs() < 1e-9,
215            "expected {expected_stress}, got {stress}"
216        );
217        assert_eq!(tangent, -50.0);
218    }
219
220    #[test]
221    fn concrete02_unloads_with_degrading_stiffness_below_initial() {
222        let m = Material::concrete02(4.0, 0.002, 3.0, 0.006, 0.1, 0.4, 50.0);
223        let ec0 = 2.0 * 4.0 / 0.002;
224
225        let m = m.commit(-0.003);
226        let (_stress, unload_tangent) = m.trial_stress_tangent(-0.002);
227        assert!(
228            unload_tangent < ec0 && unload_tangent > 0.0,
229            "reload slope should be positive but degraded below Ec0={ec0}, got {unload_tangent}"
230        );
231    }
232}
carapace_core/model/materials/concrete02.rs

carapace_core/model/materials/
concrete02.rs
