Boundary System

TrixiParticles.BoundarySPHSystem — Type

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

System for boundaries modeled by boundary particles. The interaction between fluid and boundary particles is specified by the boundary model.

Arguments

initial_condition: Initial condition (see InitialCondition)
boundary_model: Boundary model (see Boundary Models)

Keyword Arguments

movement: For moving boundaries, a BoundaryMovement can be passed.
adhesion_coefficient: Coefficient specifying the adhesion of a fluid to the surface. Note: currently it is assumed that all fluids have the same adhesion coefficient.

source

TrixiParticles.BoundaryDEMSystem — Type

BoundaryDEMSystem(initial_condition, normal_stiffness)

System for boundaries modeled by boundary particles. The interaction between fluid and boundary particles is specified by the boundary model.

Experimental Implementation

This is an experimental feature and may change in a future releases.

source

TrixiParticles.BoundaryMovement — Type

BoundaryMovement(movement_function, is_moving; moving_particles=nothing)

Arguments

movement_function: Time-dependent function returning an SVector of $d$ dimensions for a $d$-dimensional problem.
is_moving: Function to determine in each timestep if the particles are moving or not. Its boolean return value is mandatory to determine if the neighborhood search will be updated.

Keyword Arguments

moving_particles: Indices of moving particles. Default is each particle in BoundarySPHSystem.

In the example below, movement describes particles moving in a circle as long as the time is lower than 1.5.

Examples

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

movement = BoundaryMovement(movement_function, is_moving)

source

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$.

TrixiParticles.BoundaryModelDummyParticles — Type

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

Boundary model for BoundarySPHSystem.

Arguments

initial_density: Vector holding the initial density of each boundary particle.
hydrodynamic_mass: Vector holding the "hydrodynamic mass" of each boundary particle. See description above for more information.
density_calculator: Strategy to compute the hydrodynamic density of the boundary particles. See description below for more information.
smoothing_kernel: Smoothing kernel should be the same as for the adjacent fluid system.
smoothing_length: Smoothing length should be the same as for the adjacent fluid system.

Keywords

state_equation: This should be the same as for the adjacent fluid system (see e.g. StateEquationCole).
correction: Correction method of the adjacent fluid system (see Corrections).
viscosity: Slip (default) or no-slip condition. See description below for further information.

Examples

# 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))

source

Hydrodynamic density of dummy particles

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

(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.
With SummationDensity, the density is calculated by summation over the neighboring particles, and the pressure is computed from the density with the state equation.
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.
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.
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$.

TrixiParticles.AdamiPressureExtrapolation — Type

AdamiPressureExtrapolation()

density_calculator for BoundaryModelDummyParticles.

source

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).

TrixiParticles.PressureZeroing — Type

PressureZeroing()

density_calculator for BoundaryModelDummyParticles.

Note

This boundary model produces significantly worse results than all other models and is only included for research purposes.

source

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).

TrixiParticles.PressureMirroring — Type

PressureMirroring()

density_calculator for BoundaryModelDummyParticles.

Note

This boundary model requires high viscosity for stability with WCSPH. It also produces significantly worse results than AdamiPressureExtrapolation and is not more efficient because smaller time steps are required due to more noise in the pressure. We added this model only for research purposes and for comparison with SPlisHSPlasH.

source

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.

Warning

The viscosity model ArtificialViscosityMonaghan for BoundaryModelDummyParticles has not been verified yet.

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.

Warning

The no-slip conditions for BoundaryModelMonaghanKajtar have not been verified yet.

TrixiParticles.BoundaryModelMonaghanKajtar — Type

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

Boundary model for BoundarySPHSystem.

Arguments

K: Scaling factor for repulsive force.
beta: Ratio of fluid particle spacing to boundary particle spacing.
boundary_particle_spacing: Boundary particle spacing.
mass: Vector holding the mass of each boundary particle.

Keywords

viscosity: Free-slip (default) or no-slip condition. See description above for further information.

source

Open Boundaries

TrixiParticles.OpenBoundarySPHSystem — Type

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

Open boundary system for in- and outflow particles.

Arguments

boundary_zone: Use InFlow for an inflow and OutFlow for an outflow boundary.

Keywords

fluid_system: The corresponding fluid system
boundary_model: Boundary model (see Open Boundary Models)
buffer_size: Number of buffer particles.
reference_velocity: Reference velocity is either a function mapping each particle's coordinates and time to its velocity, an array where the $i$-th column holds the velocity of particle $i$ or, for a constant fluid velocity, a vector holding this velocity.
reference_pressure: Reference pressure is either a function mapping each particle's coordinates and time to its pressure, a vector holding the pressure of each particle, or a scalar for a constant pressure over all particles.
reference_density: Reference density is either a function mapping each particle's coordinates and time to its density, a vector holding the density of each particle, or a scalar for a constant density over all particles.

Experimental Implementation

