salva2d/solver/surface_tension/
wcsph_surface_tension.rs

1#[cfg(feature = "parallel")]
2use rayon::prelude::*;
3
4use crate::geometry::ParticlesContacts;
5
6use crate::math::Real;
7use crate::object::{Boundary, Fluid};
8use crate::solver::NonPressureForce;
9use crate::TimestepManager;
10
11// Surface tension of water: 0.01
12// Stable values of surface tension: up to 3.4
13// From https://cg.informatik.uni-freiburg.de/publications/2007_SCA_SPH.pdf
14/// Surface tension method introduced by the WCSPH method.
15#[derive(Clone)]
16pub struct WCSPHSurfaceTension {
17    fluid_tension_coefficient: Real,
18    boundary_tension_coefficient: Real,
19}
20
21impl WCSPHSurfaceTension {
22    /// Initializes a surface tension with the given surface tension coefficient and boundary adhesion coefficients.
23    pub fn new(fluid_tension_coefficient: Real, boundary_tension_coefficient: Real) -> Self {
24        Self {
25            fluid_tension_coefficient,
26            boundary_tension_coefficient,
27        }
28    }
29}
30
31impl NonPressureForce for WCSPHSurfaceTension {
32    fn solve(
33        &mut self,
34        _timestep: &TimestepManager,
35        _kernel_radius: Real,
36        fluid_fluid_contacts: &ParticlesContacts,
37        _fluid_boundaries_contacts: &ParticlesContacts,
38        fluid: &mut Fluid,
39        boundaries: &[Boundary],
40        _densities: &[Real],
41    ) {
42        let fluid_tension_coefficient = self.fluid_tension_coefficient;
43        let boundary_tension_coefficient = self.boundary_tension_coefficient;
44        let positions = &fluid.positions;
45        let volumes = &fluid.volumes;
46        let density0 = fluid.density0;
47
48        par_iter_mut!(fluid.accelerations)
49            .enumerate()
50            .for_each(|(i, acceleration_i)| {
51                if fluid_tension_coefficient != na::zero::<Real>() {
52                    for c in fluid_fluid_contacts
53                        .particle_contacts(i)
54                        .read()
55                        .unwrap()
56                        .iter()
57                    {
58                        if c.i_model == c.j_model {
59                            let dpos = positions[c.i] - positions[c.j];
60                            let cohesion_acc = dpos
61                                * (-fluid_tension_coefficient * c.weight * volumes[c.j] * density0
62                                    / (volumes[c.i] * density0));
63                            *acceleration_i += cohesion_acc;
64                        }
65                    }
66                }
67
68                if boundary_tension_coefficient != na::zero::<Real>() {
69                    for c in fluid_fluid_contacts
70                        .particle_contacts(i)
71                        .read()
72                        .unwrap()
73                        .iter()
74                    {
75                        let dpos = positions[c.i] - boundaries[c.j_model].positions[c.j];
76                        let mi = volumes[c.i] * density0;
77                        let cohesion_force = dpos
78                            * (boundary_tension_coefficient
79                                * c.weight
80                                * boundaries[c.j_model].volumes[c.j]
81                                * density0);
82                        *acceleration_i -= cohesion_force / mi;
83                        boundaries[c.j_model].apply_force(c.j, cohesion_force);
84                    }
85                }
86            })
87    }
88
89    fn apply_permutation(&mut self, _: &[usize]) {}
90}
salva2d/solver/surface_tension/wcsph_surface_tension.rs

salva2d/solver/surface_tension/
wcsph_surface_tension.rs
