carapace_core/model/materials/
concrete01.rs

1use super::Material;
2
3impl Material {
4    /// `fpc`/`epsc0`/`fpcu`/`epscu` are normalized to be negative
5    /// (compression), matching OpenSees' `Concrete01` constructor.
6    pub fn concrete01(fpc: f64, epsc0: f64, fpcu: f64, epscu: f64) -> Self {
7        let fpc = -fpc.abs();
8        let epsc0 = -epsc0.abs();
9        let fpcu = -fpcu.abs();
10        let epscu = -epscu.abs();
11        let ec0 = 2.0 * fpc / epsc0;
12        Material::Concrete01 {
13            fpc,
14            epsc0,
15            fpcu,
16            epscu,
17            min_strain: 0.0,
18            end_strain: 0.0,
19            unload_slope: ec0,
20            strain: 0.0,
21            stress: 0.0,
22            tangent: ec0,
23        }
24    }
25}
26
27/// Ported from `Concrete01::setTrialStrain` (`Concrete01.cpp`), the
28/// version actually driving the reload/envelope/unload state machine
29/// (`determineTrialState` is dead code in the original — never called
30/// from `setTrialStrain`).
31pub(super) fn evaluate(m: &Material, strain: f64) -> (f64, f64, Material) {
32    let Material::Concrete01 {
33        fpc,
34        epsc0,
35        fpcu,
36        epscu,
37        min_strain,
38        end_strain,
39        unload_slope,
40        strain: cstrain,
41        stress: cstress,
42        tangent: ctangent,
43    } = m
44    else {
45        unreachable!()
46    };
47    let (fpc, epsc0, fpcu, epscu) = (*fpc, *epsc0, *fpcu, *epscu);
48    let (min_strain, end_strain, unload_slope) = (*min_strain, *end_strain, *unload_slope);
49    let (cstrain, cstress) = (*cstrain, *cstress);
50
51    // `Concrete01::setTrialStrain` returns the committed state below
52    // `DBL_EPSILON`.
53    if (strain - cstrain).abs() < f64::EPSILON {
54        return (cstress, *ctangent, m.clone());
55    }
56
57    // Quick return: tension, zero stiffness, no compression-envelope
58    // history update — but strain/stress still commit to the new
59    // point.
60    if strain > 0.0 {
61        return (
62            0.0,
63            0.0,
64            Material::Concrete01 {
65                fpc,
66                epsc0,
67                fpcu,
68                epscu,
69                min_strain,
70                end_strain,
71                unload_slope,
72                strain,
73                stress: 0.0,
74                tangent: 0.0,
75            },
76        );
77    }
78
79    let temp_stress = cstress + unload_slope * strain - unload_slope * cstrain;
80
81    let (stress, tangent, min_strain, end_strain, unload_slope) = if strain < cstrain {
82        // Further into compression: `reload()`.
83        let (mut stress, mut tangent, min_strain, end_strain, unload_slope) = if strain
84            <= min_strain
85        {
86            let new_min_strain = strain;
87            let (stress, tangent) = concrete01_envelope(new_min_strain, fpc, epsc0, fpcu, epscu);
88            let (end_strain, unload_slope) =
89                concrete01_unload(new_min_strain, stress, fpc, epsc0, epscu);
90            (stress, tangent, new_min_strain, end_strain, unload_slope)
91        } else if strain <= end_strain {
92            let tangent = unload_slope;
93            let stress = tangent * (strain - end_strain);
94            (stress, tangent, min_strain, end_strain, unload_slope)
95        } else {
96            (0.0, 0.0, min_strain, end_strain, unload_slope)
97        };
98
99        if temp_stress > stress {
100            stress = temp_stress;
101            tangent = unload_slope;
102        }
103
104        (stress, tangent, min_strain, end_strain, unload_slope)
105    } else if temp_stress <= 0.0 {
106        (
107            temp_stress,
108            unload_slope,
109            min_strain,
110            end_strain,
111            unload_slope,
112        )
113    } else {
114        (0.0, 0.0, min_strain, end_strain, unload_slope)
115    };
116
117    (
118        stress,
119        tangent,
120        Material::Concrete01 {
121            fpc,
122            epsc0,
123            fpcu,
124            epscu,
125            min_strain,
126            end_strain,
127            unload_slope,
128            strain,
129            stress,
130            tangent,
131        },
132    )
133}
134
135/// `Concrete01::envelope()`: the Kent-Scott-Park backbone, evaluated at a
136/// strain at or past the previous minimum (i.e. loading further into
137/// compression along the virgin curve).
138fn concrete01_envelope(strain: f64, fpc: f64, epsc0: f64, fpcu: f64, epscu: f64) -> (f64, f64) {
139    if strain > epsc0 {
140        let eta = strain / epsc0;
141        let stress = fpc * (2.0 * eta - eta * eta);
142        let ec0 = 2.0 * fpc / epsc0;
143        let tangent = ec0 * (1.0 - eta);
144        (stress, tangent)
145    } else if strain > epscu {
146        let tangent = (fpc - fpcu) / (epsc0 - epscu);
147        let stress = fpc + tangent * (strain - epsc0);
148        (stress, tangent)
149    } else {
150        (fpcu, 0.0)
151    }
152}
153
154/// `Concrete01::unload()`: the Karsan-Jirsa degrading unload-slope formula,
155/// given the new `min_strain` and the envelope stress just computed there.
156fn concrete01_unload(
157    min_strain: f64,
158    stress_at_min: f64,
159    fpc: f64,
160    epsc0: f64,
161    epscu: f64,
162) -> (f64, f64) {
163    let mut temp_strain = min_strain;
164    if temp_strain < epscu {
165        temp_strain = epscu;
166    }
167    let eta = temp_strain / epsc0;
168
169    let mut ratio = 0.707 * (eta - 2.0) + 0.834;
170    if eta < 2.0 {
171        ratio = 0.145 * eta * eta + 0.13 * eta;
172    }
173
174    let mut end_strain = ratio * epsc0;
175    let temp1 = min_strain - end_strain;
176    let ec0 = 2.0 * fpc / epsc0;
177    let temp2 = stress_at_min / ec0;
178
179    let unload_slope;
180    if temp1 > -f64::EPSILON {
181        // Should always be negative in practice.
182        unload_slope = ec0;
183    } else if temp1 <= temp2 {
184        end_strain = min_strain - temp1;
185        unload_slope = stress_at_min / temp1;
186    } else {
187        end_strain = min_strain - temp2;
188        unload_slope = ec0;
189    }
190
191    (end_strain, unload_slope)
192}
193
194#[cfg(test)]
195mod tests {
196    use super::*;
197
198    #[test]
199    fn concrete01_matches_parabolic_envelope_on_first_loading() {
200        let m = Material::concrete01(4.0, 0.002, 3.0, 0.006);
201        // Envelope is symmetric-normalized negative; strain -0.001 is
202        // within the parabola (|strain| < |epsc0|).
203        let strain = -0.001;
204        let (stress, _tangent) = m.trial_stress_tangent(strain);
205        let eta = strain / -0.002;
206        let expected = -4.0 * (2.0 * eta - eta * eta);
207        assert!(
208            (stress - expected).abs() < 1e-9,
209            "expected {expected}, got {stress}"
210        );
211    }
212
213    #[test]
214    fn concrete01_carries_no_tension() {
215        let m = Material::concrete01(4.0, 0.002, 3.0, 0.006);
216        assert_eq!(m.trial_stress_tangent(0.001), (0.0, 0.0));
217    }
218
219    #[test]
220    fn concrete01_unloads_with_degrading_stiffness_below_initial() {
221        let m = Material::concrete01(4.0, 0.002, 3.0, 0.006);
222        let ec0 = 2.0 * 4.0 / 0.002;
223
224        // Push past peak stress (epsc0=-0.002) then commit, then partially
225        // unload — but stay on the unload line (between min_strain and
226        // end_strain), not far enough to cross back through zero stress.
227        let m = m.commit(-0.003);
228        let (_stress, unload_tangent) = m.trial_stress_tangent(-0.002);
229        assert!(
230            unload_tangent < ec0 && unload_tangent > 0.0,
231            "Karsan-Jirsa unload slope should be positive but degraded below Ec0={ec0}, got {unload_tangent}"
232        );
233    }
234}
carapace_core/model/materials/concrete01.rs

carapace_core/model/materials/
concrete01.rs
