sparkl2d_core/dynamics/models/
plasticity_nacc.rs

1use crate::math::{Matrix, Real, Vector, DIM};
2use na::vector;
3
4#[cfg(not(feature = "std"))]
5use na::ComplexField;
6
7/// The Non-Associated-Cam-Clay plasticity model.
8#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
9#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
10#[derive(Copy, Clone, Debug, PartialEq)]
11#[repr(C)]
12pub struct NaccPlasticity {
13    mu: Real,
14    kappa: Real,
15    hardening_enabled: bool,
16    hardening_factor: Real, // xi
17    cohesion: Real,         // beta
18    friction: Real,         // M
19}
20
21impl NaccPlasticity {
22    pub fn new(
23        young_modulus: Real,
24        poisson_ratio: Real,
25        cohesion: Real,
26        hardening_enabled: bool,
27        hardening_factor: Real,
28        friction_angle: Real,
29    ) -> Self {
30        let sin_f = friction_angle.sin();
31        let d = DIM as Real;
32        Self {
33            mu: crate::utils::shear_modulus(young_modulus, poisson_ratio),
34            kappa: crate::utils::bulk_modulus(young_modulus, poisson_ratio),
35            hardening_enabled,
36            hardening_factor,
37            cohesion,
38            friction: (2.0 as Real / 3.0).sqrt() * 2.0 * sin_f / (3.0 - sin_f) * d
39                / (2.0 / (6.0 - d)).sqrt(),
40        }
41    }
42
43    pub fn with_m(
44        young_modulus: Real,
45        poisson_ratio: Real,
46        cohesion: Real,
47        hardening_enabled: bool,
48        hardening_factor: Real,
49        m: Real,
50    ) -> Self {
51        Self {
52            mu: crate::utils::shear_modulus(young_modulus, poisson_ratio),
53            kappa: crate::utils::bulk_modulus(young_modulus, poisson_ratio),
54            hardening_enabled,
55            hardening_factor,
56            cohesion,
57            friction: m,
58        }
59    }
60
61    pub fn project_deformation_gradient(
62        &self,
63        deformation_gradient: Matrix<Real>,
64        mut alpha: Real,
65    ) -> (Matrix<Real>, Real) {
66        let xi = self.hardening_factor;
67        let beta = self.cohesion;
68        let m = self.friction;
69        let d = DIM as Real;
70
71        let mut svd = deformation_gradient.svd_unordered(true, true);
72
73        let square_eigv = svd.singular_values.component_mul(&svd.singular_values);
74        let square_eigv_trace = square_eigv.sum();
75
76        let p0 = self.kappa * (1.0e-5 + ((xi * (-alpha).max(0.0)).sinh()));
77        #[cfg(feature = "dim2")]
78        let j_e_tr = svd.singular_values[0] * svd.singular_values[1];
79        #[cfg(feature = "dim3")]
80        let j_e_tr = svd.singular_values[0] * svd.singular_values[1] * svd.singular_values[2];
81        let s_tr =
82            self.mu * j_e_tr.powf(-2.0 / d) * (square_eigv - Vector::repeat(square_eigv_trace / d));
83        let psi_kappa = self.kappa / 2.0 * (j_e_tr - 1.0 / j_e_tr);
84        let p_tr = -psi_kappa * j_e_tr;
85
86        if p_tr > p0 {
87            // Project to max tip of the yield surface.
88            let j_e_n1 = (-2.0 * p0 / self.kappa + 1.0).sqrt();
89            svd.singular_values.fill(j_e_n1.powf(1.0 / d));
90
91            if self.hardening_enabled {
92                alpha += (j_e_tr / j_e_n1).ln();
93            }
94            // info!("Return a: {} > {}", p_tr, p0);
95            return (svd.recompose().unwrap(), alpha);
96        }
97
98        if p_tr < -beta * p0 {
99            // Project to min tip of the yield surface.
100            let j_e_n1 = (2.0 * beta * p0 / self.kappa + 1.0).sqrt();
101            svd.singular_values.fill(j_e_n1.powf(1.0 / d));
102
103            if self.hardening_enabled {
104                alpha += (j_e_tr / j_e_n1).ln();
105            }
106            // info!("Retrun B");
107            return (svd.recompose().unwrap(), alpha);
108        }
109
110        let y0 = (1.0 + 2.0 * beta) * ((6.0 - d) / 2.0);
111        let y1 = m * m * (p_tr + beta * p0) * (p_tr - p0);
112        let y = y0 * s_tr.norm_squared() + y1;
113
114        if y < 1.0e-4 {
115            // Inside the yield surface.
116            return (deformation_gradient, alpha);
117        }
118
119        if self.hardening_enabled && p0 > 1.0e-4 && p_tr < p0 - 1.0e-4 && p_tr > -beta * p0 + 1.0e-4
120        {
121            // Hardening routine.
122            let p_c = (1.0 - beta) * p0 / 2.0;
123            let q_tr = ((6.0 - d) / 2.0).sqrt() * s_tr.norm();
124            let direction = vector![p_c - p_tr, 0.0 - q_tr];
125            let direction = direction.normalize();
126            let c = m * m * (p_c + beta * p0) * (p_c - p0);
127            let b = m * m * direction[0] * (2.0 * p_c - p0 + beta * p0);
128            let a = m * m * direction[0] * direction[0]
129                + (1.0 + 2.0 * beta) * direction[1] * direction[1];
130            let discr = (b * b - 4.0 * a * c).sqrt();
131            let l1 = (-b + discr) / (2.0 * a);
132            let l2 = (-b - discr) / (2.0 * a);
133            let p1 = p_c + l1 * direction[0];
134            let p2 = p_c + l2 * direction[0];
135            let p_x = if (p_tr - p_c) * (p1 - p_c) > 0.0 {
136                p1
137            } else {
138                p2
139            };
140            let j_e_x = (-2.0 * p_x / self.kappa + 1.0).abs().sqrt();
141
142            if j_e_x > 1.0e-4 {
143                alpha += (j_e_tr / j_e_x).ln();
144            }
145        }
146
147        // Yield surface projection.
148        let b_e_n1 = (-y1 / y0).sqrt() * (j_e_tr.powf(2.0 / d) / self.mu) * s_tr.normalize()
149            + Vector::repeat(square_eigv_trace / d);
150
151        svd.singular_values = b_e_n1.map(|e| e.sqrt());
152        // info!("Retrun D");
153        return (svd.recompose().unwrap(), alpha);
154    }
155
156    pub fn update_particle(
157        &self,
158        particle_deformation_gradient: &mut Matrix<Real>,
159        particle_nacc_alpha: &mut Real,
160    ) {
161        let (new_def, new_alpha) =
162            self.project_deformation_gradient(*particle_deformation_gradient, *particle_nacc_alpha);
163        *particle_deformation_gradient = new_def;
164        *particle_nacc_alpha = new_alpha;
165    }
166}
sparkl2d_core/dynamics/models/plasticity_nacc.rs

sparkl2d_core/dynamics/models/
plasticity_nacc.rs
