Boundary System

BoundarySPHSystem(initial_condition, boundary_model; movement=nothing, adhesion_coefficient=0.0)

BoundaryDEMSystem(initial_condition, normal_stiffness)

BoundaryMovement(movement_function, is_moving; moving_particles=nothing)

movement_function(t) = SVector(cos(2pi*t), sin(2pi*t))
is_moving(t) = t < 1.5

movement = BoundaryMovement(movement_function, is_moving)

Boundary Models

Dummy Particles

Boundaries modeled as dummy particles, which are treated like fluid particles, but their positions and velocities are not evolved in time. Since the force towards the fluid should not change with the material density when used with a TotalLagrangianSPHSystem, the dummy particles need to have a mass corresponding to the fluid's rest density, which we call "hydrodynamic mass", as opposed to mass corresponding to the material density of a TotalLagrangianSPHSystem.

Here, initial_density and hydrodynamic_mass are vectors that contains the initial density and the hydrodynamic mass respectively for each boundary particle. Note that when used with SummationDensity (see below), this is only used to determine the element type and the number of boundary particles.

To establish a relationship between density and pressure, a state_equation has to be passed, which should be the same as for the adjacent fluid systems. To sum over neighboring particles, a smoothing_kernel and smoothing_length needs to be passed. This should be the same as for the adjacent fluid system with the largest smoothing length.

In the literature, this kind of boundary particles is referred to as "dummy particles" (Adami et al., 2012 and Valizadeh & Monaghan, 2015), "frozen fluid particles" (Akinci et al., 2012) or "dynamic boundaries (Crespo et al., 2007). The key detail of this boundary condition and the only difference between the boundary models in these references is the way the density and pressure of boundary particles is computed.

Since boundary particles are treated like fluid particles, the force on fluid particle $a$ due to boundary particle $b$ is given by

\[f_{ab} = m_a m_b \left( \frac{p_a}{\rho_a^2} + \frac{p_b}{\rho_b^2} \right) \nabla_{r_a} W(\Vert r_a - r_b \Vert, h).\]

The quantities to be defined here are the density $\rho_b$ and pressure $p_b$ of the boundary particle $b$.

BoundaryModelDummyParticles(initial_density, hydrodynamic_mass,
                            density_calculator, smoothing_kernel,
                            smoothing_length; viscosity=nothing,
                            state_equation=nothing, correction=nothing)

# Free-slip condition
boundary_model = BoundaryModelDummyParticles(densities, masses, AdamiPressureExtrapolation(),
                                             smoothing_kernel, smoothing_length)

# No-slip condition
boundary_model = BoundaryModelDummyParticles(densities, masses, AdamiPressureExtrapolation(),
                                             smoothing_kernel, smoothing_length,
                                             viscosity=ViscosityAdami(nu=1e-6))

Hydrodynamic density of dummy particles

We provide five options to compute the boundary density and pressure, determined by the density_calculator:

1. (Recommended) With AdamiPressureExtrapolation, the pressure is extrapolated from the pressure of the fluid according to (Adami et al., 2012), and the density is obtained by applying the inverse of the state equation. This option usually yields the best results of the options listed here.
2. With SummationDensity, the density is calculated by summation over the neighboring particles, and the pressure is computed from the density with the state equation.
3. With ContinuityDensity, the density is integrated from the continuity equation, and the pressure is computed from the density with the state equation. Note that this causes a gap between fluid and boundary where the boundary is initialized without any contact to the fluid. This is due to overestimation of the boundary density as soon as the fluid comes in contact with boundary particles that initially did not have contact to the fluid. Therefore, in dam break simulations, there is a visible "step", even though the boundary is supposed to be flat. See also dual.sphysics.org/faq/#Q_13.
4. With PressureZeroing, the density is set to the reference density and the pressure is computed from the density with the state equation. This option is not recommended. The other options yield significantly better results.
5. With PressureMirroring, the density is set to the reference density. The pressure is not used. Instead, the fluid pressure is mirrored as boundary pressure in the momentum equation. This option is not recommended due to stability issues. See PressureMirroring for more details.

1. AdamiPressureExtrapolation

The pressure of the boundary particles is obtained by extrapolating the pressure of the fluid according to (Adami et al., 2012). The pressure of a boundary particle $b$ is given by

