3. SPH formulation
Smoothed Particle Hydrodynamics (SPH) is a Lagrangian meshless method. The technique discretises a continuum using a set of material points or particles. When used for the simulation of fluid dynamics, the discretised Navier-Stokes equations are locally integrated at the location of each of these particles, according to the physical properties of surrounding particles. The set of neighbouring particles is determined by a distance based function, either circular (two-dimensional) or spherical (three-dimensional), with an associated characteristic length or smoothing length often denoted as h. At each timestep new physical quantities are calculated for each particle, and they then move according to the updated values.
The conservation laws of continuum fluid dynamics are transformed from their partial differential form to a form suitable for particle based simulation using integral equations based on an interpolation function, which gives an estimate of values at a specific point.
Typically this interpolation or weighting function is referred to as the kernel function (W) and can take different forms, with the most common being cubic or quintic. In all cases however, it is designed to represent a function F(r) defined in r' by the integral approximation
The smoothing kernel must fulfil several properties [Monaghan, 1992; Liu, 2003], such as positivity inside a defined zone of interaction, compact support, normalization and monotonically decreasing value with distance and differentiability. For a more complete description of SPH, the reader is referred to [Monaghan, 2005; Violeau, 2013].
The function F in Eq. (1) can be approximated in a non-continuous, discrete form based on the set of particles. In this case the function is interpolated at a particle (a) where a summation is performed over all the particles that fall within its region of compact support, as defined by the smoothing length h
where the subscript denotes an individual particle, vb is the volume of a neighbouring particle (b). If vb=mb/ρb, with m and ρ being the mass and the density of particle b respectively then Eq. (2) becomes
The performance of an SPH model depends heavily on the choice of the smoothing kernel. Kernels are expressed as a function of the non-dimensional distance between particles (q), given by q=r/h , where r is the distance between any two given particles a and b and the parameter h (the smoothing length) controls the size of the area around particle a in which neighbouring particles are considered. Within DualSPHysics, the user is able to choose from one of the following kernel definitions:
- Cubic Spline [Monaghan and Lattanzio, 1985]
where αD is equal to 10/7πh2 in 2-D and 1/πh3 in 3-D.
The tensile correction method, proposed by [Monaghan, 2000], is only actively used in the cases of a kernel whose first derivative goes to zero with the particle distance q.
- Quintic [Wendland, 1995]
In the text that follows, only kernels with an influence domain of 2 h are considered.
The momentum conservation equation in a continuum is
where Γ refers to dissipative terms and g is gravitational acceleration. DualSPHysics offers different options for including the effects of dissipation.
The artificial viscosity scheme, proposed by [Monaghan, 1992], is a common method within fluid simulation using SPH due primarily to its simplicity. In SPH notation, Eq. 6 can be written as
where Pk and ρk are the pressure and density that correspond to particle k (as evaluated
at a or b). The viscosity term 
where ab r_ab=r_a-r_b and v_ab=va-vb with r_k and v_k being the particle position and velocity respectively. 
Laminar viscous stresses in the momentum equation can be expressed as [Lo and Shao, 2002]
where υo is kinematic viscosity (typically 10^(-6) for water). In SPH discrete notation this can be expressed as
The concept of the Sub-Particle Scale (SPS) was first described by [Gotoh et al., 2001] to represent the effects of turbulence in their Moving Particle Semi-implicit (MPS) model. The momentum conservation equation is defined as
where the laminar term is treated as per Eq. 9 and 
Throughout the duration of a weakly-compressible SPH simulation (as presented herein) the mass of each particle remains constant and only their associated density fluctuates. These density changes are computed by solving the conservation of mass, or continuity equation, in SPH form:
Following the work of [Monaghan, 1994], the fluid in the SPH formalism defined in DualSPHysics is treated as weakly compressible and an equation of state is used to determine fluid pressure based on particle density. The compressibility is adjusted so that the speed of sound can be artificially lowered; this means that the size of time step taken at any one moment (which is determined according to a Courant condition, based on the currently calculated speed of sound for all particles) can be maintained at a reasonable value. Such adjustment however, restricts the sound speed to be at least ten times faster than the maximum fluid velocity, keeping density variations to within less than 1%, and therefore not introducing major deviations from an incompressible approach. Following [Monaghan et al., 1999] and [Batchelor, 1974], the relationship between pressure and density follows the expression
where 
Within DualSPHysics it is also possible to apply a density diffusion term, that introduces a diffusive term to reduce density fluctuations. The state equation describes a very stiff density field, and together with the natural disordering of the particles, high-frequency low amplitude oscillations are found to populate the density scalar field. DualSPHysics uses a diffusive term in the continuity equation, now written as
with,
This represents the density diffusion formulation by [Molteni and Colagrossi, 2009], with a free parameter that needs to be attributed a suitable value. This modification can be explained as the addition of the Laplacian of the density field to the continuity equation. [Antuono et al., 2012] has presented a careful analysis of the influence of this term in the system, by decomposing the Laplacian operator, observing the converge of the operators and performing linear stability analysis to inspect the influence of the diffusive coefficient. This equation represents exactly a diffusive term in the domain bulk. The behaviour changes close to open boundaries such as free-surface. Due to truncation of the kernel (there are no particles being sampled outside of an open boundary), the first-order contributions are not null [Antuono et al., 2010], resulting in a net force applied to the particles. This effect is not considered relevant for nonhydrostatic situations, where this force is many orders of magnitude inferior to any other force involved. Corrections to this effect were proposed by [Antuono et al., 2010], but involve the solution of a renormalization problem for the density gradient, with considerable computational cost. A (δΦ) coefficient of 0.1 is recommended for most applications.
The work of [Fourtakas et al., 2020] introduces a correction in Eq. (15.2). The key idea is to use the same formulation proposed by [Molteni and Colagrossi, 2009] but substituting the dynamic density with the total one. Thus, the term ψab takes the following form,
where the superscript D denotes the dynamic density or pressure. Since ρD=ρT-ρH (where superscript T and H denote the total and hydrostatic component, respectively), Eq. (16) can be rewritten in term of the hydrostatic and total parts as,
In the context of the weakly compressible SPH, the equation of state relates the density to the total pressure at the particle location. However, only the hydrostatic part of the pressure is needed. Using the equation of state (EOS), the hydrostatic density difference can be obtained by,
The term PH is simply the hydrostatic pressure difference of particle a and b,
where zab is the vertical distance between particle a and b and
In comparison with the formulation of the diffusive term proposed by [Molteni and Colagrossi, 2009] of Eq. (15.2), Eq. (16) improves the behaviour of pressure near the wall boundaries by avoiding the troublesome dynamic pressure and avoiding to compute the normalized density gradient. However, it must be noted that the formulation of [Antuono et al., 2010] is consistent and it can be adopted for any type of flow, whereas the use of total density in the diffusive term of Eq. (16) is expected to work better than Eq. (15.2) mostly for gravity-dominated flows. In the absence of a gravity force the formulation collapses to the original formulation of [Molteni and Colagrossi, 2009].
Anisotropic particle spacing is an important stability issue in SPH as, especially in violent flows, particles cannot maintain a uniform distribution. The result is the introduction of noise in the velocity and pressure field, as well as the creation of voids within the water flow for certain cases.
To counter the anisotropic particle spacing, [Xu et al., 2009] proposed a particle shifting algorithm to prevent the instabilities. The algorithm was first created for incompressible SPH, but can be extended to the weakly compressible SPH model used in DualSPHysics [Vacondio et al., 2013]. With the shifting algorithm, the particles are moved (“shifted”) towards areas with fewer particles (lower particle concentration) allowing the domain to maintain a uniform particle distribution and eliminating any voids that may occur due to the noise.
An improvement on the initial shifting algorithm was proposed by [Lind et al., 2012] who used Fick’s first law of diffusion to control the shifting magnitude and direction. Fick’s first law connects the diffusion flux to the concentration gradient:
where J is the flux, C the particle concentration, and DF the Fickian diffusion coefficient.
Assuming that the flux, i.e. the number of particles passing through a unit surface in unit time, is proportional to the velocity of the particles, a particle shifting velocity and subsequently a particle shifting distance can be found. Using the particle concentration, the particle shifting distance δrs is given by:
where D is a new diffusion coefficient that controls the shifting magnitude and absorbs the constants of proportionality. The gradient of the particle concentration can be found through an SPH gradient operator:
The proportionality coefficient D is computed through a form proposed by [Skillen et al., 2013]. It is set to be large enough to provide effective particle shifting, while not introducing significant errors or instabilities. This is achieved by performing a Von Neumann stability analysis of the advection-diffusion equation:
where Δtmax is the maximum local time step that is permitted by the CFL condition for a given local velocity and particle spacing. The CFL condition states that:
Combining Eq. 24 and 25 we can derive an equation to find the shifting coefficient D:
where A is a dimensionless constant that is independent of the problem setup and discretization and dt is the current time step. Values in the range of [1,6] are proposed with 2 used as default.
The shifting algorithm is heavily dependent on a full kernel support. However, particles at and adjacent to the free surface cannot obtain the full kernel support, which will introduce errors in the free-surface prediction, potentially causing non-physical instabilities. Applying Fick’s law directly would result in the rapid diffusion of fluid particles from the fluid bulk, due to the large concentration gradients at the free surface.
To counter this effect, [Lind et al., 2012] proposed a free-surface correction that limits diffusion to the surface normal but allow shifting on the tangent to the free surface. Therefore, this correction is only used near the free surface, identified by the value of the particle divergence, which is computed through the following equation, first proposed by [Lee et al., 2008]:
This idea is applied to the DualSPHysics code by multiplying the shifting distance of Equation (22) with a free-surface correction coefficient AFSC.
where AFST is the free-surface threshold and AFSM is the maximum value of the particle divergence. The latter depends on the domain dimensions:
while the free surface threshold is selected for DualSPHysics as:
To identify the position of the particle relative to the free surface, the difference of the particle divergence to AFST is used. Therefore, the full shifting equation (Eq. 22) with the free surface correction is:
More information about the shifting implementation can be found in [Mokos, 2013].
In DualSPHysics two explicit time integration schemes are implemented. For brevity, the governing equations are written as
These equations are integrated in time using a computationally simple Verlet based scheme or a more numerically stable but computationally intensive two-stage Symplectic method.
The Verlet scheme [Verlet, 1967] is commonly used in molecular dynamics since it is a low computational cost scheme with a second order accurate space integrator that does not require multiple calculation steps within an iteration interval, the WCSPH variables are calculated according to
Due to the integration over a staggered time interval, the equations of density and velocity are decoupled, which may lead to divergence of the integrated values. Therefore, an intermediate step is required every N_S steps (Ns≈40 is recommended) according to
where the subscript n denotes the time step and t= nΔt
The symplectic position Verlet time integrator scheme [Leimkuhler and Matthews, 2016] is second order accurate in time. It is ideal for Lagrangian schemes as it is time reversible and symmetric in the absence of diffusive terms that preserve geometric futures. The position Verlet scheme in the absence of dissipation forces reads,
However, in the presence of viscous forces and density evolution in DualSPHysics, the velocity is required in n+1/2 step thus, a velocity Verlet half step is used to compute the required velocity for the acceleration and density evolution for 
where 
Finally, the density evolution follows the half time steps of the symplectic position Verlet scheme as follows [Parshikov et al, 2000],
where 
With explicit time integration schemes the timestep is dependent on the Courant- Friedrichs-Lewy (CFL) condition, the forcing terms and the viscous diffusion term. A variable time step Δt is calculated according to [Monaghan et al., 1999] using
where Δtf is based on the force per unit mass (|fa|), and Δtcv combines the Courant and the viscous time step controls.
In DualSPHysics, the boundary is described by a set of particles that are considered as a separate set to the fluid particles. The software currently provides functionality for solid impermeable and periodic open boundaries. Methods to allow boundary particles to be moved according to fixed forcing functions are also present.
The Dynamic Boundary Condition (DBC) is the default method provided by DualSPHysics [Crespo et al., 2007]. This method sees boundary particles that satisfy the same equations as fluid particles, however they do not move according to the forces exerted on them. Instead, they remain either fixed in position or move according to an imposed/assigned motion function (i.e. moving objects such as gates, wave-makers or floating objects).
When a fluid particle approaches a boundary and the distance between its particles and the fluid particle becomes smaller than twice the smoothing length (h), the density of the affected boundary particles increases, resulting in a pressure increase. In turn this results in a repulsive force being exerted on the fluid particle due to the pressure term in the momentum equation.
Stability of this method relies on the length of time step taken being suitably short in order to handle the highest present velocity of any fluid particles currently interacting with boundary particles and is therefore an important point when considering how the variable time step is calculated.
Note that different boundary conditions have been tested in DualSPHysics in the work of [Domínguez et al., 2015].
DualSPHysics provides support for open boundaries in the form of a periodic boundary condition. This is achieved by allowing particles that are near an open lateral boundary to interact with the fluid particles near the complimentary open lateral boundary on the other side of the domain.
In effect, the compact support kernel of a particle is clipped by the nearest open boundary and the remainder of its clipped support applied at the complimentary open boundary [Gómez-Gesteira et al., 2012a].
Within DualSPHysics it is possible to define a pre-imposed movement for a set of boundary particles. Various predefined movement functions are available as well as the ability to assign a time-dependant input file containing kinematic detail. These boundary particles behave as a DBC described in Section 3.8.1, however rather than being fixed, they move independently of the forces currently acting upon them. This provides the ability to define complex simulation scenarios (i.e. a wave-making paddle) as the boundaries influence the fluid particles appropriately as they move.
It is also possible to derive the movement of an object by considering its interaction with fluid particles and using these forces to drive its motion. This can be achieved by summing the force contributions for an entire body. By assuming that the body is rigid, the net force on each boundary particle is computed according to the sum of the contributions of all surrounding fluid particles according to the designated kernel function and smoothing length. Each boundary particle k therefore experiences a force per unit mass given by
where fka is the force per unit mass exerted by the fluid particle a on the boundary particle k, which is given by
For the motion of the moving body, the basic equations of rigid body dynamics can then be used
where M is the mass of the object, I the moment of inertia, V the velocity, Ω the rotational velocity and R0 the centre of mass. Equations 38.1 and 38.2 are integrated in time in order to predict the values of V and Ω for the beginning of the next time step. Each boundary particle within the body then has a velocity given by
Finally, the boundary particles within the rigid body are moved by integrating Eq. 39 in time. The works of [Monaghan et al., 2003] and [Monaghan, 2005] show that this technique conserves both linear and angular momentum. [Bouscasse et al., 2013] presented successful validations of nonlinear water wave interaction with floating bodies in SPH comparing with experimental data from [Hadzić et al., 2005] that includes deformations in the free-surface due to the presence of floating boxes and the movement of those objects during the experiment (heave, surge and roll displacements). Several validations using DualSPHysics are performed in [Canelas et al., 2015] that analyse the buoyancy-driven motion with solid objects larger than the smallest flow scales and with various densities. They compared SPH numerical results with analytical solutions, with other numerical methods [Fekken, 2004] and with experimental measurements.
The Dynamic Boundary Condition (DBC) presents some drawbacks such as: over dissipation; the evolution of density and pressure of the fixed boundary particles leading to unphysical values for those points, and unphysically large boundary layers. Another drawback of using DBC arises when considering the steady-state 2D channel flows such as Poiseuille and Couette flows, as the boundary layer is not accurately solved unless a high resolution is used. The approximation of no-slip by setting zero velocity at each of the boundary particles fails to capture the condition for low resolutions. This then leads to higher than expected velocities close to the boundary and throughout the flow. Despite all this, DBC can be successfully used for coastal engineering problems due the capability of discretizing complex 3D geometries without the need of implementing complex mirroring techniques or semi-analytical wall boundary conditions.
We will describe here the approach presented in [English et al. 2019] that includes a modification of DBC (mDBC). Within this implementation, the boundary particles are arranged in the same way as the boundary particles in the original DBC, with the boundary interface located half a particle spacing from the inner most layer of boundary particles. For each boundary particle a ghost node is projected into the fluid across the boundary interface in a similar procedure to [Marrone et al. 2011]. For a flat surface the ghost node is mirrored across the boundary interface along the direction of the boundary normal pointing into the fluid (Figure 3-1a) and fluid properties are found at this ghost node through a corrected SPH sum over the surrounding fluid particles only (Figure 3-1b). For a boundary particle located in a corner using the boundary normal would not work due to the presence of more than one normal. However, the boundary interfaces of each solid boundary meet to form an interface corner, and the ghost node is mirrored through the point of this corner into the fluid region (Figure 3-1c). Again, fluid properties are found at the ghost node through a corrected SPH sum of the surrounding fluid (Figure 3-1d).
Figure 3-1. Mirroring of ghost nodes (crosses) and the kernel radius around the ghost nodes for boundary particles in a flat surface (a) and a corner (c). Fluid particles (pink) included in the kernel sum around ghost nodes for boundary particles in a flat surface (b) and a corner (d).
The boundary particles receive fluid properties using the values calculated at the ghost node and an extrapolation method similar to the one used for open boundaries in [Tafuni et al. 2018]. For the density of the boundary particle ρb the density ρg and its gradient [∂xρg;∂yρg;∂zρg] are computed at the ghost node using the first order consistent SPH interpolation proposed by [Liu and Liu, 2006], which requires solving the following linear system for each ghost node g:
where:
where the volume Vj is computed as Vj=mj ⁄ ρj.
To evaluate the matrix Ag and the vector bg at each ghost node g, a new particle interaction loop needs to be calculated. The new interaction requires the inversion of the matrix Ag; however, this requires limited additional computational efforts when compared to the original DBC formulation. Note also that the matrix Ag is diagonally dominant so rarely becomes ill-conditioned, even for a very disordered particle distribution. Due to this property, the adopted formulation to compute ρg and [∂xρg;∂yρg;∂zρg] is still suitable for real engineering applications with fragmented free surfaces. Once the density and density gradient are computed at the ghost node, then the density of the boundary particle ρg is obtained by means of a Taylor series expansion with the values found using the above relations through:
where rg and rg are the position of the boundary particle and associated ghost node respectively. Using this method to find the density of the boundary particle, creates smoother and more physical pressure and density fields in the boundary, and reduces the gap between the boundary and the fluid observed when using only DBC. The velocity at the ghost node is found using a Shepard filter sum
The boundary particle then receives this velocity with the direction reversed to create a no-slip condition at the boundary interface. The no-slip condition created using this method is more accurate than using the zero-velocity approach of the original DBC as the negative velocity is more effective at reducing the velocity of the fluid. With this method it is also possible to create a free slip boundary by giving the boundary particles the exact tangential velocity found at the ghost node.
Wave generation is included in this version of DualSPHysics, which includes piston- and flap-type long-crested wave generation. In this way, the numerical model can be used to simulate a physical wave flume. Both monochromatic (regular) and random waves can be generated. Second-order correction for bound long waves (sub-harmonics) is implemented for random waves with piston-type wave maker only. For regular waves, second-order correction (super-harmonics) is implemented for both piston- and flap-type wavemaker. The following sections refer only to the piston-type wavemaker.
The Biesel transfer functions express the relation between wave amplitude and wavemaker displacement [Biesel and Suquet, 1951], under the assumption of irrotational and incompressible fluid and constant pressure at the free surface. The transfer function links the displacement of the piston-type wavemaker to the water surface elevation, under the hypothesis of monochromatic sinusoidal waves in one dimension in the x-direction:
where H is the wave height, d the water depth, x is distance and δ is the initial phase. The quantity ω=2π/T is the angular frequency and k=2π/L is the wave number with T equal to the wave period and L the wave length. The initial phase δ is given by a random number between 0 and 2π.
Eq. 44 expresses the surface elevation at infinity that Biesel defined as the far-field solution. The Biesel function can be derived for the far-field solution and for a pistontype wavemaker as:
where S0 is the piston stroke. Once the piston stroke is defined, the time series of the piston movement is given by:
The implementation of a second order wavemaker theory will prevent the generation of spurious secondary waves. The second order wave generation theory implemented in DualSPHysics is based on [Madsen, 1971] who developed a simple second-order wavemaker theory to generate long second order Stokes waves that would not change shape as they propagated. The theory proposed by [Madsen, 1971] is straightforward, controllable, computationally inexpensive with efficient results, and is accurate for waves of first and second order.
The piston stroke S0 can be redefined from Eq. 45 as S0=H/m1 where:
Following [Madsen, 1971], to generate a wave of second order, an extra term must be added to Eq. 46. This term, for piston-type wavemaker, is equal to:
Therefore, the piston displacement, for regular waves, is the summation of Eq. 46 and Eq. 48:
Madsen limited the application of this approximate second order wavemaker theory to waves that complied with the condition given by HL^2/d^3 < (8π^2)/3. A specific warning is implemented in DualSPHysics to inform the user whether or not this condition is fulfilled.
Monochromatic waves are not representative of sea states that characterise real wave storm conditions. Sea waves are mostly random or irregular in nature. Irregular wave generation is performed in DualSPHysics based on [Liu and Frigaard, 2001]. Starting from an assigned wave spectra, the Biesel transfer function (Eq. 45) is applied to each component in which the spectrum is discretised. The procedure for the generation of irregular waves is summarised as follows:
- Defining the wave spectrum through its characteristic parameters (peak frequency, spectrum shape, etc.).
- Dividing the spectrum into N parts (N>50) in the interval (fstart, fstop), where generally the values assumed by the spectrum (Sη) at the extremes of this interval are smaller than the value assumed for the peak frequency, fp: Sη(fstart)≤0.01·Sη(fp) and Sη(fstop) ≤ 0.01 · Sη(fp).
- The frequency band width is so-defined as Δf=(fstop-fstart)/N. The irregular wave is decomposed into N linear waves.
- Determining the angular frequency ωi, amplitude ai and initial phase δi (random number between 0 and 2π) of each i-th linear wave. The angular frequency ωi and amplitude ai can therefore be expressed as follows:
- Converting the time series of surface elevation into the time series of piston movement with the help of Biesel transfer function:
- Composing all the i-th components derived from the previous equation into the time series of the piston displacement as:
In DualSPHysics two standard wave spectra have been implemented and used to generate irregular waves: JONSWAP and Pierson-Moskowitz spectra. The characteristic parameters of each spectrum can be assigned by the user together with the value of N (number of parts in which the spectrum is divided).
The user can choose among four different ways to define the angular frequency. It can be determined assuming an equidistant discretization of the wave spectrum (fi=fstart+iΔf-Δf/2), or chosen as unevenly distributed between (i-0.5)Δf and (i+0.5)Δf. An unevenly distributed band width should be preferred: in fact, depending on N, an equidistant splitting can lead to the repetition of the same wave group in the time series that can be easily avoided using an unevenly distributed band width. The last two ways to determine the angular frequency of each component and its band width consist of the application of a stretched or cosine stretched function. The use of a stretched or cosine stretched function has been proven to lead the most accurate results in terms of wave height distribution and groupiness, even when the number of spectrum components N, is relatively low. If there is a certain wave group that is repeating, finally the full range of wave heights and wave periods is not reproduced and statistically the wave train is not representing a real sea state of random waves.
A phase seed is also used and can be changed in DualSPHysics to obtain different random series of δi. Changing the phase seed allows generating different irregular wave time series both with the same significant wave height (Hm0) and peak period (Tp).
When natural waves are reproduced in a laboratory or in a numerical model the wave generator is normally controlled by a first-order signal. The wave set-down (bound long waves) generates drift velocities directed landwards under the crests and seawards under the troughs [Ottesen-Hansen et al., 1980]. However, at the wavemaker, there is no flow through the paddle, therefore the natural drift velocities are compensated by identical ones with opposite signs. The latter velocities create a progressive long wave that is not bound anymore but free. This phenomenon is called parasitic long waves and it results in an exaggeration of the long wave effects. Second-order wave generation is required in order to cancel out this parasitic long wave to achieve the drift velocities that characterise the set-down. Another long wave, called displacement long wave, is caused by finite wavemaker displacements away from the mean position.
The correction of the first-order wave generation at the wavemaker must be capable of compensating both parasitic long waves and displacement long waves. The method implemented in DualSPHysics is based on the solution for the control signal of the wavemaker that is described in [Barthel et al., 1993].
Generation of solitary waves in DualSPHysics is achieved by means moving boundaries. In particular, a piston-type wavemaker is employed to generate single and multiple solitary waves [Dominguez et al., 2019].
3.9.5.1. Single solitary wave
Goring [1978] laid the groundwork for the generation of solitary waves by means of a piston-type wavemaker.His method was based on one fundamental hypothesis: the horizontal velocity of the water particles under the wave crest, u (depth-integrated velocity), is equal to the velocity of the wavemaker.
Equation (54) defines the aforementioned assumption, where xs refers to the piston displacement.
In order to give the possibility to the user to generate the solitary wave profile that would be the most appropriate for his problem and the closest the the real studied phenomenon (e.g. tsunami wave), three different theories for solitary wave generation have been implemented in DualSPHysics. For all cases, a piston-type wavemaker is employed. The three different approaches result in different wave propagation properties (i.e. different loss of amplitude and differences in generated trailing waves) [Dominguez et al., 2019]. One solution is derived from the Rayleigh solitary wave solution; the second one uses a wavemaker law-of-motion derived from the KdV solution; the third one consists of solving equation (57) iteratively using the Newton-Raphson method. The following subsections described the three different theories.
a) Rayleigh approximation to a solitary wave Following Serre (1953), the theoretical free-surface elevation expressed by equation (57) can be re-written as follows:
b) First-order shallow water model (KdV) Clamond and Germain (1999) used a paddle law of motion derived from the KdV (first-order shallow water) solution. The piston displacement can be expressed as follows:
c) Goring solution (Boussinesq) Equations (56) and (57) correspond, respectively, to the equation of the piston displacement and the theoretical water elevation using the solution proposed by Goring which requires an iterative procedure to solve Equation (56).
3.9.5.2. Multiple solitary waves
To generate a series of multiple solitary waves, the law of motion of the wavemaker is calculated for each i-th solitary wave. Any of the three theories indicated in the previous sections can be used. Then, the time lag between the i+1-th and the i-th solitary wave has to be specified as a fraction or multiple of the generation time used for the i-th solitary wave. The time lag-between the generation of the two solitary waves will be expressed by λTf, being Tf the generation time used to generated the first solitary wave. The coefficient λ is a positive real number that multiplies Tf. If λ<1, the generation of the second solitary wave starts when the displacement of the piston has not reached yet the maximum stroke that is needed to generate the first solitary wave.
The generation of multiple solitary waves is a simple but efficient way to model a wave train of several tsunamis which might have been triggered by the same tectonic event.
NewWave theory [Whittaker et al., 2017] is employed to create the time series of focused wave groups. The NewWave theory describes the most probable shape of a large wave in a given sea state. Theoretically, the focus wave group energy reaches the structure in a compact, maximized form if the focus location is close enough to the focus location. A NewWave-type focused wave group comprising N infinitesimal wave components is given by:
where S_ηη is the power spectral density, ωis the angular frequency, t is time, σ is the standard deviation of the sea state (with an associated variance σ^2=∑S_ηη (ω_i )Δω in this discretised form) and ki is the wavenumber of the i-th wave component with angular frequency ωi and related to it by the familiar linear dispersion relation ω^2=gk*tanh(kh) (where g is the acceleration due to gravity and h is the water depth), and x is the horizontal distance. All wave components come into phase at the focus location x_f and focus time tf to form a large wave with a linear focus amplitude equal to A.
A full range of focusing behaviours can be allowed by introducing the the phase angle ϕ of the group at focus (e.g. crest, trough, otherwise), while the energy concentration within the group is independent of the value of ϕ. However, the wave shape can affect its breaking patterns and therefore the impacts exerted on the structure.
The use of wave absorption allows generating long time series of sea waves in relatively short domains with negligible wave reflection. Wave absorption techniques can be categorized as passive and active, respectively. Passive wave absorbers are usually required to damp the wave energy and to reduce reflection exerted by the boundary of the model domain. In active absorption systems, the wavemaker real-time displacement is corrected to cancel out the reflected waves and to damp the re-reflection phenomenon. Active wave absorption is implemented in DualSPHysics for a piston-type wavemaker only.
A damping zone is implemented in DualSPHysics as passive absorption system. The implemented damping system consists in gradually reducing the velocity of the particles at each time step according to their location, but using quadratic decay rather than exponential. In this way, the velocity is modified following
where v0 is the initial velocity of the particle i, v is the final velocity and f(x,Δt) is the reduction function defined as
where Δt is the duration of the last time step, x is position of the particles, x0 and x1 are the initial and final position of the damping zone, respectively. It is recommended to use one wavelength, L, as the length of the damping zone. The coefficient β modifies the reduction function that is applied to the velocity.
The active wave absorption system implemented in DualSPHysics is based on the time-domain filtering technique that uses the free-surface elevation at the wavemaker position as feedback for the control of the wavemaker displacement. The target wavemaker position e(t) is corrected in real time in order to avoid reflection at the wavemaker. It is necessary to estimate the free-surface elevation of the reflected waves, ηR, to be absorbed, by comparing the target incident free-surface elevation, ηI, with the measured one in front of the wavemaker, ηSPH. This is measured at 4-10·hSPH from the wavemaker. The reflected free-surface elevation can be expressed as
The wavemaker velocity is then modified to match the velocity induced by the wave that will be absorbed. For a piston-type wavemaker the wave absorption is performed using linear long wave theory [Schaffer and Klopman, 2000; Dean and Dalrymple, 1991]. The velocity correction to absorb the reflected waves, UR, is expressed as
where ηR is the free-surface elevation of the reflected waves and g is the gravitational acceleration. The corrected wavemaker velocity UC is the subtraction of UR from the theoretical wave maker velocity UI. This UI is the derivative in time of the wavemaker displacement e(t). Here the implementation details are expressed for the regular wave case. The theoretical velocity at time t, UI(t), can be computed as:
The free-surface elevation in front of the wavemaker is measured, the target incident free-surface elevation is calculated and therefore the velocity correction UR is estimated. The corrected wavemaker velocity at the instant t+dt can be expressed as
The wavemaker position at t+dt is then corrected using the following expression
Small deviations of the mean water level from the zero can imply accumulations in time leading to a drift of the wavemaker from its initial position that will grow indefinitely. In order to prevent this drift, the wavemaker should be forced back to its zero position. Hence, a drift correction algorithm is implemented in the code. The algorithm checks when the 80% of the maximum forward or backward stroke is reached (80% was assumed as default value but it can be modified by the user). If this happens the wave board is pushed/pulled slowly back, while continuing to generate the target waves, in such a way that its average position will finally correspond to its initial zero-position (i.e. at the beginning of the simulation). A smoothed transition, in form of a power function, is used to prevent abrupt changes in the wavemaker displacement.
The discrete element method (DEM) allows for the computation of rigid particle dynamics, by considering contact laws to account for interaction forces. The coupled numerical solution, based on SPH and DEM discretisations, resolves solid-solid and solid-fluid interactions over a broad range of scales.
Forces arise whenever a particle of a solid object interacts with another. In the particular case of a solid-solid collision, the contact force is decomposed into Fn and Ft, normal and tangential components respectively. Both of these forces include viscous dissipation effects. This is because two colliding bodies undergo a deformation which will be somewhere between perfectly inelastic and perfectly elastic, usually quantified by the normal restitution coefficient
The total forces are decomposed into a repulsion force, Fr, arising from the elastic deformation of the material, and a damping force, Fd, for the viscous dissipation of energy during the deformation.
Figure 3-2 generally illustrates the proposed viscoelastic DEM mechanism between two interacting particles.
Figure 3-2. Schematic interaction between particles with viscoelastic DEM mechanism.
The normal force is given by
where the stiffness is given by
and the damping coefficient is
where 
The restitution coefficient ij e is taken as the average of the two materials coefficients, in what is the only calibration parameter of the model. It is not a physical parameter, but since the current method does not account for internal deformations and other energy losses during contact, the user is given a choice to change this parameter freely in order to control the dissipation of each contact. The reduced radius and reduced elasticity are given by
where Ri is simply the particle radius, Ei is the Young modulus and νp is the Poisson coefficient of material i, as specified in the Floating_Materials.xml.
This results in another restriction to the time-step, adding
to the existing CFL restrictions (Eq. 35), where vn is the normal relative velocity and M* is the reduced mass of the system where there is a contact.
Regarding tangential contacts, friction is modelled using the same model:
where the stiffness and damping constants are derived to be
as to insure internal consistency of the time scales between normal and tangential components. This mechanism models the static and dynamic friction mechanisms by a penalty method. The body does not statically stick at the point of contact, but is constrained by the spring-damper system. This force must be bounded above by the Coulomb friction law, modified with a sigmoidal function in order to make it continuous around the origin regarding the tangential velocity:
where μIJ is the friction coefficient at the contact of object I and object J and is simply taken as the average of the two friction coefficients of the distinct materials, indicated in the Floating_Materials.xml.
More information about DEM implementation can be found in [Canelas, 2015; Canelas et al., 2016].
Devices composed of rigid bodies interacting through frictional contacts and several nonlinear constraints are extensively used in many engineering fields, either featuring a small number of unilateral contacts or including thousands of contacts between a large number of parts. Mechanisms involving contacts and impacts between parts can be modeled in terms of multi-body systems with unilateral constraints. The simulation of rigid contacts entails the solution of non-smooth equations of motion: the dynamics are non-smooth since the non-interpenetration, collision, and adhesion constraints are discontinuous. Considering the success of SPH for fluid descriptions and non-smooth multibody solvers for mechanical systems, attempting to couple both under a generalized framework should provide new simulation possibilities, by leveraging the strengths in both methods.
The multi-body solver used is Project Chrono, an open-source high-performance library and set of applications enabling many applications such as wheeled and tracked vehicles operating on deformable terrains, robots, mechatronic systems, compliant mechanisms, and fluid solid interaction phenomena. Systems can be made of rigid and flexible/compliant parts with constraints, motors and contacts; parts can have three-dimensional shapes for collision detection.
DualSPHysics at the moment uses the rigid body, constraints and collision detection parts of the library. This allows a user to use a set of bodies described in meshes, specify restrictions applied to the bodies (from the implemented list of joints, hinges and springs) and compute interactions between bodies, in a similar fashion to the DEM implementation, but in a more stable manner. Realistic frictional behavior is also now possible, since a full Coulomb sliding/sticking/rolling model is used.
The implementation strategy is to couple both models with a message passing interface. Once quantities from the rigid bodies are computed by DualSPHysics, the time step, along with the linear and angular accelerations from each body, computed solely from fluid- body interaction, are sent to our Project Chrono module. For that time step, Chrono returns the linear and angular velocities, as well as CC position, computed by integrating the fluid contributions with the dynamic or kinematic restrictions of the system, including collisions.
Figure 3-3. Flow chart of the coupling between DualSPHysics and Project Chrono.
Section 3.8 lays the fundamental ideas behind the modelling of the fluid phase and the modelling of otherwise unconstrained floating rigid bodies. Mechanisms to cope with interactions between rigid bodies and other constraints need to be explored. The ideas laid by Potapov-2001 and later adapted by Canelas-al-2016 represent traditional regularization strategies, which model contacts and other restrictions by means of spring-dashpot systems (DEM class methods). From a SPH implementation perspective, this approach requires small code adaptations and allows the adoption of already used integrators. The high stiffness affecting the explicit integration, however, typically imposes prohibitively small time steps to ensure stability and guarantees the need to tune parameters Kruggel-Emden-2007. These issues motivate the search for methods that can deal with multiple frictional restrictions, even in cases with thousands or millions of moving bodies. To that end, much attention was drawn by time-stepping approaches that produce weak solutions of the DVI that describes the continuous time motion of rigid bodies with collision, contact, and friction, as those applied in Project Chrono. The DVI as a problem formulation was recently introduced in full generality and classified by differential index Pang-2008, though earlier numerical approaches based on DVI formulations do exist. On the following subsections, a DVI description is built by defining the notation space, writing constraints that represent physical characteristics of the system and finally writing the full equation system.
At a time , the position of the system is described by generalized coordinates
(which may include rotational coordinates that cannot be defined over a subspace homeomorphic to
, for some
), and generalized velocities
. In classical mechanics,
is continuous, and we can write
, where
is used to transform the generalized velocities into derivatives of the generalized positions. For instance, when dealing with rotations,
can be a linear mapping from three dimensional angular speeds into four dimensional time derivatives of unit quaternions. Four-dimensional unitary quaternions are adopted,
, though their space is not homeomorphic to
. This allows the position vector as
and the velocity vector as
. Therefore a system with
bodies in three dimensions is represented by
position and
speed coordinates.
Two rigid bodies should not penetrate, and, if they are in contact, there should be friction acting at the interface. To enforce the non-penetration constraint, we assume that there exists a function which we call the gap function, that satisfies
(82)
For such a function, the non-penetration constraint becomes . An example of such a mapping is the signed distance function, which is differentiable when the bodies are smooth and convex, at least up to some value of the interpenetration. For most cases, even simple ones involving the relative position of two spheres, a differentiable signed distance function cannot be defined for all values of
. The fact that
can be differentially defined only on a neighborhood of the set
can be accommodated at the cost of making the analysis substantially more involved, prompting the following assumption: any contact is described by a gap function
that is everywhere twice continuously differentiable.
The model we represent and approximate is the Coulomb friction model. If a position is feasible and the contact is active, that is,
, then at the contact we have a normal force and a tangential force. Let
be the normal at the contact, pointing toward the exterior of the body, and let
and
be the tangents at the contact such that
,
and
are mutually orthogonal vectors of unitary length and a function of the position
. Let
to refer to velocities, and the subscripts 1 and 2 to refer to quantities related to the two linearly independent tangential directions at a given contact. The reaction force is impressed on the system by means of multipliers
,
and
. The normal component of the force is
and the tangential component of the force is
. The Coulomb model consists of the following constraints:
(83)
where is the relative tangential velocity at contact. The first part of the constraint can be restated as
(84)
where is a cone in three dimensions, whose slope is
. The constraint
requires that the tangential force be opposite to the tangential velocity. This results in the reaction force being dissipative. In fact, an equivalent convenient way of expressing this constraint is by using the maximum dissipation principle
(85)
These constraints are represented by mapping vectors ,
and
from contact coordinates to generalized coordinates. For example, if we have a two-body system, then the generalized coordinates in the three-dimensional space are embedded in a twelve-dimensional space by using the coordinates
. For a three-dimensional vector
, the mapping to generalized coordinates is given by
(86)
where and
are the relative positions of the contact point with respect to the centers of mass of the two bodies. Using this mapping, we denote the generalized vector version of
,
and
by
. One unfortunate side effect of generalized coordinates mapping is that, in the new coordinates,
cease to be mutually orthogonal.
If u is the generalized velocity, the tangential velocity satisfies the following
(87)
In generalized coordinates, the Coulomb model can now be written as
(88)
The other dynamical data needed for the model are the mass matrix and the external force
. The mass matrix
is positive definite and constant.
Assuming now that we have p potential contact constraints, which are enforced by the non-penetration constraints
. The superscript i denotes the data associated to the potential contact i. The continuous model is the following differential variational inequality
(89)
System 89 must be discretized for use in a numerical model. Casting the system in Cone Complementary Problem (CCP) form and solving it with a fixed point iterative method allows for the expected performance levels, while taking advantage of the DVI capabilities in describing generalized restrictions. The full discretization can be seen in Project Chrono's reference papers, Anitescu (2010) and Tasora (2011).
Multi-scaled complex simulations, often characterised by three-dimensional domains and long duration, require huge computational resources, if only one stand-alone model is employed. Despite the significant improvements in computational time of SPH-based solvers by using GPU (Graphics Processing Unit) and multi-GPU techniques (Domínguez et al., 2013), the computational cost is still high to model realistic engineering problems. In coastal engineering, in particular, it is required modelling properly all processes from wave generation to wave-structure/beach interaction. SPH-based solvers provide a detailed modelling of local phenomena, such as interaction between sea waves and coastal defenses, e.g. Altomare et al. (2015a) and St-Germain et al. (2014). The mesh-free nature of SPH allows mimicking the wave generation system of experimental facilities, employing moving boundaries as piston- or flap-type wavemakers. But, using this schemes, waves need to be generated and to propagate over relative long domains in order to simulate properly wave transformation and to guarantee the correct wave sea state at the toe of the structure, similarly to what happens in physical model tests. To reduce the computational cost of such simulations, less computationally demanding wave propagation models (e.g. Boussinesq models, Non-Linear Shallow Water equation models) can be coupled with SPH-based solvers. Here two coupling schemes between DualSPHysics model and SWASH model are briefly described (Figure 3-4): multi-layered piston wavemaker (Altomare et al., 2015b) and relaxation zone technique (Altomare et al., 2018). Both techniques are 1-way offline coupling schemes. The SWASH model is an open-source non-hydrostatic wave-flow model (available at swash.sourceforge.net). A full description of the numerical model, boundary conditions and numerical scheme and applications are given in (Zijlema et al., 2011). The differences are following described:
- In (Altomare et al., 2015b) the real time displacement of the multi-layered piston wavemaker (consisting of boundary particles) is reconstructed on the basis of the velocity time series calculated by SWASH model in a specific point of the fluid domain. There is no reflection compensation, therefore the use of this technique is recommended for low-reflective cases or very short time series.
- In the relaxation zone technique implemented by (Altomare et al., 2018), the coupling zone is extended from one single location to an area. The orbital velocity calculated by the wave propagation model is imposed to the fluid particles in the SPH model. The particle velocity is controlled by analytical solution with a weighting C function: in this way, the reflected waves are compensated inside the relaxation zone.
Figure 3-4. Representation of the coupling techniques between DualSPHysics and SWASH.
A scheme of the multi-layered approach is depicted in Figure 3-5. For sake of simplicity, the geometry of the wave flume at Flanders Hydraulics Research (Belgium) is resembled as example. The top picture represents the whole flume domain and indicates the points where the coupling is set. The middle picture represents the DualSPHysics domain and multi-layered piston wavemaker consisting of SPH boundary particles. The bottom picture is the SWASH domain: the real flume geometry is replaced after the coupling point by horizontal bottom and absorbing boundary conditions. In this way, only incident wave characteristics will be measured at the coupling location. In the small box on the left-hand side of Figure 3-5, a scheme of the multi-layered boundary is depicted: the different layers indicate the number of layers used in SWASH simulation (8 is recommended). The velocities provided by SWASH at the coupling point for each layer are then translated in the boundary particle displacement and are passed to DualSPHysics. The final displacement for each of the boundary particles in SPH is calculated based on the interpolated and velocity between neighbor layers. The method, validated for 2DV or quasi-3D applications, has been implemented to work fully in 3D, by reading the velocity in x- and y- directions and integrated them over time to obtain the three-dimensional displacement of the moving boundaries.
Figure 3-5. Scheme of the multi-layered approach.
Relaxation zone technique is a general framework to generate waves and currents in DualSPHysics, which can be employed as stand-alone wave generation and wave absorption method, as alternative moving boundaries, or as a proper coupling scheme between DualSPHysics and a wave propagation model. The basics of the method remain the same in both cases. We refer here to implementation in 2DV or quasi-3D. Horizontal and, if necessary, vertical velocities are imposed to SPH fluid particles within a specific generation. The particle velocity distribution between the center of the relaxation zone and its extremes is controlled by with a weighting C function, the use of which prevents potential particle disorder due to the discontinuity of particle velocities at the interface between the RZ and main computational region. The velocities in the RZ can be expressed as follows:
where 
Figure 3-6. Shape of the weighting function C for different values of ψ and β.
The value assumed by the function depends on the set of values chosen for the 𝜓 and β parameters and on the distance from the center of the RZ, being zero at both extremes. The variation of the two parameters allows exploring the performance of the RZ depending on the shape of C. The basic scheme of RZ is depicted in Figure 3-7, where the RZ is defined by the area bounded by the red dotted line. At the right side of the RZ it is located the main fluid domain. Since the generated waves will freely propagate in both directions along the x-axis, a damping area is placed at the left side of the RZ (bounded by a blue dotted line). It has been verified that a damping area 0.25 times the characteristic wavelength is already quite effective.
Figure 3-7. Generic scheme of relaxation zone in DualSPHysics model.
For wave generation, non-linear effects might induce a drift in the particle displacement that can be identified actually by the mass transport predicted by Stokes at second-order solution for monochromatic waves. The drift velocity can be expressed by the following equation:
This mass transport leads to water level increase in the main fluid domain and consequent water level drop in the RZ. In order to avoid unrealistic increase of water level outside the generation zone, a correction to prevent this drift has been implemented. This correction is proportional to the velocity as predicted by Stokes. The performance of the relaxation zone technique to generate waves and to compensate wave reflection is is described in detailed in Altomare et al. (2018), where a step-by-step procedure to design properly the RZ is detailed.
MoorDyn is an open-source dynamic mooring line model that has been designed for coupling with other numerical models (http://www.matt-hall.ca/moordyn.html).
MoorDyn discretizes mooring lines as point masses (nodes i) connected by linear spring-damper segments to provide elasticity in the axial direction. If ri and ri+1 represent the absolute position vectors of two adjacent nodes, the strain in the segment connecting them (i+1/2) is calculated as:
where l is the segment unstretched length.
The tension forces acting within each segment due to material stiffness and internal damping are then calculated, respectively, as:
where diam is line diameter, E is elasticity modulus and Cint is the internal damping coefficient. For chains, which have very small internal structural damping, a small amount of structural damping is still included in the model to critically damp non-physical resonances that can occur at the individual segment natural frequencies due to the model discretization. Bending and torsional stiffnesses are neglected. Hydrodynamic damping and added mass are represented using the Morison equation applied to each node. The tangent direction at a node i is approximated as a unit vector aligned between the nodes on either side:
From the Morison equation, the drag force in the transverse (Dn) and tangential (Dt) directions applied to node i are:
where Cdn and Cdt are transverse and tangential drag coefficients, respectively. These forces are calculated assuming quiescent water; wave kinematics are neglected in MoorDyn’s calculation of hydrodynamic loads.
Added mass on each node (ai) is calculated as a matrix
where Can and Cat are transverse and tangential added mass coefficients, respectively, and Im is the identity matrix.
Vertical seabed contact forces (Bi) are modelled by a spring-damper approach using a stiffness coefficient kb and damping coefficient cb with magnitude (Hall, 2017):
that is activated whenever the elevation of a node, zi, passes below the defined seabed depth, zb. The model is only active when a node contacts the seabed (i.e. when zi ≤ zb).
The total equation of motion for each node along a mooring line, accounting for the terms already mentioned along with node mass (matrix mnode,i), and submerged weight (Wsub,i), is then
This three-by-three matrix equation is solved for each node within a constant-time-step second-order Runge-Kutta integrator.
With its simple formulation, MoorDyn has shown to be computationally efficient and reliable for common offshore renewable energy mooring scenarios, typically giving accurate results with discretizations in the order of 20 segments per mooring line (Hall and Goupee, 2015; Vissio et al., 2015; Sirnivas et al., 2016). Finer discretizations are usually only necessary when transients emerge in the seabed contact forces due to individual nodes touching down abruptly.
Information on using MoorDyn is available in the MoorDyn Users Guide.
Fundamentally, coupling MoorDyn with other codes is based on a loose coupling approach in which fairlead kinematics are passed to MoorDyn, the mooring system dynamics are calculated for one or more time steps, and then the resulting fairlead tension vectors are transferred back to the other simulation code used for the coupling.
The coupling between DualSPHysics and the MoorDyn library can split into three steps:
1st) the motions and rotations initially solved in DualSPHysics (V, Ω, R0) are passed to MoorDyn and used as input for the mooring line fairlead kinematics;
2nd) MoorDyn solves the mooring line behaviour during the time step used in the DualSPHysics model, namely MoorDyn updates the position of the segments of the mooring line and computes the forces at the fairlead connections (dV/dt, dΩ/dt), which are transferred back to DualSPHysics;
3rd) the extra forces calculated in the second step are added in DualSPHysics to obtain the final resulting force acting onto the floating structure that is used to compute the final motions and rotations of the floating structure.
The computational cost of MoorDyn at each time step is negligible compared to the execution time needed to solve one time step in the SPH model.
Figure 3-8. Flow chart of the coupling between DualSPHysics and MoorDyn.
The work of Domínguez et al., 2019 included more details about this coupling and a complete validation comparing numerical values with experimental data. Overall, the results demonstrate the accuracy of the coupling between DualSPHysics and MoorDyn to simulate the motion of a moored floating structure under the action of regular waves in terms of free-surface elevation, motions of the floater and mooring tensions.
Going forward, this modelling approach can be employed to simulate more complex floating structures such as floating wind turbines, buoys, WECs, offshore platforms, etc.
Several types of boundary conditions are available in DualSPHysics to simulate a variety of engineering scenarios. Some examples are the Dynamic Boundary Particles (Crespo et al. 2015), SPH-DEM coupling (Canelas et al. 2016), and Periodic Open Boundary Conditions (Gomez-Gesteira et al. 2012). Nevertheless, none of these formulations is appropriate when the computational problem requires specific flow conditions to be enforced at the domain boundaries. The classic example of flow past an object is one of such cases, where usually an inflow velocity needs to be prescribed at the inlet while other velocity or pressure conditions can be either prescribed or extracted from the fluid domain at the outlet. The treatment of open boundary conditions in SPH is an active area of study among SPH researchers and falls within one of the SPHERIC Grand Challenges. The work of several authors in the literature provides a number of choices to model open boundaries with both WCSPH and ISPH, see for example Ferrand et al. 2017, Federico et al. 2012, Lastiwka et al. 2009, and Vacondio et al. 2012.
To address this issue in DualSPHysics, an open boundary algorithm has been implemented based on the work presented in Tafuni et al. 2016, Tafuni et al. 2017 and Tafuni et al. 2018. Figure 3-9 summarizes the working principles of the algorithm in the generic case of a fluid flowing near a buffer area identifying an open boundary. The innermost dashed curve represents the buffer threshold boundary, i.e. the fluid-buffer interface, followed by a layer of SPH particles used to enforce certain boundary conditions. The buffer width is chosen to equal or exceed the kernel radius so to ensure full kernel support for the fluid particles in the near proximity of an inlet or outlet.
Figure 3-9. Sketch of the implemented open boundary model adapted from Tafuni et al. 2018.
Two ways of providing the information to an open boundary are considered: physical quantities are either assigned a priori or passed from the fluid domain to the buffer zones (inflow and outflow) using ghost nodes. A similar idea was used by Marrone et al. 2013 to enforce closed boundary conditions. The position of the ghost nodes is obtained by mirroring the normal distance of the boundary particles from the permeable surface into the fluid. In order to calculate fluid quantities at the ghost nodes, a standard particle interpolation would not be consistent due to the proximity of these points to an open boundary, which translates into the lack of a full kernel support. The method proposed by Liu and Liu 2006 is thereby adopted to retrieve first order kernel and particle consistency. The multi-dimensional first-order Taylor series approximation of the field function f(x) multiplied by the kernel function, Wk(x), and its first order derivatives, Wk,B(x), is given by
Beta is an index going from 1 to d, the total number of dimensions. Equations (101, 102) form a system of d + 1 equations in d + 1 unknowns, i.e. fk and fk,B. Using particle notation, the approximations of fk and fk,B translate into matrix equations that give the values of fk and fk,B, calculated as:
The above are finally used to calculate the value of fo at the open boundary (more details are shown in Tafuni et al. 2018.)
Figures 3-10 – 3-12 further elucidate the working principles of the open boundary algorithm. To define the reference permeable boundary for the placement of buffer particles and ghost nodes, a series of fixed points is created along a user-defined curve in 2-D or surface in 3-D. Buffer particles are created along the normal direction to the permeable boundary, starting from the fixed points and only for depths below the defined water level. The latter can be either user-defined or extrapolated from the ghost nodes using a procedure similar to that of the shifting algorithm presented in Lind et al. 2012. Buffer particles adjacent to the fluid domain are generated at half the particle spacing dp from the fixed points, while other are positioned multiples of dp away.
Figure 3-10. Sketch of fixed points, normal vector, and buffer threshold.
Figure 3-11. Sketch of the buffer-fluid transition and the creation of a new buffer particle. Here one of the buffer particles has become a fluid particle (light blue) and a new particle is thus created in the buffer (light yellow).
Figure 3-12. Sketch of the change in liquid depth in the buffer area. As depth is increased by the user or from knowledge of the free-surface elevation inside the fluid domain, new particles (light yellow) are created to fill the buffer area.
As shown on the bottom left of Figure 3-11, when a buffer particle crosses the permeable boundary it becomes a fluid particle and simultaneously a new boundary particle is initialized in the buffer. The process is similar in reverse, that is, when a fluid particle crosses a permeable boundary it becomes a buffer particle and follows the flow conditions specified in the buffer. In this case there is no creation of new particles following this transition. Finally, when a buffer particle crosses the domain edge it is discarded from the computational space, with its properties frozen for the remainder of the simulation. Changes in liquid depth are tackled by the algorithm as follows: the maximum height of the buffer particles is determined by the highest fixed point, while the water level can be either imposed by the user or extrapolated from the ghost nodes. In either cases, only buffer particles below the water level are activated and participate in the fluid-buffer and buffer-fluid interactions. In addition to the more general use of the buffer particles, the new algorithm allows:
- Variable Velocity and Pressure Profiles: Unsteady velocity and pressure profiles and/or pressure and velocity gradients along a given direction can now be assigned in buffer areas.
- Variable Water Depth: The water depth can be either predefined or extrapolated from the fluid.
- Backflow: This last option is related to the previous two and allows flow reversion in a buffer.
Liquid-gas flows with mixing and violent free-surface hydrodynamic interaction exist in various industrial and research fields such as coastal and nuclear engineering. They include a diverse range of problems such as overturning or breaking waves and multi-phase pipe flow and are an ideal application of SPH. To address these problems, DualSPHysics includes the multi-phase model presented in [Mokos et al., 2016]. This guide only provides a concise depiction of the multi-phase liquid-gas model implemented in the code.
The SPH formulation is based on the work of [Colagrossi and Landrini, 2003] and [Nugent and Posh, 2000] and its implementation on a GPU is relatively straightforward [Mokos et al., 2015]. The multi-phase model uses a modified version of Tait’s equation of state [Batchelor, 1974] (the same equation used in single-phase DualSPHysics) for incompressible and inviscid fluids:
where γ is the isentropic expansion factor, ρ0 is the initial density of the fluid, cs is the speed of sound, X signifies a constant background pressure, while the last term represents the cohesion forces between the particles of a single phase. As proposed by [Nugent and Posh, 2000] and [Mokos et al., 2015], the coefficient is based on the properties of the different phases and a characteristic length scale of the problem, L
where ρw and ρa are the initial densities of the two phases (water and air in this case).
The variationally consistent forms of the velocity divergence and pressure gradient used in single-phase DualSPHysics are used to update the continuity and momentum equations for the liquid phase. For the gas phase, the same SPH formulation is not applicable due to the large density discontinuity at the interface [Colagrossi and Landrini, 2003]. An extra term for the cohesion forces in the momentum equation for the lighter phase is required [Nugent and Posh, 2000]
When using high resolutions, made possible by the use of GPUs, a numerical inconsistency has emerged for gas-liquid simulations. In violent flows unphysical voids and phase separation occur ultimately leading to numerical instability. New Fickian-based particle shifting algorithms with a selectively activated free-surface correction have been developed to prevent the creation of unnatural voids and maintain numerical stability through nearly uniform distributions [Mokos et al., 2016].
The Fickian shifting algorithm is a modified version of the one developed by [Lind et al., 2012] and described in Section 3.6. The particle concentration C is estimated here using the sum of the smoothing kernel.
The gradient of the particle concentration can be found using a standard SPH gradient approximation through the smoothing kernel [Monaghan, 2005; Violeau and Rogers, 2016]
To evaluate the diffusion coefficient, the approach proposed by [Skillen et al., 2013] is used, shown in Eq. 26. A different approach is used for shifting near the free surface. The concentration gradient near the surface is controlled using the local tangent and normal vectors at the free surface [Lind et al., 2012]
where s and n are the tangent and normal vectors to the free surface in 2D, while βn is a reference concentration gradient (usually taken as the initial value). The parameter αn limits the diffusion in the direction normal to the free surface; for violent flows, such as the ones simulated in this article, αn is set equal to 0. However, for long slow flows, such as standing gravity waves, errors can potentially accumulate at the free surface due to the incompleteness of the kernel [Lind et al., 2012]. A small degree of diffusion at the normal direction is then allowed with the parameter αn set equal to 0.1.
For 3-D cases, the surface correction term needs to be slightly modified as Eq. 110 is only suitable for a two-dimensional simulation. For a 3-D case, the bi-tangent vector, b, also needs to be taken into consideration which allows shifting only on the tangent, s, and bi-tangent directions, treating the normal direction in the same way as the two-dimensional implementation
where b is a unit vector orthogonal to s and n. In both 2-D and 3-D, the free surface is identified using the method of [Lee et al., 2008], Eq. 27.
The free surface term is only applied in the liquid phase. It was found to have a significant effect on the evolution of the flow, restricting the flow expansion on the free-surface. This is an undesirable effect for the gas phase, which should be expanding freely, so the term is not applied there.
For slower flows, in which the Eötvös number has a value close or smaller than 1, a surface tension formulation has been included. The method uses the Continuum Surface Force (CSF) model [Brackbill et al., 1992] and the SPH implementation by [Hu and Adams, 2006]
The Continuum Surface Force (CSF) model uses the curvature of the surface to compute the surface force per unit interfacial area FS
where σ is the fluid surface tension coefficient of phase k, κ is the curvature, and n ̂ is the unit normal to the surface. Due to the difficulty in computing the curvature, [Lafaurie et al., 1994] proposed expressing the surface force as a tensor Π
The model was further simplified by [Wu et al., 1998] who suggested using a colour function instead of calculating the surface normal. The colour function is defined as following for a particle belonging to phase k and does not change throughout the computation
The colour function SPH gradient [Hu and Adams, 2006] exists only if there are neighbouring particles of phase l where l is a different phase from k
Using the colour gradient, the surface tension tensor can be computed and through that, the surface tension force, used in the momentum equation as an external force
Please read “DualSPHysics_v4.0_LiquidGas_GUIDE.pdf” (in doc\guides) for more information and details of the code.
The rheology of materials such as polymers, slurries, pastes and suspensions are a frequent occurrence in industrial applications and academic research. These materials are usually (but not necessarily) characterized by a yield stress below which the material has not yielded and does not flow. Further, the shear strain to shear stress relationship may not be linear as with Newtonian fluids. In general, fluids which do not obey the Newton law of viscosity may be classified as non-Newtonian fluids with viscoelastic, time dependant viscosity and non-Newtonian viscosity.
In DualSPHysics the generalized Herschel-Bulkley-Papanastasiou (HBP) model has been implemented following the work of [Fourtakas and Rogers 2016] on liquid-sediment two-phase flows. The Herschel-Bulkley-Papanastasiou model is capable of modelling a wide range of viscoplastic materials such as Bingham plastic, Bingham pseudoplastic and dilatant fluids. Further, the HBP generalised non-Newtonian model can be used to model shear thinning or thickening materials in the absence of yield strength.
The constitutive equation for the HBP model is written as,
where is the shear stress tensor, derivative of
is the symmetric shear rate tensor given by,
with and
the velocity and transpose respectively and ηapp is the apparent viscosity,
In the above τy is the yield stress of the fluid, μ is the viscosity (or consistency index usually denoted by K), n is the power law index (commonly referred to as the Herschel-Bulkley parameter) and m the exponential shear rate growth parameter (or Papanastasiou parameter) and the magnitude of the symmetric strain rate tensor is defined as,
The Papanastasiou parameter controls the growth of the stress with a finite value for small shear rates in the unyielded region and linear behaviour in the yielded region. Figure 3-13(a) shows the initial rapid growth of stress by varying m whereas Figure 3-13(b) shows the effect of the power law index n.
Note that as m → ∞ the HBP model reduces to the original Herschel-Bulkley model and when n = 1 the model reduces to a simple Bingham model. Consequently, when n = 1 and m = 0 the model reduces to a Newtonian fluid.
Note that, the HBP model does not exhibit a discontinuity at zero shear rates in contrast with a pure Bingham model which is usually modelled using a bi-viscosity model.
Figure 3-13. Initial rapid growth of stress by varying m and effect of the power law index n for the HBP model.
The implementation of the HBP model in DualSPHysics is outlined next. The momentum equation of eq. 11 is written as
where the brackets denote an SPH approximation. Two formulations are available in the multi-phase non-Newtonian formulation for the discretization of the viscous forces. The first uses the formulation of [Morris et al., 1997] that reads,
where vapp is the apparent kinematic viscosity. The second formulation takes advantage of the conservative gradient operator in SPH form such as the viscous forces are discretised by,
The continuity equation (13) has the additional terms of the density ratios of the interpolating and neighbouring particles, allowing different densities to be used in each phase of the system and is written as,
whereas the later takes the following form,
In this implementation, any combination of viscous operators and gradients is allowed in the multi-phase non-Newtonian DualSPHysics formulation. More information about the multi-phase non-Newtonian formulation can be found in [Fourtakas and Rogers 2016].
Future, implementations include a sediment phase with yielded and unyielded sediment regions as described in the aforementioned manuscript.
The coupled interactions between free-surface flows and flexible structures play an important role in a broad range of physical processes in both nature and industry. Therefore, DualSPHysics provides in-built capability for modelling flexible fluid-structure interactions. This capability is based on a unified framework, whereby both the fluid and structure phases are solved via SPH. Furthermore, well-known deficiencies with SPH-based structural modelling (linear inconsistency, tensile instability and rank deficiency) are addressed by adopting a Total Lagrangian formulation with kernel correction and zero-energy mode suppression. The fluid-structure coupling is handled via manipulation of the existing dynamic boundary condition within DualSPHysics, which does not require any geometrical knowledge of the interface. This, together with the unified SPH framework, permits straightforward implementation within DualSPHysics and retains its parallel scalability. A brief description of the approach is given below. However, for further details please see [O'Connor and Rogers, 2021].
Using the Cauchy momentum equation, as used in the fluid phase, to model the structural dynamics leads to serious issues with tensile instability. A solution to this is to reformulate the governing equations in a Total Lagrangian framework, whereby all relevant quantities and operators are measured with respect to an initial (material) reference frame [Belytschko et al., 2000]. The momentum conservation law in the material frame is given by:
(129)
where the zero subscript denotes a quantity/operator measured with respect to the material reference frame, J is the Jacobian determinant of the deformation gradient, and P is the first Piola-Kirchhoff (PK1) stress tensor. The momentum equation is given in discrete form by:
(130)
where L is a corrective (renormalisation) matrix that ensures the constant gradient of a linear field is recovered [Randles and Libersky, 1996], and fHG is a penalty force used to control the spurious zero-energy modes that arise from the rank deficiency within the SPH approximation [Ganzenmuller, 2015]. Note that the kernel derivative and correction matrix only depend on quantities measured in the initial reference frame. Therefore, they only need to be computed once in the initialisation step and can then be stored for the remainder of the simulation.
The PK1 stress is a function of the deformation gradient tensor, which is given by:
(131)
where x denotes the current configuration and X denotes the initial undeformed configuration (material coordinates). In discretised form, the deformation gradient is given by:
(132)
The relevant strain measure in the Total Lagrangian frame is the Green-Lagrange strain, given by:
(133)
where I is the identity matrix. For a hyperelastic material, the Saint Venant-Kirchhoff constitutive model relates the Green-Lagrange strain to the second Piola-Kirchhoff stress tensor via:
(134)
where λ and μ are the Lamé coefficients, given by:
(135)
where E is the Young's modulus and ν is the Poisson ratio. Finally, the PK1 stress is obtained via:
(136)
The standard dynamic boundary condition [Crespo et al., 2007] is used to couple the fluid and structural phases. In this approach, when a fluid particle interacts with a structure particle it sees the structure particle as a regular dynamic boundary particle, but with an arbitrary velocity. Likewise, when a structure particle is interacting with a fluid particle, the interaction is nearly exactly the same as if it were a dynamic boundary particle interacting with a fluid particle. The only difference is now the force from the fluid particle contributes to the acceleration of the structure particle, which is then used to update its position and velocity. Considering this, the acceleration of a fluid particle is thus given by the contributions from all neighbouring fluid, structure and boundary particles. Neglecting body forces for clarity, the acceleration of fluid particle a is:
(137)
where bf, bs and bb are sets containing the neighbouring fluid, structure and boundary particles for particle a. Similarly, the acceleration of structure particle a is:
(138)
As with the fluid phase, two different approaches to time integration are adopted for the structure phase. When using the Symplectic predictor-corrector scheme for the fluid phase, this scheme is also applied to the structure phase. However, the Verlet scheme used for the fluid phase was found to be unstable when applied to the structure phase. Therefore, when using the Verlet scheme for the fluid phase, a one-step semi-implicit Euler scheme is applied to the structure phase:
(139)
Finally, the structure phase also introduces an additional restriction on the timestep. The new timestep criterion becomes:
(140)
where:
(141)