sparkl3d_core/dynamics/models/
eos_monaghan_sph.rs

1use crate::dynamics::models::ActiveTimestepBounds;
2use crate::math::{DecomposedTensor, Matrix, Real, Vector};
3
4#[cfg(not(feature = "std"))]
5use na::ComplexField;
6
7#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
8#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
9#[derive(Copy, Clone, Debug, PartialEq)]
10#[repr(C)]
11pub struct MonaghanSphEos {
12    pub pressure0: Real,
13    pub gamma: u32,
14    pub viscosity: Real,
15    pub max_neg_pressure: Real,
16}
17
18impl MonaghanSphEos {
19    pub fn new(pressure0: Real, gamma: u32, viscosity: Real) -> Self {
20        Self {
21            pressure0,
22            gamma,
23            viscosity,
24            max_neg_pressure: 1.0,
25        }
26    }
27
28    pub fn pressure(
29        &self,
30        particle_mass: Real,
31        particle_volume0: Real,
32        particle_density_fluid: Real,
33    ) -> Real {
34        let density0 = particle_mass / particle_volume0;
35        (self.pressure0 * ((particle_density_fluid / density0).powi(self.gamma as i32) - 1.0))
36            .max(-self.max_neg_pressure)
37    }
38
39    pub fn is_fluid(&self) -> bool {
40        true
41    }
42
43    pub fn kirchhoff_stress(
44        &self,
45        particle_mass: Real,
46        particle_volume0: Real,
47        particle_density_fluid: Real,
48        particle_fluid_deformation_gradient_det: Real,
49        particle_velocity_gradient: &Matrix<Real>,
50    ) -> Matrix<Real> {
51        let mut stress = (-self.pressure(particle_mass, particle_volume0, particle_density_fluid)
52            * particle_fluid_deformation_gradient_det)
53            * Matrix::identity();
54
55        if self.viscosity != 0.0 {
56            let strain_rate = crate::utils::strain_rate(particle_velocity_gradient);
57            let strain_rate = DecomposedTensor::decompose(&strain_rate);
58            stress += (2.0 * self.viscosity * particle_fluid_deformation_gradient_det)
59                * strain_rate.deviatoric_part;
60        }
61
62        stress
63    }
64
65    pub fn active_timestep_bounds(&self) -> ActiveTimestepBounds {
66        ActiveTimestepBounds::CONSTITUTIVE_MODEL_BOUND
67            | ActiveTimestepBounds::PARTICLE_VELOCITY_BOUND
68            | ActiveTimestepBounds::PARTICLE_DISPLACEMENT_BOUND
69            | ActiveTimestepBounds::SINGLE_PARTICLE_STABILITY_BOUND
70    }
71
72    pub fn timestep_bound(
73        &self,
74        particle_fluid_deformation_gradient_det: Real,
75        particle_mass: Real,
76        particle_volume0: Real,
77        particle_density_fluid: Real,
78        particle_velocity: &Vector<Real>,
79        cell_width: Real,
80    ) -> Real {
81        // Single-particle criteria.
82        let j = particle_fluid_deformation_gradient_det;
83        let density0 = particle_mass / particle_volume0;
84        let k = 6.0; // For quadratic splines.
85        let d = crate::math::DIM as Real;
86        let pressure = -self.pressure(particle_mass, particle_volume0, particle_density_fluid);
87
88        let single_particle_dt =
89            (cell_width / j) * (density0 * (j - 1.0) / (k * pressure * d)).sqrt();
90
91        // CFL
92        let density_fluctuation = 0.1;
93        let c_square = particle_velocity.norm_squared().max(1.0) / density_fluctuation;
94        let cfl_dt = cell_width / c_square.sqrt();
95
96        // Result
97        single_particle_dt.min(cfl_dt)
98        // cfl_dt
99    }
100}
sparkl3d_core/dynamics/models/eos_monaghan_sph.rs

sparkl3d_core/dynamics/models/
eos_monaghan_sph.rs