This is an experimental feature and may change in any future releases.

source

TrixiParticles.InFlow — Type

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

Inflow boundary zone for OpenBoundarySPHSystem.

The specified plane (line in 2D or rectangle in 3D) will be extruded in upstream direction (the direction opposite to flow_direction) to create a box for the boundary zone. There are three ways to specify the actual shape of the inflow:

Don't pass initial_condition or extrude_geometry. The boundary zone box will then be filled with inflow particles (default).
Specify extrude_geometry by passing a 1D shape in 2D or a 2D shape in 3D, which is then extruded in upstream direction to create the inflow particles.
- In 2D, the shape must be either an initial condition with 2D coordinates, which lies on the line specified by plane, or an initial condition with 1D coordinates, which lies on the line specified by plane when a y-coordinate of 0 is added.
- In 3D, the shape must be either an initial condition with 3D coordinates, which lies in the rectangle specified by plane, or an initial condition with 2D coordinates, which lies in the rectangle specified by plane when a z-coordinate of 0 is added.
Specify initial_condition by passing a 2D initial condition in 2D or a 3D initial condition in 3D, which will be used for the inflow particles.

Note

Particles outside the boundary zone box will be removed.

Keywords

plane: Tuple of points defining a part of the surface of the domain. The points must either span a line in 2D or a rectangle in 3D. This line or rectangle is then extruded in upstream direction to obtain the boundary zone. In 2D, pass two points $(A, B)$, so that the interval $[A, B]$ is the inflow surface. In 3D, pass three points $(A, B, C)$, so that the rectangular inflow surface is spanned by the vectors $\widehat{AB}$ and $\widehat{AC}$. These two vectors must be orthogonal.
flow_direction: Vector defining the flow direction.
open_boundary_layers: Number of particle layers in upstream direction.
particle_spacing: The spacing between the particles (see InitialCondition).
density: Particle density (see InitialCondition).
initial_condition=nothing: InitialCondition for the inflow particles. Particles outside the boundary zone will be removed. Do not use together with extrude_geometry.
extrude_geometry=nothing: 1D shape in 2D or 2D shape in 3D, which lies on the plane and is extruded upstream to obtain the inflow particles. See point 2 above for more details.

Examples

# 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)

Experimental Implementation

This is an experimental feature and may change in any future releases.

source

TrixiParticles.OutFlow — Type

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

Outflow boundary zone for OpenBoundarySPHSystem.

The specified plane (line in 2D or rectangle in 3D) will be extruded in downstream direction (the direction in flow_direction) to create a box for the boundary zone. There are three ways to specify the actual shape of the outflow:

Don't pass initial_condition or extrude_geometry. The boundary zone box will then be filled with outflow particles (default).
Specify extrude_geometry by passing a 1D shape in 2D or a 2D shape in 3D, which is then extruded in downstream direction to create the outflow particles.
- In 2D, the shape must be either an initial condition with 2D coordinates, which lies on the line specified by plane, or an initial condition with 1D coordinates, which lies on the line specified by plane when a y-coordinate of 0 is added.
- In 3D, the shape must be either an initial condition with 3D coordinates, which lies in the rectangle specified by plane, or an initial condition with 2D coordinates, which lies in the rectangle specified by plane when a z-coordinate of 0 is added.
Specify initial_condition by passing a 2D initial condition in 2D or a 3D initial condition in 3D, which will be used for the outflow particles.

Note

Particles outside the boundary zone box will be removed.

Keywords

plane: Tuple of points defining a part of the surface of the domain. The points must either span a line in 2D or a rectangle in 3D. This line or rectangle is then extruded in downstream direction to obtain the boundary zone. In 2D, pass two points $(A, B)$, so that the interval $[A, B]$ is the outflow surface. In 3D, pass three points $(A, B, C)$, so that the rectangular outflow surface is spanned by the vectors $\widehat{AB}$ and $\widehat{AC}$. These two vectors must be orthogonal.
flow_direction: Vector defining the flow direction.
open_boundary_layers: Number of particle layers in downstream direction.
particle_spacing: The spacing between the particles (see InitialCondition).
density: Particle density (see InitialCondition).
initial_condition=nothing: InitialCondition for the outflow particles. Particles outside the boundary zone will be removed. Do not use together with extrude_geometry.
extrude_geometry=nothing: 1D shape in 2D or 2D shape in 3D, which lies on the plane and is extruded downstream to obtain the outflow particles. See point 2 above for more details.

Examples

# 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)

Experimental Implementation

This is an experimental feature and may change in any future releases.

source

Open Boundary Models

Method of characteristics

TrixiParticles.BoundaryModelLastiwka — Type

BoundaryModelLastiwka()

Boundary model for OpenBoundarySPHSystem. This model uses the characteristic variables to propagate the appropriate values to the outlet or inlet and have been proposed by Lastiwka et al. (2009). For more information about the method see description below.

source

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.