\[p_b = \frac{\sum_f (p_f + \rho_f (\bm{g} - \bm{a}_b) \cdot \bm{r}_{bf}) W(\Vert r_{bf} \Vert, h)}{\sum_f W(\Vert r_{bf} \Vert, h)},\]

where the sum is over all fluid particles, $\rho_f$ and $p_f$ denote the density and pressure of fluid particle $f$, respectively, $r_{bf} = r_b - r_f$ denotes the difference of the coordinates of particles $b$ and $f$, $\bm{g}$ denotes the gravitational acceleration acting on the fluid, and $\bm{a}_b$ denotes the acceleration of the boundary particle $b$.

AdamiPressureExtrapolation()

4. PressureZeroing

This is the simplest way to implement dummy boundary particles. The density of each particle is set to the reference density and the pressure to the reference pressure (the corresponding pressure to the reference density by the state equation).

PressureZeroing()

5. PressureMirroring

Instead of calculating density and pressure for each boundary particle, we modify the momentum equation,

\[\frac{\mathrm{d}v_a}{\mathrm{d}t} = -\sum_b m_b \left( \frac{p_a}{\rho_a^2} + \frac{p_b}{\rho_b^2} \right) \nabla_a W_{ab}\]

to replace the unknown density $\rho_b$ if $b$ is a boundary particle by the reference density and the unknown pressure $p_b$ if $b$ is a boundary particle by the pressure $p_a$ of the interacting fluid particle. The momentum equation therefore becomes

\[\frac{\mathrm{d}v_a}{\mathrm{d}t} = -\sum_f m_f \left( \frac{p_a}{\rho_a^2} + \frac{p_f}{\rho_f^2} \right) \nabla_a W_{af} -\sum_b m_b \left( \frac{p_a}{\rho_a^2} + \frac{p_a}{\rho_0^2} \right) \nabla_a W_{ab},\]

where the first sum is over all fluid particles and the second over all boundary particles.

This approach was first mentioned by Akinci et al. (2012) and written down in this form by Band et al. (2018).

PressureMirroring()

No-slip conditions

For the interaction of dummy particles and fluid particles, Adami et al. (2012) impose a no-slip boundary condition by assigning a wall velocity $v_w$ to the dummy particle.

The wall velocity of particle $a$ is calculated from the prescribed boundary particle velocity $v_a$ and the smoothed velocity field

\[v_w = 2 v_a - \frac{\sum_b v_b W_{ab}}{\sum_b W_{ab}},\]

where the sum is over all fluid particles.

By choosing the viscosity model ViscosityAdami for viscosity, a no-slip condition is imposed. It is recommended to choose nu in the order of either the kinematic viscosity parameter of the adjacent fluid or the equivalent from the artificial parameter alpha of the adjacent fluid ($\nu = \frac{\alpha h c }{2d + 4}$). When omitting the viscous interaction (default viscosity=nothing), a free-slip wall boundary condition is applied.

References

• S. Adami, X. Y. Hu, N. A. Adams. "A generalized wall boundary condition for smoothed particle hydrodynamics". In: Journal of Computational Physics 231, 21 (2012), pages 7057–7075. doi: 10.1016/J.JCP.2012.05.005
• Alireza Valizadeh, Joseph J. Monaghan. "A study of solid wall models for weakly compressible SPH". In: Journal of Computational Physics 300 (2015), pages 5–19. doi: 10.1016/J.JCP.2015.07.033
• Nadir Akinci, Markus Ihmsen, Gizem Akinci, Barbara Solenthaler, Matthias Teschner. "Versatile rigid-fluid coupling for incompressible SPH". ACM Transactions on Graphics 31, 4 (2012), pages 1–8. doi: 10.1145/2185520.2185558
• A. J. C. Crespo, M. Gómez-Gesteira, R. A. Dalrymple. "Boundary conditions generated by dynamic particles in SPH methods" In: Computers, Materials and Continua 5 (2007), pages 173-184. doi: 10.3970/cmc.2007.005.173
• Stefan Band, Christoph Gissler, Andreas Peer, and Matthias Teschner. "MLS Pressure Boundaries for Divergence-Free and Viscous SPH Fluids." In: Computers & Graphics 76 (2018), pages 37–46. doi: 10.1016/j.cag.2018.08.001

Repulsive Particles

Boundaries modeled as boundary particles which exert forces on the fluid particles (Monaghan, Kajtar, 2009). The force on fluid particle $a$ due to boundary particle $b$ is given by

\[f_{ab} = m_a \left(\tilde{f}_{ab} - m_b \Pi_{ab} \nabla_{r_a} W(\Vert r_a - r_b \Vert, h)\right)\]

with

\[\tilde{f}_{ab} = \frac{K}{\beta^{n-1}} \frac{r_{ab}}{\Vert r_{ab} \Vert (\Vert r_{ab} \Vert - d)} \Phi(\Vert r_{ab} \Vert, h) \frac{2 m_b}{m_a + m_b},\]

where $m_a$ and $m_b$ are the masses of fluid particle $a$ and boundary particle $b$ respectively, $r_{ab} = r_a - r_b$ is the difference of the coordinates of particles $a$ and $b$, $d$ denotes the boundary particle spacing and $n$ denotes the number of dimensions (see (Monaghan, Kajtar, 2009, Equation (3.1)) and (Valizadeh, Monaghan, 2015)). Note that the repulsive acceleration $\tilde{f}_{ab}$ does not depend on the masses of the boundary particles. Here, $\Phi$ denotes the 1D Wendland C4 kernel, normalized to $1.77$ for $q=0$ (Monaghan, Kajtar, 2009, Section 4), with $\Phi(r, h) = w(r/h)$ and

\[w(q) = \begin{cases} (1.77/32) (1 + (5/2)q + 2q^2)(2 - q)^5 & \text{if } 0 \leq q < 2 \\ 0 & \text{if } q \geq 2. \end{cases}\]

The boundary particles are assumed to have uniform spacing by the factor $\beta$ smaller than the expected fluid particle spacing. For example, if the fluid particles have an expected spacing of $0.3$ and the boundary particles have a uniform spacing of $0.1$, then this parameter should be set to $\beta = 3$. According to (Monaghan, Kajtar, 2009), a value of $\beta = 3$ for the Wendland C4 that we use here is reasonable for most computing purposes.

The parameter $K$ is used to scale the force exerted by the boundary particles. In (Monaghan, Kajtar, 2009), a value of $gD$ is used for static tank simulations, where $g$ is the gravitational acceleration and $D$ is the depth of the fluid.

The viscosity $\Pi_{ab}$ is calculated according to the viscosity used in the simulation, where the density of the boundary particle if needed is assumed to be identical to the density of the fluid particle.

No-slip condition

By choosing the viscosity model ArtificialViscosityMonaghan for viscosity, a no-slip condition is imposed. When omitting the viscous interaction (default viscosity=nothing), a free-slip wall boundary condition is applied.

BoundaryModelMonaghanKajtar(K, beta, boundary_particle_spacing, mass;
                            viscosity=nothing)

References

• Joseph J. Monaghan, Jules B. Kajtar. "SPH particle boundary forces for arbitrary boundaries". In: Computer Physics Communications 180.10 (2009), pages 1811–1820. doi: 10.1016/j.cpc.2009.05.008
• Alireza Valizadeh, Joseph J. Monaghan. "A study of solid wall models for weakly compressible SPH." In: Journal of Computational Physics 300 (2015), pages 5–19. doi: 10.1016/J.JCP.2015.07.033

Open Boundaries

OpenBoundarySPHSystem(boundary_zone::Union{InFlow, OutFlow}; sound_speed,
                      fluid_system::FluidSystem, buffer_size::Integer,
                      boundary_model,
                      reference_velocity=nothing,
                      reference_pressure=nothing,
                      reference_density=nothing)

InFlow(; plane, flow_direction, density, particle_spacing,
       initial_condition=nothing, extrude_geometry=nothing,
       open_boundary_layers::Integer)

# 2D
plane_points = ([0.0, 0.0], [0.0, 1.0])
flow_direction=[1.0, 0.0]

inflow = InFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0,
                open_boundary_layers=4)

# 3D
plane_points = ([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0])
flow_direction=[0.0, 0.0, 1.0]

inflow = InFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0,
                open_boundary_layers=4)

# 3D particles sampled as cylinder
circle = SphereShape(0.1, 0.5, (0.5, 0.5), 1.0, sphere_type=RoundSphere())

inflow = InFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0,
                extrude_geometry=circle, open_boundary_layers=4)

OutFlow(; plane, flow_direction, density, particle_spacing,
        initial_condition=nothing, extrude_geometry=nothing,
        open_boundary_layers::Integer)

# 2D
plane_points = ([0.0, 0.0], [0.0, 1.0])
flow_direction = [1.0, 0.0]

outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0,
                  open_boundary_layers=4)

# 3D
plane_points = ([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0])
flow_direction = [0.0, 0.0, 1.0]

outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0,
                  open_boundary_layers=4)

# 3D particles sampled as cylinder
circle = SphereShape(0.1, 0.5, (0.5, 0.5), 1.0, sphere_type=RoundSphere())

outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0,
                  extrude_geometry=circle, open_boundary_layers=4)

Open Boundary Models

Method of characteristics

BoundaryModelLastiwka()

The difficulty in non-reflecting boundary conditions, also called open boundaries, is to determine the appropriate boundary values of the exact characteristics of the Euler equations. Assuming the flow near the boundaries is normal to the boundary and free of shock waves and significant viscous effects, it can be shown that three characteristic variables exist:

• $J_1$, associated with convection of entropy and propagates at flow velocity,
• $J_2$, downstream-running characteristics,
• $J_3$, upstream-running characteristics.

Giles (1990) derived those variables based on a linearized set of governing equations:

\[J_1 = -c_s^2 (\rho - \rho_{\text{ref}}) + (p - p_{\text{ref}})\]

\[J_2 = \rho c_s (v - v_{\text{ref}}) + (p - p_{\text{ref}})\]

\[J_3 = - \rho c_s (v - v_{\text{ref}}) + (p - p_{\text{ref}})\]

where the subscript "ref" denotes the reference flow near the boundaries, which can be prescribed.

Specifying the reference variables is not equivalent to prescription of $\rho$, $v$ and $p$ directly, since the perturbation from the reference flow is allowed.

Lastiwka et al. (2009) applied the method of characteristic to SPH and determined the number of variables that should be prescribed at the boundary and the number which should be propagated from the fluid domain to the boundary:

• For an inflow boundary:

  • Prescribe downstream-running characteristics $J_1$ and $J_2$
  • Transmit $J_3$ from the fluid domain (allow $J_3$ to propagate upstream to the boundary).

• For an outflow boundary:

  • Prescribe upstream-running characteristic $J_3$
  • Transmit $J_1$ and $J_2$ from the fluid domain.

Prescribing is done by simply setting the characteristics to zero. To transmit the characteristics from the fluid domain, or in other words, to carry the information of the fluid to the boundaries, Negi et al. (2020) use a Shepard Interpolation

\[f_i = \frac{\sum_j^N f_j W_{ij}}{\sum_j^N W_{ij}},\]

where the $i$-th particle is a boundary particle, $f$ is either $J_1$, $J_2$ or $J_3$ and $N$ is the set of neighboring fluid particles.

To express pressure $p$, density $\rho$ and velocity $v$ as functions of the characteristic variables, the system of equations from the characteristic variables is inverted and gives

\[ \rho - \rho_{\text{ref}} = \frac{1}{c_s^2} \left( -J_1 + \frac{1}{2} J_2 + \frac{1}{2} J_3 \right),\]

\[u - u_{\text{ref}}= \frac{1}{2\rho c_s} \left( J_2 - J_3 \right),\]

\[p - p_{\text{ref}} = \frac{1}{2} \left( J_2 + J_3 \right).\]

With $J_1$, $J_2$ and $J_3$ determined, we can easily solve for the actual variables for each particle.

References

• M. B. Giles. "Nonreflecting boundary conditions for Euler equation calculations". In: AIAA Journal, 28.12 pages 2050–2058. doi: 10.2514/3.10521
• M. Lastiwka, M. Basa, N. J. Quinlan. "Permeable and non-reflecting boundary conditions in SPH". In: International Journal for Numerical Methods in Fluids 61, (2009), pages 709–724. doi: 10.1002/fld.1971
• P. Negi, P. Ramachandran, A. Haftu. "An improved non-reflecting outlet boundary condition for weakly-compressible SPH". In: Computer Methods in Applied Mechanics and Engineering 367, (2020), pages 113–119. doi: 10.1016/j.cma.2020.113119