# Smoothed particle hydrodynamics for free-surface and multiphase flows: a review

David Le Touzé, Andrea Colagrossi

► To cite this version:

David Le Touzé, Andrea Colagrossi. Smoothed particle hydrodynamics for free-surface and multiphase flows: a review. Reports on Progress in Physics, 2025, 88 (3), pp.037001. 10.1088/1361-6633/ada80f.hal-05240706

HAL Id: hal-05240706

https://hal.science/hal-05240706v1

Submitted on 4 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
REVIEW • OPEN ACCESS

# Smoothed particle hydrodynamics for free-surface and multiphase flows: a review

To cite this article: David Le Touzé and Andrea Colagrossi 2025 *Rep. Prog. Phys.* **88** 037001

View the [article online](#) for updates and enhancements.

You may also like

- Understanding the behaviour of pulsed laser dry and wet micromachining processes by multi-phase smoothed particle hydrodynamics (SPH) modelling
N Muhammad, B D Rogers and L Li

- Improved kernel gradient free-smoothed particle hydrodynamics and its applications to heat transfer problems
Juan-Mian Lei, Xue-Ying Peng et al.

- Numerical modelling of wire arc additive manufacturing: methods, status, trends, and opportunities
Yun H Kim, Dayalan R Gunasegaram, Paul W Cleary et al.

<figure><img src="image_2.png" /></figure>
# Review

# Smoothed particle hydrodynamics for free-surface and multiphase flows: a review

David Le Touzé<fnref n="1" />,*, ID and Andrea Colagrossi<fnref n="2" /> ID

E-mail: david.letouze@ec-nantes.fr and andrea.colagrossi@cnr.it

Received 6 May 2024, revised 13 November 2024
Accepted for publication 9 January 2025
Published 11 February 2025

Corresponding editor: Dr Masud Mansuripur

<figure><img src="image_3.png" /></figure>

CrossMark

## Abstract

The smoothed particle hydrodynamics (SPH) method is expanding and is being applied to more and more fields, particularly in engineering. The majority of current SPH developments deal with free-surface and multiphase flows, especially for situations where geometrically complex interface configurations are involved. The present review article covers the last 25 years of development of the method to simulate such flows, discussing the related specific features of the method. A path is drawn to link the milestone articles on the topic, and the main related theoretical and numerical issues are investigated. In particular, several SPH schemes have been derived over the years, based on different assumptions. The main ones are presented and discussed in this review underlining the different contexts and the ways in which they were derived, resulting in similarities and differences. In addition, a summary is provided of the recent corrections proposed to increase the accuracy, stability and robustness of SPH schemes in the context of free-surface and multiphase flows. Future perspectives of development are identified, placing the method within the panorama of Computational Fluid Dynamics.

Keywords: smoothed particle hydrodynamics (SPH), free-surface flow, multiphase flow, particle method, meshfree method

## Contents

1. Introduction

2. Governing equations for free-surface and multiphase flows

2

4

2.1. Boundary conditions

2.2. Equations of state (EOS)

2.3. NSEs in arbitrary-Lagrangian-Eulerian (ALE) framework

3. SPH interpolation

3.1. Mollification

3.1.1. Mollification in presence of interface

3.2. Smoothed spatial differential operators

3.2.1. Smoothed spatial differential operators in presence of interface

3.3. Quadrature

4. Solid boundary modeling

5. Particle approximation

6

6

7

8

8

8

9

10

10

12

13

<figure><img src="image_4.png" /></figure>

<footnote n="1">¹ Nantes Université, École Centrale Nantes, CNRS, LHEEA, UMR 6598, F-44000 Nantes, France</footnote>

<footnote n="2">² CNR-INM, INstitute of Marine engineering, Rome, Italy</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
<table>
  <tbody>
    <tr>
      <td>5.1. Density field approximation</td>
      <td>13</td>
    </tr>
    <tr>
      <td>5.2. Density field approximation in presence of interface</td>
      <td>14</td>
    </tr>
    <tr>
      <td>5.3. Density field approximation in presence of a free surface</td>
      <td>14</td>
    </tr>
    <tr>
      <td>6. Derivation of SPH equations: from a system of interacting macroscopic particles</td>
      <td>16</td>
    </tr>
    <tr>
      <td>6.1. Lagrangian thermo-mechanics</td>
      <td>16</td>
    </tr>
    <tr>
      <td>6.2. SPH pressure forces modeling</td>
      <td>17</td>
    </tr>
    <tr>
      <td>6.2.1. SPH pressure forces modeling in presence of a free surface</td>
      <td>17</td>
    </tr>
    <tr>
      <td>6.3. Viscous forces modeling</td>
      <td>18</td>
    </tr>
    <tr>
      <td>6.3.1. Heat conduction and viscous forces modeling in presence of an interface</td>
      <td>19</td>
    </tr>
    <tr>
      <td>6.4. SPH multiphase formulations</td>
      <td>20</td>
    </tr>
    <tr>
      <td>6.5. Conservation properties of the SPH scheme</td>
      <td>20</td>
    </tr>
    <tr>
      <td>6.6. Surface tension modeling from a microscopical description</td>
      <td>21</td>
    </tr>
    <tr>
      <td>6.7. Link between the SPH equations obtained and the NSEs</td>
      <td>21</td>
    </tr>
    <tr>
      <td>7. Derivation of SPH equations: from local NSEs</td>
      <td>22</td>
    </tr>
    <tr>
      <td>7.1. Consistency of the continuity equation in presence of free surface</td>
      <td>22</td>
    </tr>
    <tr>
      <td>7.2. Energy and momenta conservation</td>
      <td>22</td>
    </tr>
    <tr>
      <td>7.3. Viscous force modeling</td>
      <td>23</td>
    </tr>
    <tr>
      <td>7.4. Enforcement of the dynamic free-surface boundary condition</td>
      <td>24</td>
    </tr>
    <tr>
      <td>7.5. Heat diffusion modeling</td>
      <td>24</td>
    </tr>
    <tr>
      <td>7.6. Strong formulation of the SPH equations</td>
      <td>24</td>
    </tr>
    <tr>
      <td>7.7. SPH multiphase formulation</td>
      <td>24</td>
    </tr>
    <tr>
      <td>7.8. Surface tension modeling from a continuous description</td>
      <td>25</td>
    </tr>
    <tr>
      <td>7.8.1. CSF method</td>
      <td>25</td>
    </tr>
    <tr>
      <td>7.8.2. Reconstruction of a continuous interface from the particle arrangement</td>
      <td>25</td>
    </tr>
    <tr>
      <td>7.8.3. Surface tension through jump condition at the interfaces</td>
      <td>26</td>
    </tr>
    <tr>
      <td>7.8.4. Surface tension forces through minimization of the interface energy</td>
      <td>26</td>
    </tr>
    <tr>
      <td>8. Derivation of SPH equations: from continuous conservation laws</td>
      <td>26</td>
    </tr>
    <tr>
      <td>8.1. SPH multiphase formulation</td>
      <td>28</td>
    </tr>
    <tr>
      <td>9. Differences between the derived SPH models in sections 6–8</td>
      <td>28</td>
    </tr>
    <tr>
      <td>10. Particle arrangement within their Lagrangian evolution</td>
      <td>30</td>
    </tr>
    <tr>
      <td>10.1. Lagrangian trajectories, particle arrangement and interpolation accuracy</td>
      <td>30</td>
    </tr>
    <tr>
      <td>10.2. Remeshing</td>
      <td>30</td>
    </tr>
    <tr>
      <td>10.3. Particle Shifting Technique (PST)</td>
      <td>30</td>
    </tr>
    <tr>
      <td>10.4. Consistent particle shifting in the ALE framework</td>
      <td>31</td>
    </tr>
    <tr>
      <td>10.5. Volume conservation issue</td>
      <td>32</td>
    </tr>
    <tr>
      <td>10.6. Detection of free-surface particles</td>
      <td>32</td>
    </tr>
    <tr>
      <td>11. Stability and convergence of the SPH model</td>
      <td>33</td>
    </tr>
    <tr>
      <td>11.1. Time integration of the SPH schemes</td>
      <td>33</td>
    </tr>
  </tbody>
</table>

<table>
  <tbody>
    <tr>
      <td>11.2. Stability of the scheme adding numerical dissipation</td>
      <td>33</td>
    </tr>
    <tr>
      <td>11.3. Tensile instability</td>
      <td>35</td>
    </tr>
    <tr>
      <td>11.4. Pairing instability</td>
      <td>36</td>
    </tr>
    <tr>
      <td>11.5. Convergence of the SPH schemes</td>
      <td>36</td>
    </tr>
    <tr>
      <td>12. Modeling liquid phase: weakly-compressible vs. incompressible approach</td>
      <td>37</td>
    </tr>
    <tr>
      <td>12.1. Explicit vs. semi-implicit solution</td>
      <td>37</td>
    </tr>
    <tr>
      <td>12.2. Differences in the physical modeling</td>
      <td>37</td>
    </tr>
    <tr>
      <td>12.3. Liquid impacts</td>
      <td>38</td>
    </tr>
    <tr>
      <td>13. SPH simulation of liquid-gas flows</td>
      <td>38</td>
    </tr>
    <tr>
      <td>13.1. The choice of the speeds of sound</td>
      <td>39</td>
    </tr>
    <tr>
      <td>13.2. Simulation of highly-compressible gas phases</td>
      <td>39</td>
    </tr>
    <tr>
      <td>14. Examples of free-surface and multiphase flows solved by SPH models</td>
      <td>39</td>
    </tr>
    <tr>
      <td>15. Conclusion and perspectives</td>
      <td>43</td>
    </tr>
    <tr>
      <td>15.1. Conclusion</td>
      <td>43</td>
    </tr>
    <tr>
      <td>15.2. Perspectives</td>
      <td>44</td>
    </tr>
    <tr>
      <td>Data availability statement</td>
      <td>45</td>
    </tr>
    <tr>
      <td>Acknowledgments</td>
      <td>45</td>
    </tr>
    <tr>
      <td>References</td>
      <td>45</td>
    </tr>
  </tbody>
</table>

# 1. Introduction

In various engineering fields, the fluid dynamics of free-surface and multiphase flows plays an important role. Multiphase flows are also often referred to as interfacial flows, with a very similar meaning.

In marine engineering, examples cover the dynamics of coastal ocean waves, the dynamics of waves generated by ships and high-speed boats, or the interaction of large sea waves with marine vehicles, off-shore or coastal structures. In all of these cases, the air-water interface can be characterized by large deformation including breaking, fragmentation and reconnection phenomena, generally with a fast dynamics, i.e. short characteristic times.

Besides ocean waves, these kinds of geometrically complex free-surface fast-dynamics flows occur in other hydrodynamic impact situations: e.g. (i) in the aerospace and aeronautical field with the ditching of space aircraft, airplanes or helicopters, (ii) in the hydraulic field with dam breaks and in Pelton turbines, (iii) in industrial processes with waterjet cutting, etc.

Another example of such flows concerns the various situations involving dynamic motions of free-surface liquids in tanks. These situations can be found typically in the transportation of fuel, by trucks, inside ships, or in airplanes where the fuel is stored in the wings. In all these cases the induced sloshing can lead to severe loads and destabilizing motions because of external conditions to which the vehicle is subjected. This phenomenon is also a critical issue in reservoirs in the event of an earthquake, especially in the nuclear and oil-and-gas industries.

More generally, complex interface or free-surface flows with fast dynamics are found in many engineering fields such as automotive engineering (hydroplaning, gear box flow, wading), environmental engineering (landslides, oil spills, lava
flows), coastal (flooding) and hydraulic engineering (dam breaking, debris flows, Pelton turbines, fish passes), or industrial processes (casting, mixing, printing, additive manufacturing), among many other examples.

Physically speaking, when a liquid-gas interface breaks, the gas entrapment, the intrinsic three-dimensional nature of the flow and the many scales involved continue to make it a very complex phenomenon to investigate, both experimentally and numerically.

The first attempts to numerically tackle these types of flows date back to the 80s-90s when the volume-of-fluid (VoF) and then the level-set (LS) algorithms were developed to track liquid-gas interfaces within mesh-based solvers such as finite volume, finite difference and finite element methods (FEMs), see, e.g. [116, 225]. These VoF and LS methods have several accuracy and conservation issues, which have been improved in recent decades, e.g. through conservation corrections or through adaptive mesh refinement (AMR) techniques which can automatically increase the spatial resolution in the proximity of the interface, see, e.g. [42, 80, 241, 281]. Conversely, in Lagrangian meshless methods, often referred to as particle methods, such as the smoothed particle hydrodynamics (SPH) method, free surfaces and interfaces: (i) are intrinsically followed in their motion whatever their geometrical complexity, and (ii) remain perfectly sharp since interfaces exist between each pair of adjacent discrete elements of different species (e.g. between a gas element and a liquid one).

Besides the necessity of correctly solving the deformation of the fluid interfaces, another issue arises in mesh-based methods when solid surfaces have structural deformations or when (possibly many) solid bodies have relative motions with contacts, as in, e.g. a hydroplaning situation, the lubrication of a gear-box, debris flows, etc. In such conditions mesh-based methods require complex adaptations (e.g. mesh creation/destruction, penalty methods, augmented Lagrangian, immersed boundaries...) to handle mesh interpenetration and contact, and fail to accurately address the most complex situations. Conversely, these problems can be tackled straightforwardly with particle methods thanks to their meshless nature, where contact or solid body deformation induce no additional complexity. As a consequence, in a fluid-structure problem it is also easier to couple the particle method used to solve the fluid with, e.g. the FEM for structural deformations or the Discrete Element Method for multi-body interactions including contact.

It thus clearly appears that particle methods in general, and the SPH method in particular, present interesting assets to simulate geometrically-complex problems when this complexity originates from the fluid interfaces configuration and/or from the solid bodies in interaction with the fluids. In this context, the present review article is dedicated to how the SPH method was developed and used for simulating free-surface and multiphase flows. In particular, the milestone works that have allowed the development of the SPH method are identified and commented on, and the present status of the method regarding the simulation of this kind of flows is described.

It is worth noting that the study and application of the SPH method have received continuously growing interest in the past

decades. Twenty years ago, the total number of publications on 'Smoothed Particle Hydrodynamics' was less than 1000; by 2024, this number had grown to about 10 000. In the same space of time the number of citations grew by a factor of 10, as shown in figure 1 (data from Web of Science), reaching 30 000 citations per year. Importantly, the majority of the published papers on the SPH method today deal with free-surface and multiphase flows.

The diffusion of the SPH method in the last decades is also due to the development of various open-source and commercial software tools. Among the former a historically important one is GADGET [232] which deals with cosmological N-body/SPH simulations on massively parallel computers with distributed memory. Note that its updated version GADGET-2 was used for running the famous 'Millennium Simulation' in which more than 10 billion particles were used to trace the evolution of a portion of the Universe [230]. In the same field of astrophysics and cosmology another open-source solver is SWIFT [223] which is known for a famous simulation of a plausible scenario of formation of the moon as the result of an impact during the early days of our Solar System [128]. For free-surface and multiphase flows the most popular open-source solvers are DualSPHysics [60] and GPUSPH [26], which are able to run SPH simulations on Graphics Processing Units (GPUs) with a high number of particles (in the order of $10^6-10^9$).

Other open-source SPH software include: AQUAgpusph [39] based on the open cross-platform parallel programming OpenCL™, SPHinXsys [277] which adopts a multi-physics and multi-resolution library, the research code SPHERA [7] and PySPH [217] which uses a Python-based framework. Recently developed commercial SPH software are PreonLab and Simcenter-SPH-Flow [268].

In the literature there already exist other extended reviews dedicated to the SPH method, e.g. [25, 183, 231] in the astrophysical context where the method was initially developed. Monaghan [184] extended the method for the first time to free-surface flows, and also published two more reviews [187, 188] dedicated to the theory and application of the SPH method in the numerous areas where it was successfully applied.

In the review by Liu et al [153] the theory of the so-called 'particle approximation' was stressed, as well as solid boundary treatment techniques which have become a key issue in the last 20 years in enabling the extension of the SPH method to wider engineering fields.

In the last decade review articles on the SPH method focused mainly on free-surface and multiphase flows, and industrial applications related to these kinds of flows, see, e.g. [100, 104, 159, 226, 264]. Still in the last decade some others specifically addressed multiphase flows in the SPH context, e.g. [95, 247, 266], where the different solutions proposed in the literature for handling large density ratios and modeling surface tension were reviewed. Books on the SPH method are also present in the literature, see [152, 261]. In particular, in the latter the modeling of turbulent flows within the SPH method is discussed.

The most recent review articles are dedicated to the advanced versions of the SPH method in which the particle
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Number of publications and related citations from 1971 to 2022 on the topic ‘Smoothed Particle Hydrodynamics’ (data from Web of Science).</figcaption></figure>

shifting technique (PST) is embedded in the model to reduce
the disorder of the displacement of the particles in the fluid
domain, see, e.g. [271, 150]. In [150] the coupling of SPH with
other methods and the relevance of massive parallelization of
SPH schemes are also addressed. Other recent review articles
are dedicated to modeling and application to fluid-structure
interaction problems, either fully through the SPH method for
the fluid and the structure, or through coupling with the FEM
for the structure, see, e.g. [154, 276, 278].

When considering the specific context of free-surface and
multiphase flows, special attention needs to be paid to how
their different particularities are modeled within the SPH
schemes:

(1) presence of a contact discontinuity,

(2) use or not of a single-phase approximation,

(3) representation of liquid impacts at the interface within
a weakly-compressible or an incompressible approxima-
tion,

(4) accurate representation of the Lagrangian kinematics of
the interface,

(5) consistent enforcement of interface conditions including
energy exchanges, etc.

These aspects are not covered in detail in the review articles
mentioned above, whereas they are addressed in several recent
articles. In the present review these key articles are considered
and discussed to highlight the relevant steps forward which led
to the success and diffusion of the SPH method in the specific
free-surface hydrodynamics context, where gas-liquid inter-
faces play a relevant role in many real-life applications.

In the majority of SPH articles dedicated to free-surface or
multiphase flows, numerous theoretical and numerical confu-
sions are found, concerning both the above mentioned physical
and numerical particularities of the studied flows, and the core
of the SPH method itself. The combination of the Lagrangian

and meshless nature of the SPH method, indeed, makes its
numerical development and analysis very difficult.

Nowadays, several complementary points of view and
avenues of development exist:

(i) the SPH equations can be derived as a Hamiltonian system, or through the discretization of classical partial differential equations (PDEs) where the flux exchanges can be solved using Riemann solvers,

(ii) the liquid phase can be treated with a weakly-
compressible assumption or with an incompressible one,

(iii) the particles can be advected using a pure Lagrangian framework or an arbitrary Lagrangian Eulerian (ALE) one.

The present review article reports the main milestones
related to the different complementary points of view lis-
ted above, which we consider relevant to the current under-
standing of the method. Significant literature from the last
25 years is coherently linked providing a clearer understand-
ing of the SPH method when applied to free-surface and
multiphase flows, potentially aiding more meaningful future
developments in this specific framework.

Note that the overview of the SPH method given in the
present review mainly relies on the historic development of
compressible and weakly-compressible SPH models. More
recently, incompressible SPH models have been more and
more developed as well, which is reported in section 12.
However, the present review cannot be considered as compre-
hensive regarding this specific topic of incompressible SPH
models which may deserve a dedicated review.

## 2. Governing equations for free-surface and multiphase flows

A fluid domain $\Omega$ is here considered. Generally, it can be com-
posed of fluids of different natures and in different states:
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Sketch of the problem: the domain $\Omega$ is confined by an outer solid surface, different fluids a, b, c are present with an empty region e, confined by free surfaces $\partial\Omega_F$. Solid bodies can also be present inside the domain and their surfaces are indicated with $\partial\Omega_B$. All the interfaces between the fluids are oriented with normal vectors.</figcaption></figure>

liquids and gases. The fluid domain is confined by a boundary $\partial\Omega$ which can be composed of solid surfaces $\partial\Omega_B$ and free surfaces $\partial\Omega_F$. Hereinafter the term free surface will refer to a liquid-gas interface where the gas phase is not considered, i.e. an interface between liquid and void. The term interface refers to any liquid-liquid or liquid-gas interface where both fluids are considered. Figure 2 shows a sketch of the fluid

domain. The whole domain $\Omega$ is here confined by a solid surface; however, for some applications it can also be composed with inflow-outflow boundaries. Phase changes and miscibility phenomena are not considered.

The flow evolution is governed by the compressible Navier-Stokes equations (NSEs):

$$ 
\left\{ 
\begin{aligned}
\frac{D\rho}{Dt} &= -\rho \nabla \cdot \mathbf{u}, \\
\rho \frac{De}{Dt} &= -p \nabla \cdot \mathbf{u} + \mathbf{V} : \mathbf{D} - \nabla \cdot \mathbf{q},
\end{aligned}
\right. 
\qquad
\frac{D\mathbf{u}}{Dt} = -\frac{\nabla p}{\rho} + \frac{\nabla \cdot \mathbf{V}}{\rho} + \mathbf{f}, \\
\frac{Dr}{Dt} = \mathbf{u}, \quad p = p(\rho, e), 
\qquad (1)
$$

where $\mathbf{u}$ is the flow velocity, $\mathbf{r}$ the position of the material points, $p$ the fluid pressure, $\rho$ the fluid density, $e$ the specific internal energy, $\mathbf{V}$ the viscous stress tensor, $\mathbf{D}$ the rate-of-strain tensor, $\mathbf{f}$ a generic specific body force and finally $\mathbf{q}$ is the heat flux. Volumetric heat sources are not considered in the present review. Note that bold fonts are used for vectors and bar-bold fonts are used for higher-order tensors.

The derivative $D/Dt$ in equations (1) is the Lagrangian derivative of a generic field, namely:

$$ \frac{D(\bullet)}{Dt} := \frac{\partial(\bullet)}{\partial t} + \nabla(\bullet) \cdot \mathbf{u}. \qquad (2) $$

The fluids in $\Omega$ are considered to be Newtonian, so that the stress tensor $\mathbf{T}$ can be written as:

$$ \mathbf{T} = -p\mathbf{I} + \mathbf{V}, \quad \mathbf{V} = \lambda \text{tr}\mathbf{D}\mathbf{I} + 2\mu\mathbf{D}, \qquad (3) $$

where $\mu$ and $\lambda$ are the dynamic viscosity coefficients of the specific fluid, while $\mathbf{I}$ is the unit tensor. An additional hypothesis is that the viscous coefficients $\mu$ and $\lambda$ are assumed to

be constant in space and time (i.e. they are not linked with the temperature field).

The ratio between the infinitesimal mass, $\delta m$, and volume, $\delta V$, of the material point at position $\mathbf{r}$ is finite and equal, by definition, to the density field $\rho = \delta m/\delta V$. In particular, the mass of the material point is conserved along the Lagrangian trajectory, i.e. $D(\delta m)/Dt = 0$. Combining the latter relation with the continuity equation in equations (1), we obtain:

$$ \frac{D(\delta m)}{Dt} = \frac{D(\rho \delta V)}{Dt} = 0 \quad \Rightarrow \quad \frac{1}{\delta V} \frac{D \delta V}{Dt} = \nabla \cdot \mathbf{u}, \qquad (4) $$

sometimes referred to as the *volumetric strain-rate equation*, see, e.g. [139].

The Lagrangian framework is very useful when considering particle methods since in these models the fluid domain $\Omega$ is subdivided into small elements of fluid, referred to as *particles*, whose trajectories follow the velocity field, see section 5.
## 2.1. Boundary conditions

On solid surfaces $\partial\Omega_B$ the kinematic condition reads:

$$u = u_B, \quad \forall r \in \partial\Omega_B, \qquad (5)$$

due to the adhesion of the fluid on these boundaries.

When the Reynolds number is high, as is often the case in SPH applications, a perfect fluid assumption can be adopted to avoid the direct simulation of thin boundary-layer regions. The boundary condition is then a free-slip condition:

$$u \cdot n_B = u_B \cdot n_B, \quad \forall r \in \partial\Omega_B, \qquad (6)$$

where $n_B$ is the outward normal to the solid boundary. This perfect fluid assumption at the boundary is only valid for problems where the viscous and turbulent effects at the boundary have a negligible influence on the flow evolution.

Considering two generic fluids, $\mathcal{X}$ and $\mathcal{Y}$, a kinematic and a dynamic boundary condition must be satisfied on the interface between them $\partial\Omega_{\mathcal{X}\mathcal{Y}}$. Due to the hypotheses made above, of immiscible fluids and no phase changes, $\partial\Omega_{\mathcal{X}\mathcal{Y}}$ has zero thickness, thus containing zero mass, momentum and energy, see, e.g. [208]. It follows that the kinematic boundary condition requires the adhesion between material points belonging to the interface, i.e.:

$$u_{\mathcal{X}} = u_{\mathcal{Y}}, \quad \forall r \in \partial\Omega_{\mathcal{X}\mathcal{Y}}. \qquad (7)$$

In the Lagrangian framework the continuity of the normal velocity at the interface is intrinsic as well as its displacements. The tangential velocity is then assumed to be continuous at the interface. This no-slip condition assumption is not an obvious fact; indeed, it was once the subject of a long debate [208].

When $\mathcal{Y}$ is a void region, the interface $\partial\Omega_{\mathcal{X}\mathcal{Y}}$ is a free surface, and the kinematic condition reduces to only the continuity of the normal velocity.

Concerning the dynamic condition it expresses the stress balance at the interface between $\mathcal{X}$ and $\mathcal{Y}$:

$$(T_{\mathcal{X}} - T_{\mathcal{Y}}) \cdot n_{\mathcal{X}\mathcal{Y}}(\mathbf{r}) = \sigma_{\mathcal{X}\mathcal{Y}}(\mathbf{r}) \kappa_{\mathcal{X}\mathcal{Y}}(\mathbf{r}) n_{\mathcal{X}\mathcal{Y}}(\mathbf{r}) \\ + \nabla_{\mathcal{X}\mathcal{Y}}\sigma_{\mathcal{X}\mathcal{Y}}(\mathbf{r}), \quad \forall \mathbf{r} \in \partial\Omega_{\mathcal{X}\mathcal{Y}}, \qquad (8)$$

$\sigma_{\mathcal{X}\mathcal{Y}}(\mathbf{r})$ is the local surface tension coefficient between the generic fluids $\mathcal{X}$ and $\mathcal{Y}$ at the interface point $\mathbf{r}$, where $\mathcal{Y}$ can be the void when the fluid $\mathcal{X}$ is bounded by $\partial\Omega_F$. $\kappa_{\mathcal{X}\mathcal{Y}}(\mathbf{r})$ is the curvature at $\mathbf{r}$, $n_{\mathcal{X}\mathcal{Y}}(\mathbf{r})$ the corresponding unit normal vector pointing from $\mathcal{X}$ to $\mathcal{Y}$, and $\nabla_{\mathcal{X}\mathcal{Y}}$ the surface gradient on the interface.

Alternatively, the boundary condition (8) can be reshaped to see the surface tension as a concentrated volumic force $\rho f^\sigma$ acting on interfaces of $\Omega$, including free surfaces $\Omega_F$:

$$\rho f^\sigma(\mathbf{r}) = [\sigma_{\mathcal{X}\mathcal{Y}}(\mathbf{r}) \kappa_{\mathcal{X}\mathcal{Y}}(\mathbf{r}) n_{\mathcal{X}\mathcal{Y}}(\mathbf{r}) + \nabla_{\mathcal{X}\mathcal{Y}}\sigma_{\mathcal{X}\mathcal{Y}}(\mathbf{r})] \delta_{\mathcal{X}\mathcal{Y}}, \\ \kappa_{\mathcal{X}\mathcal{Y}}(\mathbf{r}) = -\nabla \cdot n_{\mathcal{X}\mathcal{Y}}(\mathbf{r}), \qquad (9)$$

where $\delta_{\mathcal{X}\mathcal{Y}}$ is the surface delta-function. The volumic force field $\rho f^\sigma(\mathbf{r})$ is then added to the right-hand side of the

momentum equation. This is useful from a numerical point of view for smearing the interface into a region of finite thickness, see sections 6.6 and 7.8. The last term in equations (8) and (9) linked to the local gradients in surface tension is connected to the Marangoni effect (see, e.g. [208, 255]) which is not considered in the present review.

As stated above, when $\mathcal{Y}$ is a void region the above relations are valid for free surfaces $\partial\Omega_F$. In the case of free surface and when $\sigma_{\mathcal{X}\mathcal{Y}} = 0$, it will be shown in section 7 how the SPH models can implicitly satisfy this condition in an integral sense, so that an explicit detection of the particles belonging to $\partial\Omega_F$ is not needed, neither is the enforcement of equation (8). This is a peculiar characteristic of this method concerning other numerical methods which implies an inherent ability to solve complex free-surface problems.

## 2.2. Equations of state (EOS)

To close the system (1) a link between pressure, density and internal energy is needed.

Differently from the gaseous phase, for the liquid phase in Standard Temperature and Pressure conditions the effect of the entropy on the pressure field can be ignored, therefore liquids can be seen as barotropic, i.e. $p = p(\rho)$ (which does not imply that the flow is isothermal).

For the gaseous phase, the typical EOS used are the stiffened and ideal gas ones, especially for flows with strong compressibility effects, see, e.g. [214, 237]. When limited compressibility occurs, both liquid and gaseous phases can be considered as barotropic. In particular the following polytropic EOS is used:

$$p = P_{0\mathcal{Y}} \left[ \left( \frac{\rho}{\rho_{0\mathcal{Y}}} \right)^{\gamma\mathcal{Y}} - 1 \right] + p_0, \quad P_{0\mathcal{Y}} = \frac{c_{0\mathcal{Y}}^2 \rho_{0\mathcal{Y}}}{\gamma\mathcal{Y}}, \qquad (10)$$

where $\rho_{0\mathcal{Y}}$ is the density at rest of the generic gas $\mathcal{Y}$, and $\gamma\mathcal{Y} > 0$ is its adiabatic index. Finally, $P_{0\mathcal{Y}}$ is the gas ambient pressure, while the constant $p_0$ is just used to re-scale the pressure field; indeed this can be useful for numerical reasons, see section 12. It is worth noting that the speed of sound $c_{0\mathcal{Y}}$ of the gas, and consequently its compressibility, is related to the ambient pressure $P_{0\mathcal{Y}}$. Besides, because of the hypothesis made on the nature of the fluids, in equations (1) the energy equation is decoupled from the momentum equation.

Following the idea in [183], equation (10) can also be used for the generic liquid $\mathcal{X}$:

$$p = \frac{c_{0\mathcal{X}}^2 \rho_{0\mathcal{X}}}{\gamma\mathcal{X}} \left[ \left( \frac{\rho}{\rho_{0\mathcal{X}}} \right)^{\gamma\mathcal{X}} - 1 \right] + p_0, \qquad (11)$$

where the $\gamma_{\mathcal{X}}$ parameter is quite high compared to that of the gaseous phases; for example for water it ranges between 6 and 7. Again $c_{0\mathcal{X}}$ is the speed of sound and $\rho_{0\mathcal{X}}$ is a reference density. The latter can be the density on the free surface where the pressure is constant and equal to $p_0$. Differently to a gas, the sound speed $c_{0\mathcal{X}}$ is not linked to the ambient pressure. Furthermore, except in rare situations, a liquid is always in a weakly-compressible regime, i.e. density variations are always smaller than 1% with respect to the reference density
$\rho_{0x}$. Within this regime a simple linearization of equation (11) can be used, see, e.g. [12]:

$$p = c_{0x}^2 (\rho - \rho_{0x}) + p_0. \quad (12)$$

In some conditions, the speeds of sound $c_0$ of the liquids and gases can be modified for computational convenience, see section 13.1. However, in conditions where, e.g. the gas is entrapped by the liquid during an impact event, the real speed of sound of the gas $c_{0y}$ has to be used since in that case, the compressibility effects on the entrapped gas pocket can be relevant, and these effects are mainly linked to the $c_{0y}$ value, which needs therefore to be evaluated correctly. Conversely, for the liquid phases, $c_{0x}$ can be reduced up to the limit of remaining in the weakly-compressible regime, a condition ensured by the following constraint:

$$c_{0x} \ge 10 U_{\text{ref}}, \\ U_{\text{ref}} := \max \left( \Delta U_{\text{max}}, \sqrt{\frac{\Delta P_{\text{max}}}{\rho_{0x}}}, \sqrt{\frac{\sigma_{XY}}{\rho} |\kappa_{\text{max}}|} \right), \quad (13)$$

where $|\kappa_{\text{max}}|$, $\Delta U_{\text{max}}$, $\Delta P_{\text{max}}$ are the maximum expected interface curvature (in absolute value), velocity and pressure variations during the time evolution. In other words, the Mach number defined as $Ma = U_{\text{ref}}/c_{0x}$ has to stay below 0.1 since the density variation is a function of $Ma^2$, see, e.g. [139, 164].

Another possible strategy when the compressibility of the fluid can be neglected is to consider the flow as incompressible. In that case, the velocity field is divergence free and the NSEs (1) can then be recombined to get a pressure Poisson equation (PPE), see, e.g. [139] (see also section 12).

## 2.3. NSEs in arbitrary-Lagrangian-Eulerian (ALE) framework

Conversely to the Lagrangian framework discussed above, in the most recent SPH models the governing equations (1) are rewritten in an ALE framework, introduced in [115]. In this framework any fluid element is advected with a modified velocity $(\mathbf{u} + \hat{\mathbf{u}})$, $\hat{\mathbf{u}}$ being an arbitrary velocity field. Fluid elements are thus no longer material elements, but elementary volumes $\delta V$ with mass fluxes through their boundaries $\partial(\delta V)$. Applying the Leibniz-Reynolds transport theorem to $\delta V$ reads:

$$\frac{d}{dt} \int_{\delta V} f \rho \, d\mathbf{r} = \int_{\delta V} \frac{\partial (\rho f)}{\partial t} \, d\mathbf{r} + \int_{\partial(\delta V)} (\rho f) (\mathbf{u} + \hat{\mathbf{u}}) \cdot \mathbf{n} \, d\mathbf{r}. \quad (14)$$

Using the continuity equation and the divergence theorem it can be re-shaped as:

$$\frac{d}{dt} \int_{\delta V} f \rho \, d\mathbf{r} = \int_{\delta V} \frac{Df}{Dt} \rho \, d\mathbf{r} + \int_{\partial(\delta V)} (\rho f) \hat{\mathbf{u}} \cdot \mathbf{n} \, d\mathbf{r}. \quad (15)$$

The last term in the right-hand side of equation (15) is the flux across the surface $\partial(\delta V)$ of the elementary volume due to the non-Lagrangian velocity $\hat{\mathbf{u}}$. After using the divergence theorem on this surface term, the integral equation (15) can be

collapsed assuming that $\delta V$ is small enough for the fields $\rho, f, \mathbf{r}, \mathbf{u}, \hat{\mathbf{u}}$ to be considered constant in it. Under such a hypothesis equation (15) becomes:

$$\left\{ \begin{aligned} \frac{\hat{D}(\rho f \delta V)}{Dt} &= \frac{Df}{Dt} \rho \delta V + \nabla \cdot (\rho f \hat{\mathbf{u}}) \delta V, \\ \frac{\hat{D} \mathbf{r}}{Dt} &= \mathbf{u} + \hat{\mathbf{u}}. \end{aligned} \right. \quad (16)$$

It is worth noting that in this equation the advective derivative $\hat{D}(\bullet)/Dt$ is with respect to the modified advection velocity $(\mathbf{u} + \hat{\mathbf{u}})$, i.e.:

$$\frac{\hat{D}f}{Dt} := \frac{\partial f}{\partial t} + \nabla f \cdot (\mathbf{u} + \hat{\mathbf{u}}) = \frac{Df}{Dt} + \nabla f \cdot \hat{\mathbf{u}}. \quad (17)$$

Substituting $f = 1/\rho$, $f = 1$, $f = \mathbf{u}$ and $f = e$ in the first equation of equations (16), the time variations of, respectively, the element volume $\delta V$, mass $\delta m$, momentum $(\delta m \mathbf{u})$ and internal energy $(\delta m e)$ are obtained. Finally, combining the latter with equations (1), it gives:

$$\left\{ \begin{aligned} \frac{\hat{D}(\delta V)}{Dt} &= \nabla \cdot \mathbf{u} \delta V + \nabla \cdot \hat{\mathbf{u}} \delta V, \\ \frac{\hat{D}(\delta m)}{Dt} &= \nabla \cdot (\rho \hat{\mathbf{u}}) \delta V, \\ \frac{\hat{D}(\delta m \mathbf{u})}{Dt} &= -\nabla p \delta V + (\nabla \cdot \mathbf{V}) \delta V + f \delta m \\ &\quad + f \sigma \delta m + \nabla \cdot (\rho \mathbf{u} \otimes \hat{\mathbf{u}}) \delta V, \\ \frac{\hat{D}(\delta m e)}{Dt} &= -p (\nabla \cdot \mathbf{u}) \delta V \\ &\quad + (\mathbf{V} : \mathcal{D}) \delta V - (\nabla \cdot \mathbf{q}) \delta V + \nabla \cdot (\rho e \hat{\mathbf{u}}) \delta V, \\ \frac{\hat{D} \mathbf{r}}{Dt} &= \mathbf{u} + \hat{\mathbf{u}}, \quad \rho = \frac{\delta m}{\delta V}, \quad p = p(\rho, e). \end{aligned} \right. \quad (18)$$

Equations (18) become purely Lagrangian if $\hat{\mathbf{u}}$ is null (i.e. $\hat{D}f/Dt = Df/Dt$ and equations (1) are found again) or purely Eulerian if $\hat{\mathbf{u}} = -\mathbf{u}$ (i.e. $\hat{D}f/Dt = \partial f/\partial t$).

Within the ALE framework $\hat{\mathbf{u}}$ is arbitrary and this can induce large variations in elementary masses and volumes. Consequently, the numerical schemes derived by equations (18) need to be carefully handled to avoid instabilities and inaccuracies. In the recently enhanced SPH models, the velocity $\hat{\mathbf{u}}$ is generally small in comparison to $\mathbf{u}$ and is just used to regularize the particle distributions as discussed in section 10, i.e. a quasi-Lagrangian formalism.

As a final comment, since in the ALE framework the advocation of elementary volumes is made with velocity $\mathbf{u} + \hat{\mathbf{u}}$, to verify the kinematic boundary condition on the fluid interfaces $\partial\Omega_{XY}$ (see section 2.1), the $\hat{\mathbf{u}}$ component needs also to fulfil $\hat{\mathbf{u}} \cdot \mathbf{n} = 0$ on $\partial\Omega_{XY}$ (see also section 10.3).
### 3. SPH interpolation

#### 3.1. Mollification

A generic field $f$ is considered with the equality:

$$f(\mathbf{r}) = \int_{\Omega} f(\mathbf{r}') \delta(\mathbf{r} - \mathbf{r}') \, d\mathbf{r}', \quad (19)$$

where $\delta(\mathbf{r} - \mathbf{r}')$ is the delta Dirac function and the convolution integral is extended over the fluid domain $\Omega$. In the SPH method $f$ is approximated through:

$$\langle f \rangle (\mathbf{r}) = \int_{\Omega} f(\mathbf{r}') W(\mathbf{r} - \mathbf{r}'; h) d\mathbf{r}'. \quad (20)$$

In the latter the $\delta$ distribution of equation (19) is approximated with a smoothing function, $W(\mathbf{r} - \mathbf{r}'; h)$. The latter is also called kernel function and different choices are possible for it. For the sake of simplicity in the present work the function $W$ has the following characteristics:

(1) It is a $C^k$ function, with $k \ge 2$, defined on a spherical compact support $\Omega_W(\mathbf{r})$ centered at the position $\mathbf{r}$. The characteristic length of the support is $h$ which is called smoothing length which is here assumed to be constant in space and time. The radius of $\Omega_W(\mathbf{r})$, $R_W$, is equal to $h$ multiplied by a number typically equal to 2 or 3. The smoothing length can be regarded as a measure of the domain of influence of the elementary volume in $\mathbf{r}$.

(2) $W$ is positive and decreases monotonously with $\tilde{r} := ||\mathbf{r} - \mathbf{r}'||$ to zero on the border $\partial\Omega_W(\mathbf{r})$ of its support.

(3) It converges in the weak sense to the Dirac measure when its characteristic length $h$ tends to zero.

(4) Its integral on the domain $\Omega$ is unitary, i.e.:

$$\Gamma(\mathbf{r}) := \int_{\Omega} W(\mathbf{r} - \mathbf{r}'; h) \, d\mathbf{r}' = 1. \quad (21)$$

(5) The kernel function is radially centered, so that:

$$\nabla W(\mathbf{r} - \mathbf{r}'; h) = -(\mathbf{r} - \mathbf{r}') F(\tilde{r}; h), \\ F(\tilde{r}; h) := \frac{\nabla W \cdot (\mathbf{r}' - \mathbf{r})}{\tilde{r}^2}, \quad (22)$$

where $\nabla$ is the spatial gradient with respect to the variable $\mathbf{r}$. It follows that $W$ is symmetric with respect to the variables $\mathbf{r}, \mathbf{r}'$, so that:

$$W(\mathbf{r} - \mathbf{r}'; h) = W(\mathbf{r}' - \mathbf{r}; h) \\ \text{and} \quad \nabla W(\mathbf{r} - \mathbf{r}'; h) = - \nabla' W(\mathbf{r} - \mathbf{r}'; h), \quad (23)$$

where $\nabla'$ is the spatial gradient with respect to the variable $\mathbf{r}'$. The function $F(\tilde{r}; h)$ is, like $W$, a positive function.

Different choices for the kernel function can be made. Recently, the Wendland function $\psi_{31}$, see [267], has become popular in the SPH community because of its benefits from the numerical point of view [67] (see also section 11.4). The $\psi_{31}$ function is also named the C2 Wendland kernel for its

continuous second derivatives. Following, e.g. [262], the C2 Wendland kernel reads as:

$$W(q) = \beta \begin{cases} (2-q)^4 (1+2q) & \text{for } 0 \le q \le 2 \\ 0 & \text{for } q > 2 \end{cases}, \\ F(q) = \frac{10\beta}{h^2} \begin{cases} (2-q)^3 & \text{for } 0 \le q \le 2 \\ 0 & \text{for } q > 2 \end{cases}, \quad (24)$$

where $q = \tilde{r}/h$ while $\beta$ is a normalization constant equal to $7/(64\pi h^2)$ in a 2D and $21/(256\pi h^3)$ in a 3D framework. $R_W = 2h$ for this specific kernel function.

In a more general context the smoothing length can be a function of time and space as reported for example in [30, 112, 187]. Indeed, the possibility to use variable $h$ in time and space is useful when considering a highly-compressible flow, for example in the astrophysical context, or for integrating shallow water equations. The main subsequent difficulties are linked to the restoring of the symmetry properties listed above since the radius of the kernel support changes point by point.

For the purpose of the present review these kinds of models are not treated here. However, having a variable $h$ in space can also be useful to have different spatial resolutions in $\Omega$. From the recent literature it was shown that it is efficient and accurate to subdivide $\Omega$ into domains each with constant $h$ which exchange information through overlapping regions, see section 3.3. This approach, proposed in [45] and further developed in [221], avoids the convolution between points with different $h$ which implies a more complex formulation when equation (20) is derived in space, involving space derivatives of $h$. Hereinafter, in the kernel function notation we implicitly assume its dependency on $(\mathbf{r} - \mathbf{r}')$ and on $h$, and also assume that $h$ is constant over the considered sub-domain.

Considering the finite value of $h$ and $\Omega_W(\mathbf{r}) \subset \Omega$ it also follows from the $W$ properties that:

$$\mathcal{I}_1 := 1 - \int_{\Omega} W d\mathbf{r}' = 0, \quad \mathcal{I}_2 := \mathbf{r} - \int_{\Omega} \mathbf{r}' W d\mathbf{r}' = 0, \\ \mathcal{I}_3 := \int_{\Omega} \nabla W d\mathbf{r}' = 0, \quad \mathcal{I}_4 := \mathcal{I} - \int_{\Omega} \mathbf{r}' \otimes \nabla W d\mathbf{r}' = 0. \quad (25)$$

For a regular enough function $f$, using the properties $\mathcal{I}_1$ and $\mathcal{I}_2$ and expanding the smoothed representation of $\langle f \rangle$ with a Taylor series, it is possible to show that, see, e.g. [173]:

$$\langle f \rangle (\mathbf{r}) = f(\mathbf{r}) + \mathcal{O}(h^2), \quad \forall \mathbf{r} \in (\Omega \setminus \partial\Omega). \quad (26)$$

It follows that if $f$ is a linear function $\langle f \rangle \equiv f$.

#### 3.1.1. Mollification in presence of interface.
It must be noted that the relation (26) is not true if $\mathbf{r} \in \partial\Omega$. In this case it becomes:

$$\langle f \rangle (\mathbf{r}) = f(\mathbf{r}) \bar{\Gamma} + \mathcal{O}(h), \quad \forall \mathbf{r} \in \partial\Omega, \quad \bar{\Gamma} := \lim_{h \to 0} \Gamma, \quad (27)$$

with $\bar{\Gamma}$ that is always less than one, as sketched in figure 3. If $\partial\Omega$ is a regular surface, $\bar{\Gamma} = 1/2$.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Sketch of the intersection between the kernel support with the fluid domain $\Omega$. To exemplify the convergence of this intersection two kernel support radii are drawn: $2h_1$ and $2h_2$, with $2h_1 > 2h_2$. Left: point close to the frontier $\partial\Omega$. Right: point belonging to the frontier $\partial\Omega$.</figcaption></figure>

If the field $f$ is null along the boundary the convergence is recovered. More generally, the non-convergence of equation (27) can be removed by normalizing the approximation $\langle f \rangle$ with the $\Gamma$ value:

$$ \langle f \rangle^S(\mathbf{r}) = \int_{\Omega} f(\mathbf{r}') W^S d\mathbf{r}', \quad W^S(\mathbf{r}-\mathbf{r}') := \frac{W(\mathbf{r}-\mathbf{r}')}{\Gamma(\mathbf{r})}. \qquad (28) $$

The above correction is known in the literature as the Shepard correction see, e.g. [93], and it allows for recovery of the first-order convergence of the mollification procedure since the new kernel $W^S$ satisfies the $\mathcal{I}_1$ property of equations (25). The second-order convergence can be recovered using the moving least square (MLS) correction described, e.g. in [93], through the definition of a new kernel $W^{\text{MLS}}$ which is built to satisfy all the properties of equations (25):

$$ \left\{ \begin{array}{ll} \langle f \rangle^{\text{MLS}}(\mathbf{r}) = \int_{\Omega} f(\mathbf{r}') W^{\text{MLS}} d\mathbf{r}', & W^{\text{MLS}} := \mathbb{A}^{-1}(\mathbf{r}) e_1 \cdot a W, \quad e_1^T := [1, 0, 0, 0], \\ a^T := [1, (x-x'), (y-y'), (z-z')], & \mathbb{A}(\mathbf{r}) := \int_{\Omega} (a \otimes a) W d\mathbf{r}', \end{array} \right. \qquad (29) $$

where $x, y, z$ and $x', y', z'$ are the coordinates of the points $\mathbf{r}$ and $\mathbf{r}'$, respectively. Since the corrected kernels $W^S$ and $W^{\text{MLS}}$ are based on procedures that depend on the specific point $\mathbf{r}$ treated, they do not satisfy the symmetry properties (23) and cannot therefore be used straightforwardly in the SPH schemes discussed in the following sections.

It is worth noting that the lack of consistency of $\langle f \rangle$ due to the invalidity of properties of equations (25) can also take place inside $\Omega$ if the function $f$ presents discontinuities inside the fluid domain. This can be the case of the density field when different phases or fluids are present in $\Omega$, see section 5.2.

## 3.2. Smoothed spatial differential operators

In order to approximate the differential operators present in the governing equations the mollification (20) is applied to the spatial derivative of the function $f$:

$$ \langle \nabla f \rangle (\mathbf{r}) = \int_{\Omega} \nabla' f(\mathbf{r}') W d\mathbf{r}', \qquad (30) $$

using the identities:

$$ \nabla f W = \nabla(f W) - f \nabla W, \quad \nabla(f W) = \nabla \cdot (f W \mathbb{I}), \qquad (31) $$

and applying the divergence theorem, equation (30) can be reshaped as:

$$ \begin{aligned} \langle \nabla f \rangle (\mathbf{r}) &= \int_{\Omega} \nabla' f(\mathbf{r}') W d\mathbf{r}' \\ &= \int_{\Omega} f(\mathbf{r}') \nabla W d\mathbf{r}' + \int_{\partial\Omega} f(\mathbf{r}') W n' d\mathbf{r}', \end{aligned} \qquad (32) $$

where $n'$ is a normal vector to $\partial\Omega$ pointing outside $\Omega$, while $\nabla$ indicates the spatial derivatives with respect to the variable $\mathbf{r}$ (the symmetric kernel property (23) has been used to change the spatial derivative of the kernel function $W$ from the variable $\mathbf{r}'$ to the variable $\mathbf{r}$).

The right-hand side of equation (32) allows the evaluation of the gradient of $f$ from the knowledge of the function itself, which is the key point of the SPH method as well as of all meshless methods. The surface integral in equation (32) is zero for all points $\mathbf{r}$ inside $\Omega$ which have their kernel supports which do not intersect the fluid domain boundary ($W$ is zero on $\partial\Omega_W$ as stated before).

Following [193] a mollification for the Laplacian operator can be derived as:

$$ \langle \nabla^2 f \rangle (\mathbf{r}) = 2 \int_{\Omega} [f(\mathbf{r}') - f(\mathbf{r})] F d\mathbf{r}', \qquad (33) $$
where the kernel derivative function F has been defined in equation (22).

Again if f is regular enough, using the properties (25) it is possible to show that:

$$ \begin{aligned} \langle \nabla f \rangle(\mathbf{r}) &= \nabla f(\mathbf{r}) + \mathcal{O}(h^2), \\ \langle \nabla^2 f \rangle(\mathbf{r}) &= \nabla^2 f(\mathbf{r}) + \mathcal{O}(h^2), \quad \forall \mathbf{r} \in (\Omega \setminus \partial\Omega). \end{aligned} \quad (34) $$

### 3.2.1. Smoothed spatial differential operators in presence of interface.

Considering equation (27) also for the mollified gradient we get an inconsistency for the points belonging to $\partial\Omega$, i.e.:

$$ \langle \nabla f \rangle(\mathbf{r}) = \nabla f(\mathbf{r}) \bar{\Gamma} + \mathcal{O}(h), \quad \forall \mathbf{r} \in \partial\Omega. \quad (35) $$

Furthermore, as shown in [49] for all the points belonging to the boundary $\partial\Omega$ the surface integral of equation (32) is no longer zero and diverges as $\mathcal{O}(1/h)$.

In the case where $\nabla f$ is zero on $\partial\Omega$, equation (35) converges as $\mathcal{O}(h)$. Further, rewriting equation (32) as:

$$ \langle \nabla f \rangle(\mathbf{r}) = \int_{\Omega} [f(\mathbf{r}') - f(\mathbf{r})] \nabla W d\mathbf{r}' + \int_{\partial\Omega} [f(\mathbf{r}') - f(\mathbf{r})] W n' d\mathbf{r}', \quad (36) $$

it is shown in [49] that the second term is $\mathcal{O}(h)$ so that, if omitted, it does not compromise the convergence of the first term. Summarizing, where $\nabla f$ is zero on $\partial\Omega$:

$$ \langle \nabla f \rangle(\mathbf{r}) = \int_{\Omega} [f(\mathbf{r}') - f(\mathbf{r})] \nabla W d\mathbf{r}', \quad (37) $$

is $\mathcal{O}(h^2)$ when $\mathbf{r} \in (\Omega \setminus \partial\Omega)$ and $\mathcal{O}(h)$ when $\mathbf{r} \in \partial\Omega$.

Regarding the Laplacian operator, similarly we can conclude that, on the boundary $\partial\Omega$, equation (33) diverges as $\mathcal{O}(1/h)$ due to the absence of surface terms. If $\nabla f$ is zero on $\partial\Omega$, $\langle \nabla^2 f \rangle$ converges to the limit $\bar{\Gamma}$. If in addition the Laplacian of $f$ is zero on the domain boundary, $\langle \nabla^2 f \rangle$ converges as $\mathcal{O}(h)$.

To overcome such inconsistencies a corrected kernel gradient $\nabla W^{\text{MLS}}$ can be used. Indeed, equation (29) can be extended to evaluate the corrected gradient, see [71]. An equivalent correction is also derived in [218] combining equation (37) with a correction of the kernel, which reads:

$$ \left\{ \begin{aligned} \langle \nabla f \rangle^L(\mathbf{r}) &:= \int_{\Omega} [f(\mathbf{r}') - f(\mathbf{r})] \nabla^L W d\mathbf{r}', & \nabla^L W &:= \mathbb{L}(\mathbf{r}) \cdot \nabla W, \\ \mathbb{L}(\mathbf{r}) &:= \left[ \int_{\Omega} (\mathbf{r}' - \mathbf{r}) \otimes \nabla W d\mathbf{r}' \right]^{-1} = \left[ \int_{\Omega} (\mathbf{r}' - \mathbf{r}) \otimes (\mathbf{r}' - \mathbf{r}) F(||\mathbf{r}' - \mathbf{r}||) d\mathbf{r}' \right]^{-1}, \end{aligned} \right. \quad (38) $$

where the kernel property (22) is used in the last term of the second equation which shows that the tensor $\mathbb{L}$ is symmetric (being $F(||\mathbf{r}' - \mathbf{r}||) = F(||\mathbf{r} - \mathbf{r}'||)$).

### 3.3. Quadrature

In this section the convolution integrals introduced in the previous sections are discretized, with subsequent discretization errors that combine with the previous ones related with the mollification approximation, i.e. the ones depending on $h$. The combination of these two sources of errors is briefly discussed here.

The domain $\Omega$ is discretized into elementary volumes, the generic $i$th of them is centered at position $\mathbf{r}_i$ and has a volume $V_i$. The simplest discretization of equations (20), (37), (33) is then:

$$ \begin{aligned} \langle f \rangle_i &= \sum_j f_j W_{ij} V_j, & \langle \nabla f \rangle_i &= \sum_j (f_j - f_i) \nabla W_{ij} V_j, \\ \langle \nabla^2 f \rangle_i &= 2 \sum_j (f_j - f_i) F_{ij} V_j, \end{aligned} \quad (39) $$

where a compact notation is adopted hereinafter, indicating the dependency on two points $\mathbf{r}_i$ and $\mathbf{r}_j$ just with an index notation. Note that in the present article indices always refer to

positions and/or particle labels and thus never to vector or matrix components.

In equation (39) the summations apply to all the $j$ points inside the support of the kernel centered in $\mathbf{r}_i$. It is important to stress that in this kind of approximation the calculation point positions and the volumes are not directly linked, see section 5. From the elementary volume $V_i$ the average distance:

$$ \Delta r_i := (V_i)^{1/n}, \quad (40) $$

is extracted to associate a size to the element centered in $\mathbf{r}_i$, where $n$ is the dimension of the considered physical space.

In the context of interest of this work (limited compressibility of the flow) the points are distributed in configurations where the ratio $\Delta r/h$ is subjected to low variations within the subdomains where $h$ is constant. In other words, in those subdomains, all the points have almost the same $\Delta r$. For a highly-compressible flow, this condition is no longer true so that $h$ should vary if one maintains the ratio $\Delta r/h$ close to a constant, or elementary volumes need to be created or deleted in the subdomains during their expansion or contraction motions, see section 13.2.

The discrete interpolation errors in equation (39) are strongly linked to the spatial organization of the interpolation
points, as stressed in, e.g. [48, 215, 10]. These errors can be measured at a given point i by calculating how the four sum-mations:

$$I_1 := 1 - \sum_j W_{ij} V_j, \qquad I_2 := \mathbf{r}_i - \sum_j \mathbf{r}_j W_{ij} V_j, \qquad (41)$$

$$I_3 := \sum_j \nabla W_{ij} V_j, \qquad I_4 := \mathbb{I} - \sum_j \mathbf{r}_j \otimes \nabla W_{ij} V_j,$$

are different from zero. Equations (41) are the discrete version of equations (25). The level of those errors is dependent on two parameters: (i) the number of points used in the quadrature, i.e. the number of elements in the sums, directly linked to the $\Delta r/h$ ratio, and (ii) the regularity of the spacing of these quadrature points to properly reproduce the integrated function variations, i.e. the spatial organization of the interpolation points. Of course, the more varying the function, the greater the error, so that the $I_4$ error will a priori be higher than, in order, $I_3$ to $I_1$ for a given distribution of the points in the sums. If all the four I-terms are exactly zero the discrete interpolations are exact for linear function (first order completeness, see [22]).

The global error of SPH approximations in equation (39) thus depends on three parameters: $h$ for the mollification error in continuous space, and $\Delta r/h$ and $\xi$ for the discretization error, where $\xi$ is a dimensionless measure of the non-uniformity of the calculation points distribution. Since the main space operators in the Euler equations are a gradient and a divergence, studying this error is especially important for the first-order derivative of functions. This was studied by [215] who proved that, in a one-dimensional problem:

$$\sum_j f_j \nabla W_{ij} V_j - f_i' = \underbrace{\frac{f_i}{h} \left[ \xi \mathcal{O} \left( \left( \frac{\Delta r}{h} \right)^3 \right) + \frac{1}{2} \left( \xi^2 + \frac{1}{12} \right) \mathcal{O} \left( \left( \frac{\Delta r}{h} \right)^4 \right) \right]}_{\textcircled{1}} + f_i' \underbrace{\left[ \xi \mathcal{O} \left( \left( \frac{\Delta r}{h} \right)^3 \right) + \mathcal{O} \left( \left( \frac{\Delta r}{h} \right)^4 \right) \right]}_{\textcircled{2}} + f_i'' h \underbrace{\left[ \xi \mathcal{O} \left( \frac{\Delta r}{h} \right) + \mathcal{O} \left( \left( \frac{\Delta r}{h} \right)^4 \right) \right]}_{\textcircled{3}} + \underbrace{\frac{h^2}{4} f_i'''' + \dots}_{\textcircled{4}} \qquad (42)$$

Extension to two or three dimensions yields similar results. If $\xi$ and $\Delta r/h$ are sufficiently small the I-term errors are small enough for the convergence with $h$ to be dominated by the last term ④ in equation (42). After which, for smaller values of $h$, the first (non convergent or divergent) term ① becomes dominant. Based on this observation SPH practitioners generally fix the number of neighbors (i.e. the ratio $h/\Delta r$) high enough and seek a distribution of calculation points that is as regular as possible (low $\xi$, see section 10), in order to remain convergent with $h$ for the range of space resolutions they use.

Figure 4 shows graphs from [215] of the interpolation errors using a sinusoidal function as a test function. On the left is plotted the interpolation error on a Cartesian lattice, i.e. $\xi = 0$ in equation (42). On the right plot the error is shown for different values of the particle spacing perturbation $\sigma$ (an empirical measure of the parameter $\xi$) and fixing the ratio $h/\Delta r$.

Note that this condition worsens with the order $m$ of derivation of the function since in equivalent developments as of equation (42) the first term ① would go with $h^{-m}$.

Note also that, if we use the gradient approximation $\langle \nabla f \rangle_i$ of equations (39), i.e. with $(f_j - f_i)$ instead of just $f_j$ in equation (42), one can prove that one order of convergence can be gained. In fact, considering the function $f=1$ in equation (42) it follows that:

$$\sum_j \nabla W_{ij} V_j = \frac{1}{h} \left[ \xi \mathcal{O} \left( \left( \frac{\Delta r}{h} \right)^3 \right) + \frac{1}{2} \left( \xi^2 + \frac{1}{12} \right) \mathcal{O} \left( \left( \frac{\Delta r}{h} \right)^4 \right) \right], \qquad (43)$$

which, multiplied by $f_i$ and subtracted from equation (42) yields:

$$\sum_j (f_j - f_i) \nabla W_{ij} V_j - f_i' = \textcircled{2} + \textcircled{3} + \textcircled{4} + \dots, \qquad (44)$$

where the divergent term ① disappears.

If the errors linked to the I-terms are too high, or if the support of the kernel intersects the domain boundary $\partial\Omega$, the consistency can be recovered through the corrected kernels $W^S$ or $W^{MLS}$ introduced in the previous section 3.1.1, which read in discrete form as:

$$W_{ij}^{S} := W_{ij} / \Gamma_{i}, \qquad \Gamma_{i} = \sum_{j} W_{ij} V_{j}; \qquad (45)$$

$$\left\{ \begin{array}{l} W_{ij}^{\text{MLS}} := \mathbb{A}_i^{-1} \mathbf{e}_1 \cdot \mathbf{a}_{ij} W_{ij}, \qquad \mathbf{e}_1^\text{T} := [1, 0, 0, 0], \\ \mathbf{a}_{ij}^\text{T} := [1, (x_i - x_j), (y_i - y_j), (z_i - z_j)], \qquad \mathbb{A}_i := \sum_j (\mathbf{a}_{ij} \otimes \mathbf{a}_{ij}) W_{ij} V_j. \end{array} \right. \qquad (46)$$
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Analytically and empirically calculated L2-norm error of SPH estimates of the first derivative of A(x) = A₀ sin(2πx/λ). Left: as a function of smoothing length, for various values of particle spacing. Right: as a function of smoothing length, for various particle spacing perturbations σ, computed with Δx/h = 0.364. [215] John Wiley & Sons. Copyright © 2005 John Wiley & Sons, Ltd.</figcaption></figure>

Equation (45) makes it possible to recover the property $I_1 = 0$ while with equation (46) all the I-terms are null by construction. Analogously for the gradient function the corrected kernel (38) can be rewritten in discrete form as:

$$ 
\left\{ 
\begin{array}{l} 
\langle \nabla f \rangle_i^L := \mathbb{L}_i \left[ \sum_j f_{ji} \nabla W_{ij} V_j \right], \quad \mathbb{L}_i := \left[ \sum_j (\mathbf{r}_{ij} \otimes \mathbf{r}_{ij}) F_{ij} V_j \right]^{-1}, \\ 
f_{ji} := (f_j - f_i), \quad \mathbf{r}_{ij} := (\mathbf{r}_i - \mathbf{r}_j), \quad F_{ij} := F(||\mathbf{r}_{ij}||). 
\end{array} 
\right. 
\quad (47) 
$$

The above renormalized gradient is exact when f is a linear function. Note that since $F_{ij} = F_{ji}$ it follows that $\mathbb{L}_i$ is symmetric.

A similar correction can be derived for the Laplacian approximation, see [87]. More sophisticated possibilities inspired by mesh-based techniques can also be used to increase the order of the interpolation, with a different perspective. In fact, in the context of the solution of compressible flows with the use of Riemann solvers (see section 11.2), increasing the method order can be obtained by increasing the accuracy of the reconstruction of the variables entering the Riemann problem at the interface between two calculation points (seen as a one dimensional problem). This principle is used in, e.g. [18, 257].

## 4. Solid boundary modeling

The difficulty linked to boundary modeling with the SPH method is related to two main elements. First, in a meshless method there is no surface discretization of the domain boundaries associated with the volume discretization of the interior of the domain, contrary to mesh-based methods where both go together. Second, in the presence of a boundary the kernel support is truncated when evaluated at or in the vicinity of the boundary. This kernel support truncation has two consequences on the SPH interpolation in the vicinity of the boundary (see section 3), which are already present in continuous space, i.e. after mollification but before discretization: (i)

the $I_1$ condition in equations (25) is no longer respected, with an error tending to 1/2 when approaching the boundary and, (ii) the boundary integral of equation (32) is non-zero since the truncated support is no longer compact, and must thus be accounted for. To address this difficulty three main ways arose in the SPH literature, historically:

(i) The first way is to simply introduce boundary forces through 'repulsive particles' placed on the boundary, thus avoiding having to account for the latter two sources of error. These forces partly counteract the error introduced by the kernel support truncation, especially preventing the fluid particles from escaping the domain, see, e.g. [185, 190]. The advantage of this technique is its ease of implementation for complex geometries but its major drawback is that it is not conservative and of poor accuracy.

(ii) Later, authors followed another way which is to fill the missing kernel support part with fictitious particles inside the boundary. These fictitious particles are named in different ways depending on the authors. This straightforwardly cancels the aforementioned sources of error but introduces the difficulty of placing and assigning physical properties to these fictitious particles. Different variants were proposed in which the fictitious particles are either moving or fixed inside the boundary and with different ways to extract fictitious particle properties from fluid particle information, see, e.g. [2, 54, 61, 79, 90, 148, 167, 272]. This technique has better conservation properties but has the drawback of requiring arbitrary algorithmic choices close to geometrical singularities. Imposing boundary conditions onto a volume rather than at the boundary is also challenging. The left and middle plots of figure 5 depict two sketches of the so-called Ghost-Particles approach taken from [168]. Fictitious particles (ghosts) are positioned within the solid volume and their pressures and velocities are given in such a way that the solid boundary conditions can be enforced on the fluid particles. An extension of the above approach is
<figure><img src="image_9.png" /><figcaption>**Figure 5.** SPH solid boundary modeling. Left and middle: ghost particles approach. Reprinted from [168], Copyright (2013), with permission from Elsevier. Right: boundary integral term approach. Reproduced from [138], with permission from Springer Nature.</figcaption></figure>

given in [15] through the Clone-Particles technique where the main idea relies on a local definition of the normal and tangent vectors to the solid body along with a local mirroring of the flow field for each fluid particle.

(iii) The most recent, advanced and consistent way is to discretize the boundary integral term of equation (32) and to compensate the missing kernel support part by a correction as in equation (45). However, this introduces the difficulty of evaluating a boundary integral in the context of a meshless method. Different authors have proposed variants to calculate the correction and the surface integral, analytical or discrete, with time integration or not, with panel discretization of the boundary or not, and with different ways to extrapolate pressure information at the boundary, see, e.g. [66, 70, 88, 138, 165]. This technique solves the drawbacks of the two previous ones but introduces an increased algorithmic complexity, leading to the need for efficient implementation for complex geometries [43]. The right plot of figure 5 shows a sketch taken from [138] for the evaluation of the boundary integral term of equation (32).

Note that some authors followed slightly different paths such as, e.g. introducing an immersed boundary technique [113] inspired from mesh-based methods. Note also that modeling boundary presence and conditions remains an open and active branch of research on the SPH method [251].

## 5. Particle approximation

In particle methods for continuous media, the main idea is to consider the medium $\Omega$ as a collection of small macroscopic pieces of material named particles, see [137]. The whole set of particles represents a partition of the material domain. The particles move according to the velocity field $\mathbf{u}$, for example the trajectory of the ith particle is integrated through $dr_i/dt = \mathbf{u}_i$, and each particle carries a quantity $\mathcal{F}$ along its trajectory.

For example, in two dimensional Particle Vortex Methods the continuous vorticity field $\omega(\mathbf{r}, t)$ is approximated using a set of discrete elementary vortices, each having a circulation

$\Gamma_i$ which is transported along its Lagrangian trajectory. The link between the vorticity and the circulations is then given through a convolution with a suitable kernel function $W$:

$$ \langle \omega \rangle (\mathbf{r}, t) = \sum_{j} \Gamma_{j} W(\mathbf{r} - \mathbf{r}_{j}), \quad \text{Particle Vortex Method.} \qquad (48) $$

Similarly, in the SPH method the continuous density field $\rho(\mathbf{r}, t)$ is approximated using a set of discrete elementary particles, each having a mass $m_i$ which is transported along its Lagrangian trajectory:

$$ \langle \rho \rangle (\mathbf{r}, t) = \sum_{j} m_{j} W(\mathbf{r} - \mathbf{r}_{j}), \quad \text{SPH Method.} \qquad (49) $$

The latter equation corresponds to the SPH interpolation (39) where $m_j = \rho_j V_j$.

Thanks to this approach the particle masses do not change in time in the SPH method and mass conservation is intrinsically guaranteed. This is one of the main advantages of this method with respect to other numerical methods, especially in the presence of interfaces in a multiphase context, since in that context mass conservation is guaranteed in each phase by the SPH method.

Equation (49) is the historical way of assessing the density in the SPH method in [98, 158]. The evaluation of the density is actually a central point for the method which determines then the particle evolution equations and their conservation properties, see section 6. From there, further developments and variants have then been proposed by different authors. The present section discusses this topic.

### 5.1. Density field approximation

When particle densities are obtained by equation (49) and their masses are known from the initial conditions, their volumes can be simply evaluated as:

$$ \rho_i(t) = \sum_j m_j W_{ij} \Rightarrow V_i(t) = \frac{m_i}{\rho_i(t)} = \frac{m_i}{\sum_j m_j W_{ij}}, \qquad (50) $$
where the time dependency is implicitly linked to the particle positions $r_i$ and $r_j$. In this relation the particle volume $V_i$ is not evaluated by a geometrical procedure but through an approximation weighting the particle positions and masses in the support $\Omega_W(r_i)$ of the kernel function $W$. Through this approximation it is thus not possible to take into account how the particle volumes are deformed by the flow and this lack can lead to inconsistencies between particle positions, masses and densities. This drawback is discussed more in detail in section 10. It is also worth noting that this representation is not a 'partition of unity', namely the sum of all $V_j$ in $\Omega$ is not equal to the volume of the fluid domain, which is further discussed in the next subsection.

The density field approximation in equation (50) is useful to solve astrophysical problems involving gas dynamics in highly-compressible regime, for which the SPH method was originally developed, see [98, 158]. In that context particle distribution largely varies in space, even at the scale of the kernel support. For this reason, as stressed for example in [112], the smoothing length $h$ cannot be constant.

Note that alternative definitions of the volume can be used instead of the one in equation (50), see [38], so as, e.g. to take into account non-uniform mass distributions, to suppress spurious surface tension at contact discontinuities, or to build a preconditioner to define consistent initial conditions.

## 5.2. Density field approximation in presence of interface

In the context of the present article the density approximation (50) has the drawback that it smooths the interfaces between fluids having different densities. In other words equation (50) is not able to maintain density discontinuities within the fluid domain $\Omega$ and thus, cannot be used in such a context as remarked in [54].

Beside astrophysics, the SPH method was also used for stochastic simulations for studying the dynamic and rheological properties of fluids at the mesoscopic scale. Within this field of application, where the compressibility effects are limited, a different formula was derived in [81] for linking the particle densities with their mass spatial distribution:

$$ \rho_i(t) = m_i \sum_j W_{ij} \Rightarrow V_i(t) = \frac{m_i}{\rho_i(t)} = \frac{1}{\sum_j W_{ij}}. \quad (51) $$

This second formula is based on the assumption that the particle masses do not vary much within the kernel support.

It is interesting to see that in equation (51) the particle volumes are given by kernel summation only and are therefore just a function of the particle positions. Again, it is worth noting that this representation is not a 'partition of unity', and that the error is linked to the ratio $\Delta r/h$. In other words the errors on the volumes are always linked to the quadrature error and not to the size of the kernel support.

For example, if an $L \times L$ two-dimensional square domain is considered and discretized by $N$ particles arranged in a

<table><caption>**Table 1.** Relative error $\epsilon_V$ on the volume of a square domain $L \times L$, varying the ratio $\Delta r/h$, i.e. the number of particles present in the kernel support $\Omega_W$.</caption><thead><tr><td>Δr/h</td><td>No. of neighbors (2D)</td><td>ε<sub>V</sub></td></tr></thead><tbody><tr><td>1</td><td>13</td><td>3.4%</td></tr><tr><td>1/2</td><td>50</td><td>0.12%</td></tr><tr><td>1/4</td><td>201</td><td>0.0039%</td></tr></tbody></table>

Cartesian lattice configuration, using equation (51) and definition (41) it follows that:

$$ \begin{aligned} V_i &= \frac{\Delta r^2}{\sum_j W_{ij} \Delta r^2} = \frac{L^2/N}{(1 - I_1)}, \\ \epsilon_V &= \frac{L^2 - \sum_i V_i}{L^2} = \frac{-I_1}{(1 - I_1)} \approx -I_1, \end{aligned} \quad (52) $$

where $\epsilon_V$ is the relative error on the evaluation of the whole domain volume, i.e. $L^2$, through particle volumes summation. Error $\epsilon_V$ is given by the quadrature error $I_1$ which is primarily linked to the ratio $\Delta r/h$; hence if the latter is constant the error $\epsilon_V$ remains constant even when increasing the number of particles $N$. Table 1 reports the error $\epsilon_V$ for three different $\Delta r/h$ ratios; the second column refers to the corresponding number of particles present in the kernel support $\Omega_W$ in a 2D framework.

For a more general case, where the domain shape is more complex and the particle distribution is not simply Cartesian, errors linked to the spatial resolution $L/\Delta r$, i.e. to the number of particles $N$ in the domain, are also present. However, in the limit of large enough $N$ the error will be governed by the quadrature error, i.e. by the ratio $\Delta r/h$, and this is a general rule for the SPH interpolation as already stressed in section 3.

Since the volumes $V_i$ are linked to the kernel and not to a geometrical tessellation of $\Omega$, in [81] these volumes are called 'thermodynamic volumes'; indeed, in section 6 it is shown that equation (51) determines the particle interaction forces.

The main advantage of the density approximation in equation (51) with respect to that in equation (50) is that it is able to preserve density discontinuities in $\Omega$. In fact, when the kernel support intersects an interface between fluids, equation (51) only considers the geometrical contribution of the particles located on the other side of the interface, not their mass. Thanks to this property the SPH model presented in [81] is used in [120, 121] to simulate multiphase flows with large density ratios.

## 5.3. Density field approximation in presence of a free surface

Contrary to astrophysical applications where the density decays to zero at the domain boundaries, the density approximations (50) and (51) cannot be used for free-surface flows since, when the support of the kernel intersects $\partial\Omega_F$, the kernel summation induces a non-physical decay of the particle densities and volumes.

To overcome this limit for simulating free-surface flows, in [183] a solution was given by deriving equation (50) in time.
From the ordinary time derivative of the kernel function, see [51]:

$$ \frac{dW_{ij}}{dt} = \mathbf{u}_{ij} \cdot \nabla W_{ij} = -\mathbf{u}_{ij} \cdot \mathbf{r}_{ij} F_{ij}, \quad \mathbf{u}_{ij} := (\mathbf{u}_i - \mathbf{u}_j), \quad (53) $$

it follows that the time derivative of equation (50) reads:

$$ \frac{d\rho_i}{dt} = -\sum_j m_j \mathbf{u}_{ij} \cdot \mathbf{r}_{ij} F_{ij}, \quad (54) $$

which can be seen as a discrete continuity equation, making it possible to evaluate the particle densities by integrating them in time. The advantage of this approach is that the right-hand side of equation (54) in the continuous limit $\Delta r/h \to 0$, i.e.:

$$ \lim_{\Delta r/h \to 0} \frac{d\rho_i}{dt} = - \int_{\Omega} \rho(\mathbf{r}') [\mathbf{u}(\mathbf{r}_i) - \mathbf{u}(\mathbf{r}')] \cdot (\mathbf{r}_i - \mathbf{r}') F \, d\mathbf{r}', \quad (55) $$

is an approximation of $(\rho \nabla \cdot \mathbf{u})$ consistent in the limit $h \to 0$ even on the free surface. Therefore the problem of the kernel truncation at the free surface can be overcome thanks to this property.

In order to demonstrate such a property it is enough to substitute the following approximations for the velocity and density:

$$ \begin{aligned} \mathbf{u}(\mathbf{r}') &= \mathbf{u}(\mathbf{r}_i) + \nabla \mathbf{u}(\mathbf{r}_i) \cdot [\mathbf{r}' - \mathbf{r}_i] + \mathcal{O}(\bar{r}^2), \\ \rho(\mathbf{r}') &= \rho(\mathbf{r}_i) + \nabla \rho(\mathbf{r}_i) \cdot [\mathbf{r}' - \mathbf{r}_i] + \mathcal{O}(\bar{r}^2), \end{aligned} \quad (56) $$

into equation (55). Using the kernel properties it follows that:

$$ \lim_{\Delta r/h \to 0} \frac{d\rho_i}{dt} = - (\rho \nabla \mathbf{u})_i : \int_{\Omega} (\mathbf{r}' - \mathbf{r}_i) \otimes (\mathbf{r}' - \mathbf{r}_i) F d\mathbf{r}' + \mathcal{O}(h). \quad (57) $$

Considering a point $\mathbf{r}_i$ on the free surface and the properties of the kernel function (see section 3), the above equation becomes:

$$ \lim_{\Delta r/h \to 0} \frac{d\rho_i}{dt} = -\bar{\Gamma} (\rho \nabla \cdot \mathbf{u})_i + \mathcal{O}(h), \quad \forall \mathbf{r}_i \in \partial \Omega_F. \quad (58) $$

Since the free surface is an isobaric surface and the fluids considered are barotropic with polytropic laws (see section 2), $\partial \Omega_F$ is also an iso-density surface. The velocity divergence is zero on it and it is not necessary to normalize equation (55) with the factor $\bar{\Gamma}$.

Note that equation (54) can be rewritten as:

$$ \frac{d\rho_i}{dt} = -\sum_j \rho_j \mathbf{u}_{ij} \cdot \mathbf{r}_{ij} F_{ij} V_j = -\langle \rho \nabla \cdot \mathbf{u} \rangle_i, \quad (59) $$

which is the SPH version of the continuity equation in the governing equations (1).

It is worth noting that doing this comparison between equation (59) and the continuity equation in equations (1) the ordinary time derivative $d/dt$ is associated to the Lagrangian time derivative $D/Dt$.

The relation (59) is a key equation in the SPH model, since:

* We showed that it is consistent with the governing equation in the limits $h \to 0$, $\Delta r/h \to 0$ even on the free surface boundary.

* It is Galilean invariant.

* It links the kernel function $F$ with the velocity gradients and this property will be used in section 6.3 for linking $F$ also to the viscous forces model.

Even if equation (59) can be used for simulating free-surface flows, it cannot be adopted in problems where fluids of different densities are present in $\Omega$. To overcome this problem in [102] the time derivative is applied to equation (51) defined by [81] as:

$$ \frac{dV_i}{dt} = V_i^2 \sum_j (\mathbf{u}_{ij} \cdot \mathbf{r}_{ij}) F_{ij}. \quad (60) $$

In such a case the discrete version of the volumetric strain-rate equation (4) is used instead of the continuity equation. It follows that a second smoothed differential operator for the velocity divergence can be written as:

$$ \frac{dV_i}{dt} = V_i \langle \nabla \cdot \mathbf{u} \rangle_i, \quad \langle \nabla \cdot \mathbf{u} \rangle_i = V_i \sum_j (\mathbf{u}_{ij} \cdot \mathbf{r}_{ij}) F_{ij}. \quad (61) $$

Equation (60) can be used for both free-surface and multiphase flows and shares the three properties of equation (59) listed above.

Once the particle volumes are known through equation (60) the densities can be calculated directly with the ratio $\rho_i(t) = m_i/V_i(t)$. However, in [102] the particle densities are evaluated using the Shepard kernel defined in section 3.3:

$$ \rho_i = \sum_j m_j W_{ij}^S, \quad W_{ij}^S = W_{ij} / \sum_{k \in P_i} W_{ik} V_k, \quad (62) $$

where the summation on the right equation is extended only to the particles belonging to the same fluid as that of the ith particle (i.e. the set $P_i$). Equation (62) can be seen as the generalization of equation (50) where, thanks to the Shepard kernel, the drawback of the kernel truncation is overcome. However, it is worth noting that equation (62), conversely to equation (50), needs a priori knowledge not only of the particle positions but also of the particle volumes.

Thanks to the combination of equations (60) and (62) the SPH model presented in [102] was successfully applied to multiphase-flow problems where a free surface is also present, see section 8.1.

Note that in the presence of solid boundaries $\partial \Omega_B$ the kernel support also intersects the solid boundary for particles nearby, thus needing specific treatments if one wants to use any of the density approximations mentioned in this section, see section 4.
# 6. Derivation of SPH equations: from a system of interacting macroscopic particles

## 6.1. Lagrangian thermo-mechanics

In the present section, the particle motion equations are derived using the particle approximation presented in section 5 combined with the Lagrangian of the particle system.

Let us consider a macroscopic system of N fluid particles (fluid elements) in a physical space of n dimensions. The generic ith particle has a position $\mathbf{r}_i$, velocity $\mathbf{u}_i$, internal energy per unit mass $e_i$, entropy per unit mass $s_i$ and mass $m_i$. The latter is assumed to be constant in time. The fluid particles are large enough for thermal agitation not to be observed, and ignored. This perspective permits us to consider each particle as a thermodynamic system in itself, and thus thermodynamic relations are valid at the particle level. It is further assumed that only central forces between pairs of particles, directed along the line connecting their centers of mass, act as interparticle forces. No coupling between the translational motion of the particles and the possible vibrational or rotational degrees of freedom is hence considered.

The particle system considered can be described in the framework of Lagrangian and Hamiltonian mechanics. The Lagrangian $\mathcal{L}$ is a state function of its N generalized coordinates $(\mathbf{r}_i, \mathbf{u}_i)$ defined, at a certain time $t$, as:

$$ \mathcal{L} = \sum_i \left[ m_i \frac{\mathbf{u}_i^2}{2} - m_i U(t, \mathbf{r}_i) - m_i e(\rho_i, s_i) \right], \quad (63) $$

where $U$ is a generic potential energy per unit of mass whose gradient gives the force field $\mathbf{f} = -\partial U/\partial \mathbf{r}$, such as, e.g. gravity. The last term in the right-hand side contains the specific internal energy $e$ related to the thermodynamic part of the particle-particle interaction. The specific internal energy is assumed to be a function of the density and entropy fields through a suitable EOS $e = e(\rho, s)$.

The dynamics of the particle system is given by the equation:

$$ \frac{d}{dt} \left( \frac{\partial \mathcal{L}}{\partial \mathbf{u}_i} \right) - \frac{\partial \mathcal{L}}{\partial \mathbf{r}_i} = Q_i^V + Q_i, \quad (64) $$

where $Q_i^V$ are the generalized dissipative forces acting between the particles and $Q_i$ are the non-conservative forces plus the non-dissipative forces due to the Galilean constraints. The latter are linked to the boundary conditions enforced on the surface of the solid boundaries $\partial \Omega_B$, see section 2.1. Different techniques have been developed to model such kinematic constraints for the SPH models as summarized in section 4.

In order to close the system of ODEs (64) we need to specify the thermodynamic behavior of the system through the specific internal energy of a particle $e_i$. In accordance with the first law of thermodynamics the pressure $p_i$ of a particle is obtained through the variations of the specific internal energy as a consequence of variations of the density field, while in accordance with the second law of thermodynamics the temperature $T_i$ of a particle is obtained as the variation of the

specific internal energy as a consequence of variations of the entropy field:

$$ p_i = \rho_i^2 \frac{\partial e_i}{\partial \rho_i} \Big|_s, \quad T_i = \frac{\partial e_i}{\partial s_i} \Big|_V. \quad (65) $$

As stressed in section 2 the fluids considered are all barotropic in this work and the momentum equation (64) is decoupled from the energy equation, i.e.:

* the temperature field is not linked with the pressure,

* the temperature field changes only under the action of the viscous forces,

* the viscous forces do not depend on the temperature field.

Under the above assumptions, substituting equation (63) into equation (64), and using the first equation of equations (65) we get:

$$ \begin{gathered} m_i \frac{d\mathbf{u}_i}{dt} - m_i \mathbf{f}_i - m_i \mathbf{f}_i^P = m_i \mathbf{f}_i^V + m_i \mathbf{f}_i^\sigma + m_i \mathbf{f}_i^S, \\ m_i \mathbf{f}_i^P := -\sum_j m_j \frac{p_j}{\rho_j^2} \frac{\partial \rho_j}{\partial \mathbf{r}_i}, \quad m_i \mathbf{f}_i^V := Q_i^V, \quad m_i \mathbf{f}_i^S := Q_i, \end{gathered} \quad (66) $$

where the terms $\mathbf{f}^P, \mathbf{f}^V$ and $\mathbf{f}^\sigma$ represent the specific pressure, viscous and surface tension force interactions, respectively, while $\mathbf{f}$ is the specific body force due to the action of the potential field $U$ and $\mathbf{f}^S$ are the specific forces related to the constraints on the solid boundaries. In order to make the $\mathbf{f}^P$ term explicit, we need to define how the density field is connected with the particle spatial distribution, while the terms $\mathbf{f}^V$ and $\mathbf{f}^\sigma$ will be treated respectively in sections 6.3 and 6.6. In section 5 we showed that different equations can be used for the density evaluation. This implies that different SPH equations can be derived from equations (66), as discussed in section 6.2.

In association with the Lagrangian (63) of the system, the energy of the particle system is defined as:

$$ \mathcal{E} = \sum_i \left[ m_i \frac{\mathbf{u}_i^2}{2} + m_i U(t, \mathbf{r}_i) + m_i e(\rho_i, s_i) \right]. \quad (67) $$

Following [53], differentiating equation (67) with respect to the time and considering that the time derivative of the specific internal energy is expressed as:

$$ \frac{\mathrm{d}e_i}{\mathrm{d}t} = \left( \frac{p_i}{\rho_i^2} \right) \frac{\mathrm{d}\rho_i}{\mathrm{d}t} + T_i \frac{\mathrm{d}s_i}{\mathrm{d}t}, \quad (68) $$

the equation governing the evolution of the entropy of each particle can be derived:

$$ m_i T_i \frac{ds_i}{dt} = -m_i \mathbf{u}_i \cdot \mathbf{f}_i^V + m_i \frac{dq_i}{dt}, \quad (69) $$

where $q_i$ is the heat per unit mass exchanged by particle $i$ with the surroundings. In the modeling of $q_i$, see section 6.3.1, considering an isolated system $dq_i/dt$ should satisfy the condition $\sum_i m_i dq_i/dt = 0$ in order to maintain the global energy conservation.
Note that, unlike in a dissipative particle dynamics (DPD) model, see, e.g. [82, 83], in the SPH context the fluid particles are large enough for thermal agitation not to be observed, as already stated. Hence, the case where the terms $Q_i$ and $Q_i^Y$ are set to zero effectively corresponds to a discrete model of an inviscid fluid.

Additional considerations on energy aspects of the Lagrangian framework in particle methods can be found in [27, 85, 261].

## 6.2. SPH pressure forces modeling

In this section the specific pressure force interaction term $f_p$ of equations (66) is discussed. As stated before, this term can be made explicit by using the density field approximations discussed in sections 5.1 and 5.2. From these two different SPH density approximations the following specific pressure force terms are derived:

$$ egin{cases} 
\begin{array}{l} 
\text{SPH schemes for flows without a free surface} 
\end{array} \\
\begin{array}{l} 
S_1: \rho_i = \sum_j m_j W_{ij} \\ 
S_2: \rho_i = \frac{m_i}{V_i}, \quad V_i = \frac{1}{\sum_j W_{ij}} 
\end{array} 
\Rightarrow 
\begin{array}{l} 
m_i f_i^p = -m_i \sum_j \left( \frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2} \right) \nabla W_{ij} m_j, \\ 
m_i f_i^p = -\sum_j (p_i V_i^2 + p_j V_j^2) \nabla W_{ij}. 
\end{array} 
\end{cases} 
\quad (70) $$

As stated in section 5: (i) the scheme $S_1$ can be used in the astrophysical context; it is suitable for compressible flow while it cannot be used in the context of multiphase flow; (ii) conversely, the scheme $S_2$ can be adopted for simulating multiphase flow but cannot be used in the context of compressible flow; (iii) neither of these two schemes can be used if a free surface is present at the boundary of the flow domain.

### 6.2.1. SPH pressure forces modeling in presence of a free surface.
To deal with free-surface flow the density approximations described in section 5.3, based on the time integration of the continuity equation, are needed. In order to use formulas (54) and (60) the total power of the pressure forces is used:

$$ \sum_i m_i f_i^P \cdot \mathbf{u}_i = - \sum_i \sum_j m_j \frac{p_j}{\rho_j^2} \frac{\partial \rho_j}{\partial \mathbf{r}_i} \cdot \mathbf{u}_i = \sum_i m_i \frac{p_i}{\rho_i^2} \left( \frac{d\rho_i}{dt} \right), \quad (71) $$

where in the last term the time derivative of the density can be recognized as:

$$ \frac{d\rho_i}{dt} = - \sum_j \frac{\partial \rho_i}{\partial \mathbf{r}_j} \cdot \mathbf{u}_j. \quad (72) $$

Substituting equations (54) and (60) into equation (71) two other SPH schemes can be derived:

$$ egin{cases} 
\begin{array}{l} 
\text{SPH schemes for flows with a free surface} 
\end{array} \\
\begin{array}{l} 
S_3: \frac{d\rho_i}{dt} = \sum_j m_j \mathbf{u}_{ij} \cdot \nabla W_{ij} \\ 
S_4: \frac{dV_i}{dt} = -V_i^2 \sum_j \mathbf{u}_{ij} \cdot \nabla W_{ij} 
\end{array} 
\Rightarrow 
\begin{array}{l} 
m_i f_i^p = -m_i \sum_j \left( \frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2} \right) \nabla W_{ij} m_j, \\ 
m_i f_i^p = -\sum_j (p_i V_i^2 + p_j V_j^2) \nabla W_{ij}. 
\end{array} 
\end{cases} 
\quad (73) $$
<table><caption>**Table 2.** The derived SPH schemes and the kind of flow they can solve.</caption><thead><tr><th>SPH scheme</th><th>Compressible flow</th><th>Multiphase flow</th><th>Free-surface flow</th></tr></thead><tbody><tr><td>S<sub>1</sub></td><td>Yes</td><td>No</td><td>No</td></tr><tr><td>S<sub>2</sub></td><td>No</td><td>Yes</td><td>No</td></tr><tr><td>S<sub>3</sub></td><td>Yes</td><td>No</td><td>Yes</td></tr><tr><td>S<sub>4</sub></td><td>No</td><td>Yes</td><td>Yes</td></tr></tbody></table>

As expected, the expressions of the pressure forces $f^p$ in the above equations (73) are the same as those of equations (70).

A similar derivation to the one presented in this section can be found in [29] where a variational procedure on the energy of the fluid particles was used.

Note that since the scheme $S_4$ is based on the particle volume distribution instead of the particle mass distribution, it can be adopted for simulating multiphase flow with large density ratio and where a free-surface boundary is also present, see [102]. It is important to stress that the summations in scheme $S_4$ comprise the particles belonging to different fluids, close to the interfaces. Thus, even though this scheme is able to preserve density discontinuity, the velocity divergence and the pressure gradient are smoothed across the interfaces. A different approach is used in [280] where those summations apply only to the particles belonging to the same fluid and to fictitious particles on the other side of the interface with extrapolated pressure and velocity. However, this second approach, even if more correct, can only be adopted in simple cases in practice, where the interface remains regular enough.

Summarizing, in table 2 the four derived SPH schemes are classified with respect to the flow condition they can apply to. For the topic of the present review the scheme $S_4$ is therefore the most suitable.

### 6.3. Viscous forces modeling

Following the work of Rayleigh, dissipative forces $f^N$ can be derived from a variational principle involving the so-called dissipation function $\Phi_D$. The dissipation function may be expressed as a quadratic form of the relative particle velocities, $\mathbf{u}_{ij} = \mathbf{u}_i - \mathbf{u}_j$:

$$ \Phi_D = \frac{1}{2} \sum_i \sum_j \mathbf{u}_{ij} \cdot \mathbf{A}_{ij}(r_{ij}) \cdot \mathbf{u}_{ij}, \quad (74) $$

where $\mathbf{A}_{ij}$ is a positive symmetric tensor. Moreover, the symmetry is extended to the particle indices, i.e. $\mathbf{A}_{ij} = \mathbf{A}_{ji}$, to ensure the action-reaction principle.

As shown at the end of this subsection $\Phi_D$ must be a positive function to comply with the requirements of the second principle of thermodynamics. According to this assumption, the specific viscous force $f^N$ takes the form:

$$ m_i f_i^N = - \frac{\partial \Phi_D}{\partial \mathbf{u}_i} = - \sum_j \mathbf{A}_{ij}(r_{ij}) \cdot \mathbf{u}_{ij}. \quad (75) $$

In the derivation described in [53, 260] it is shown that in order to guarantee:

* the preservation of the angular momentum, i.e. $\sum_i \mathbf{r}_i \times m_i \mathbf{f}_i^N = 0$,

* zero value of $f^N$ when the velocity field is a rigid body motion,

* finite value of $f^N$ when the inter-particle distance goes to zero,

* consistency with the dimension $n$ of the physical space,

the tensor $\mathbf{A}$ needs to take the form:

$$ \mathbf{A}_{ij} = g_{ij}(r_{ij}) \left( \frac{\mathbf{r}_{ij} \otimes \mathbf{r}_{ij}}{r_{ij}^2} \right), \quad g_{ij}(r_{ij}) = g_{ji}(r_{ij}) > 0, \quad (76) $$

where the condition of positiveness of $g_{ij}$ guarantees the positiveness of the tensor $\mathbf{A}$. Combining equations (75) and (76) one gets:

$$ m_i f_i^N = - \sum_j g_{ij}(r_{ij}) \pi_{ij} \mathbf{r}_{ij}, \quad \pi_{ij} := \left( \frac{\mathbf{u}_{ij} \cdot \mathbf{r}_{ij}}{r_{ij}^2} \right). \quad (77) $$

In order to derive the function $g_{ij}$ a limit process towards the governing equations (1) is used in [53], leading to:

$$ g_{ij}(r_{ij}) = K \mu_{ij} F(r_{ij}) V_i V_j, \quad K := 2(n+2), \quad (78) $$

where $\mu_{ij}$ is a dynamic viscosity coefficient. When particle $i$ and $j$ belong to the generic fluid $\mathcal{X}$, $\mu_{ij} = \mu_{\mathcal{X}}$. Conversely, the condition when $i$ and $j$ belong to different fluids will be treated in the next subsection 6.3.1.

Summarizing, for the SPH scheme viscous forces take this form:

$$ m_i f_i^N = -K \sum_j \mu_{ij} \pi_{ij} F_{ij} \mathbf{r}_{ij} V_i V_j. \quad (79) $$

Differently from the SPH pressure forces approximation presented above, which is linked to the density approximation adopted, the SPH viscous forces approximation derived here is independent from this choice and can be combined with all the four schemes discussed in the previous section 6.2.

Figure 6 shows a comparison, taken from [56], between an SPH model using the viscous formula (79) and a Finite Volume Method solution for the simulation of a viscous flow past a cylinder close to a free surface. The comparison is performed in terms of vorticity field, which highlights the complex wake interaction between the shear layers and the deformed free surface.

The power dissipated by the viscous forces (79) can be evaluated as:

$$ \mathcal{P}^N = - \sum_i m_i f_i^N \cdot \mathbf{u}_i. \quad (80) $$

Substituting equation (79) into (80) we can reshape it as:

$$ \mathcal{P}^N = \frac{K}{2} \sum_i \sum_j \mu_{ij} \pi_{ij} F_{ij} V_i V_j \left[ (\mathbf{u}_i + \mathbf{u}_j) \cdot \mathbf{r}_{ij} + (\mathbf{u}_i - \mathbf{u}_j) \cdot \mathbf{r}_{ij} \right]. \quad (81) $$
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Viscous flow past a cylinder close to a free surface. Vorticity field comparison between: SPH simulation (left) and Finite Volume Method coupled with a level-set function for tracking the free-surface evolution (right). Reprinted from [56], Copyright (2019), with permission from Elsevier.</figcaption></figure>

Considering that $\mu_{ij}$, $\pi_{ij}$, $F_{ij}$, $V_i V_j$ and $(\mathbf{u}_i + \mathbf{u}_j)$ are all symmetric terms with respect to swapping the indices $i$ and $j$ while $r_{ij}$ is anti-symmetric, it follows that:

$$ \sum_i \sum_j \mu_{ij} \pi_{ij} F_{ij} V_i V_j (\mathbf{u}_i + \mathbf{u}_j) \cdot \mathbf{r}_{ij} = 0. \quad (82) $$

Therefore, equation (81) can be simplified as:

$$ \mathcal{P}^V = \frac{K}{2} \sum_i \sum_j \mu_{ij} \pi_{ij} F_{ij} (\mathbf{u}_{ij} \cdot \mathbf{r}_{ij}) V_i V_j. \quad (83) $$

Considering the definition (77) of $\pi_{ij}$ and the definition (74) of the dissipation function it follows that:

$$ \mathcal{P}^V = \Phi_D = \frac{K}{2} \sum_i \sum_j \mu_{ij} (\pi_{ij} r_{ij})^2 F_{ij} V_i V_j > 0, \quad (84) $$

$F_{ij}$ being strictly positive $\Phi_D$ is too. Then using the entropy production equation (69) and summing all the particle contributions one finally gets:

$$ \sum_i m_i T_i \frac{ds_i}{dt} = - \sum_i m_i \mathbf{u}_i \cdot \mathbf{f}_i^V = \Phi_D > 0. \quad (85) $$

It follows that the dissipation function $\Phi_D$ defined by equation (84) satisfies the second law of thermodynamics for the particle system.

Finally, the entropy production of the $i$th particle given by equation (69) can be rewritten as:

$$ m_i T_i \frac{ds_i}{dt} = \frac{K}{2} \sum_j \mu_{ij} (\pi_{ij} r_{ij})^2 F_{ij} V_i V_j + m_i \frac{dq_i}{dt}. \quad (86) $$

### 6.3.1. Heat conduction and viscous forces modeling in presence of an interface.
In the present subsection formula (79) is examined in the case of the interaction of particles belonging to different fluids. In [47] SPH equations are developed for heat conduction problems where materials with different

thermal conductivities are present and, to tackle this discontinuity, a harmonic mean is used in the heat equation, i.e.:

$$ m_i \frac{dq_i}{dt} = 2 \sum_j k_{ij} (T_j - T_i) F_{ij} V_j V_i, \quad k_{ij} := \frac{2 k_i k_j}{k_i + k_j}, \quad (87) $$

where $k_i$ is the thermal conductivity of the generic particle $i$.

In the right-hand side of equation (87) the approximation of the Laplacian operator (39) can be recognized. In [47] heat conduction problems with discontinuous $k$, and with $k$ varying with $T$, were accurately solved. The structure presented in equation (87) was then extended to other diffusive processes. For example, in [187] the same equation is used for the time variation of salt concentration in water and for the viscosity diffusion process where equation (79) is reshaped as:

$$ m_i \mathbf{f}_i^V = -K \sum_j \mu_{ij} \pi_{ij} F_{ij} \mathbf{r}_{ij} V_i V_j, \quad \mu_{ij} := \frac{2 \mu_i \mu_j}{\mu_i + \mu_j}. \quad (88) $$

This equation is also used in [103] where it is shown that using the arithmetic mean between the viscosity of two different fluids, $\mathcal{X}$ and $\mathcal{Y}$, implies that the $\mathcal{X}$-particles close to the interface are strongly influenced by the $\mathcal{Y}$-particles, because of the SPH interpolation. Conversely, using the harmonic mean, the $\mathcal{X}$ particles close to the interface are comparatively less influenced by the $\mathcal{Y}$ particles, since the main contributions are those of the neighbor particles belonging to fluid $\mathcal{X}$. Concluding, the use of the harmonic mean formulation allows a better global accuracy with respect to using an arithmetic mean. In [103] a spatial convergence close to second order is obtained when solving a Poiseuille flow with two fluids having a viscosity ratio $\mu_{\mathcal{X}}/\mu_{\mathcal{Y}} = 100$.

It is worth noting that the use of the harmonic mean in the diffusive processes smooths the discontinuities of the fields $k$ and $\mu$ at the interfaces $\partial\Omega_{\mathcal{X}\mathcal{Y}}$, differently from the density field. Besides, the formulas (87) and (88) are based on symmetric exchanges of thermal energy and momentum, hence, the global conservation of the particle system is guaranteed.
## 6.4. SPH multiphase formulations

Considering the above subsections one can finally write the complete scheme for the SPH model. Different schemes can be

written considering the four choices discussed in section 6.2. As an example we consider here the combination which leads to an SPH scheme able to treat both the presence of multiple phases as well as of a free surface:

$$ 
\left\{
\begin{aligned}
\mathbf{S}_4 \text{ : SPH scheme for multiphase and free-surface flows} \\
\frac{dV_i}{dt} &= V_i^2 \sum_j (\mathbf{u}_{ij} \cdot \mathbf{r}_{ij}) F_{ij}, \quad \rho_i(t) = \frac{m_i}{V_i(t)}, \quad p_i = p(\rho_i, e_i), \\
m_i \frac{d\mathbf{u}_i}{dt} &= \sum_j (p_i V_i^2 + p_j V_j^2) \mathbf{r}_{ij} F_{ij} - K \sum_j \mu_{ij} \pi_{ij} F_{ij} \mathbf{r}_{ij} V_i V_j + m_i (\mathbf{f}_i + \mathbf{f}_i^\sigma + \mathbf{f}_i^S), \\
m_i \frac{de_i}{dt} &= -p_i V_i^2 \sum_j (\mathbf{u}_{ij} \cdot \mathbf{r}_{ij}) F_{ij} + \frac{K}{2} \sum_j \mu_{ij} (\pi_{ij} r_{ij})^2 F_{ij} V_i V_j + 2 \sum_j k_{ij} (T_j - T_i) F_{ij} V_i V_j, \\
\frac{d\mathbf{r}_i}{dt} &= \mathbf{u}_i.
\end{aligned}
\right.
\quad (89) $$

In order to derive the third equation of (89) the first law of thermodynamics (68) has been reshaped as:

$$ m_i \frac{de_i}{dt} = -p_i \frac{dV_i}{dt} + m_i T_i \frac{ds_i}{dt}, \quad (90) $$

where the time derivative $d\rho_i/dt$ is rewritten in terms of time derivative of the particle volume $V_i$. In this way, the continuity equation, i.e. the first equation of (89), can be used for the first term on the right-hand side of equation (90), while for the second term the entropy production equation (86) is adopted.

With the last equation of (89) the particles are advected in the fluid domain $\Omega$ by their material velocity.

In order to integrate in time the system of equations (89) the initial conditions on particle volumes, positions, velocities and internal energy have to be prescribed. The mass of the generic $i$th particle, $m_i$, is evaluated by the initial condition at time $t = t_0$:

$$ m_i = \rho_{i0} V_{i0}, \quad (91) $$

where $\rho_{i0}$ is the initial density field evaluated from the given initial pressure field using the EOS of the specific fluid, whereas the particle volumes $V_{i0}$ are known from the initial particle positions inside the fluid domain. Usually, at $t = t_0$ the particles are regularly positioned in the domain. Therefore, at the initial instant, the particle volumes coincide with the geometrical ones, and can be initialized as such. One can also use the particle packing technique described in [51]; in that case the volume can be assigned from the obtained irregular tessellation. Differently from the density and volume of the particles, the mass $m_i$ remains constant during the time evolution.

It is worth stressing that in all the summations of the scheme (89) particle $j$ can belong to a different fluid than that of particle $i$. As discussed in section 6.3.1, the discontinuities

of the viscosity and thermal conductivity fields are handled through the harmonic means between the particles. Only the discontinuity of the density field is explicitly treated; once the particle densities are prescribed in the initial condition, they are then advected in time through the particle motions.

## 6.5. Conservation properties of the SPH scheme

Regarding the SPH schemes derived in the previous subsections, such as the scheme $\mathbf{S}_4$ of equations (89), different conservation properties are guaranteed. First, since the particle masses $m_i$ are set in the initial condition and do not change during the time evolution, mass conservation is ensured.

Second, the pressure and viscous forces $f^P$ and $f^V$ are given by summation of components collinear with $\mathbf{r}_{ij}$, with all $ij$ terms that are symmetric with respect to swapping the indices; hence, the conservation of linear and angular momenta is also guaranteed. The above result is also predicted by the Noether theorem for which any symmetry in the Lagrangian results in the conservation of a physical quantity, see, e.g. [261, 135]. However, note that these conservative formulas for the pressure forces $f^P$ are not invariant by change of pressure reference; this drawback and its link with the momenta conservation were first pointed out in [71, 72].

Third, as already stressed, the first and second laws of thermodynamics are respected by the particle system. Finally, the equations are Galilean invariant and they do not change when passing from an inertial frame of reference to a non-inertial one. In the second case the non-inertial forces simply need to be added to the non-conservative force field $Q_i$, see, e.g. [52].

Beside these conservation properties, as highlighted in section 5, the conservation of the volume of the fluid domain $\Omega$ cannot be ensured. In other words, the consistency between particle positions and volumes given initially is lost during the
time evolution with errors which again depend on the two para-meters h and Δr, see section 10.5.

The second problem that can limit the derived SPH equations is that the energy conservation does not prevent non-physical transfers from mechanical energy to internal energy, as discussed in section 11.

## 6.6. Surface tension modeling from a microscopical description

The first general idea for modeling surface tension in SPH methods is to start from a mollification of the Dirac δ-function of the surface tension force (9) in the vicinity of the interface, Ω<sub>XY</sub>, so that it becomes a local volumic force (in absence of the Marangoni effect), i.e. for a particle i close to the interface:

$$ \rho f_i^\sigma = \sigma_{XY} \kappa_i n_i \delta_i^\sigma, \qquad (92) $$

where δ<sub>i</sub><sup>σ</sup> is a positive smooth approximation of the surface δ-function within a compact thickness around the interface. This thickness is typically proportional to the kernel radius, so that in the limit of h → 0 the consistency for f<sup>σ</sup> is recovered.

The second general idea is that the normal vector in equation (92) can be evaluated easily in this close-to-the-interface region by calculating the kernel gradient integral I<sub>3</sub> of equation (41) using only neighboring particles of the same fluid χ. More generally, this kernel gradient can be replaced by any approximate form of it and different variants are found in the related literature, e.g.:

$$ (a) \quad n_i = -\beta_i \sum_{j \in \chi} \nabla_i \tilde{W}_{ij} V_j, $$

$$ (c) \quad n_i = \beta_i \sum_{j \in \chi} r_{ij} \tilde{W}_{ij} V_j, $$

$$ (b) \quad n_i = \beta_i \sum_{j \in \chi} r_{ij} \tilde{W}_{ij}, $$

$$ (d) \quad n_i = \beta_i \sum_{j \in \chi} r_{ij} \tilde{W}_{ij} m_j, \qquad (93) $$

where β<sub>i</sub> is simply the normalizing factor for n<sub>i</sub> to be a unit vector, and ˜W<sub>ij</sub> is a compact kernel which has different shapes depending on authors. Note that β<sub>i</sub> has different physical dimensions in the different formulas. As discussed later, it is worth noting that, differently to the kernel W<sub>ij</sub>, the kernel ˜W<sub>ij</sub> can also have negative values on its support.

From this starting point which can be found in the pioneering work by [117], itself inspired from [36], two ways of modeling surface tension then arose. The straightforward way is to impose the physical value of σ<sub>XY</sub> and to find a way to calculate the curvature κ<sub>i</sub>, see section 7.8. However, since the latter calculation is not easy to perform on a set of particles, other authors drew inspiration from the microscopic generation of surface tension, trying to use the microscopic interaction principle at the macroscopic scale, and then calibrating this force so that its local macroscopic magnitude σ<sub>XY</sub> κ<sub>i</sub>, which depends on the scale and the nature of the problem solved, is recovered. These works inspired by the microscopic nature of surface tension are reported in this section.

For single-phase flows with a free surface, the next important work after [117] was [202] where the authors showed that using a Van der Waals fluid, a surface tension force directly comes from the pressure gradient of the attraction part of the Van der Waals EOS. As already noted in [117] however, these authors needed an interface thickness of twice the kernel radius, which increased the computational cost. To overcome this difficulty, an additional cohesion/repulsion force acting within the kernel radius was introduced in [248]:

$$ \begin{gathered} \rho f_i^\sigma = \sum_j s_{ij} r_{ij} \tilde{W}_{ij}, \\ \tilde{W}_{ij} = \begin{cases} \frac{1}{r_{ij}} \cos\left(\frac{\pi}{2} \frac{r_{ij}}{h}\right) & \text{for } r_{ij} \le h, \\ 0 & \text{for } r_{ij} > h \end{cases}, \quad (94) \end{gathered} $$

where s<sub>ij</sub> is the strength of the force acting between particles i and j.

Surface tension then arises from both the EOS pressure gradient and from this 'particle-particle' interaction force later on called 'pairwise force'. Following this path other authors, such as in [4, 6, 21, 133], then used different shapes for this interaction force (cosine, kernel gradients, polynomials...) combining it or not with a pressure-gradient part, in both weakly-compressible and incompressible SPH contexts. In all these works, the numerical surface tension strength has to be calibrated against known analytical solutions (liquid drop oscillations, Laplace-Young law, capillary rise in a tube, etc) to recover its physical value. Through these pairwise force models, authors are also able to model adhesion with contact angles through the use of wall boundary particles and calibrated force exchanged with the latter. Finally, it can be noted that some authors, such as in [54], used non-physically calibrated surface tension through the EOS pressure gradient so as to keep smooth interfaces in multiphase SPH simulations, using a summation per phase.

Recently, a significant step forward was made in [246] where, by using molecular dynamics techniques, the authors were able to analytically relate, for a given pairwise force shape, the imposed coefficients of this force to the physical surface tension and static contact angle values, with convincing results on the cases presented. Another major change with respect to previous literature is that this 'self-calibrated' pairwise force model is derived in the context of multiphase flows this time, not in that of single-phase free-surface flows. Subsequently, this proposed model was used in the context of single-phase free-surface incompressible SPH (iSPH) in [197].

## 6.7. Link between the SPH equations obtained and the NSEs

The system of equations S<sub>4</sub> (89), obtained starting from the equations governing a thermo-mechanical system of particles, is clearly a discretization of the NSEs (1). In the next section
the opposite path is followed, starting from the NSEs (1) and discretizing directly their differential operators.

## 7. Derivation of SPH equations: from local NSEs

In this section the NSEs (1) are rewritten, substituting the differential operators with their smoothed (mollified) version:

$$ 
\left\{ 
\begin{aligned}
\frac{d\rho}{dt} &= -\rho \langle \nabla \cdot \mathbf{u} \rangle, \\
\rho \frac{de}{dt} &= -p \langle \nabla \cdot \mathbf{u} \rangle + \langle \mathbf{V} : \mathbf{D} \rangle - \langle \nabla \cdot \mathbf{q} \rangle,
\end{aligned}
\right. 
\qquad \rho \frac{d\mathbf{u}}{dt} = -\langle \nabla p \rangle + \langle \nabla \cdot \mathbf{V} \rangle + \rho \mathbf{f}, \\
\frac{d\mathbf{r}}{dt} = \mathbf{u}, \quad p = p(\rho, e). \tag{95}
$$

The above equations are the non-local model of the NSEs, see [76] for more details, in which the smoothed differential operators need to be properly formulated.

The initial target is to derive the smoothed operators $\langle \nabla \cdot \mathbf{u} \rangle$ and $\langle \nabla p \rangle$. To this end, we start by considering an inviscid fluid, obeying the so-called Euler equations, excluding the viscous tensor $\mathbf{V}$ and the heat fluxes $\mathbf{q}$. The smoothed velocity gradient operator and its transpose can be obtained using equation (37) and applying it to the vector field $\mathbf{u}$:

$$ 
\left\{ 
\begin{aligned}
\langle \nabla \mathbf{u} \rangle (\mathbf{r}) &= \int_{\Omega} [\mathbf{u}(\mathbf{r}') - \mathbf{u}(\mathbf{r})] \otimes \nabla W \, d\mathbf{r}', \\
\langle \nabla \mathbf{u} \rangle^{\mathrm{T}} (\mathbf{r}) &= \int_{\Omega} \nabla W \otimes [\mathbf{u}(\mathbf{r}') - \mathbf{u}(\mathbf{r})] \, d\mathbf{r}',
\end{aligned}
\right. 
\qquad \rho \frac{de}{dt} = -p \langle \nabla \cdot \mathbf{u} \rangle + \langle \nabla p \rangle \cdot \mathbf{u}, \quad e_{\mathrm{tot}} := \left( e + \frac{u^2}{2} + U \right), \tag{96}
$$

and for the smoothed velocity divergence this gives:

$$ \langle \nabla \cdot \mathbf{u} \rangle (\mathbf{r}) = \int_{\Omega} [\mathbf{u}(\mathbf{r}') - \mathbf{u}(\mathbf{r})] \cdot \nabla W d\mathbf{r}'. \quad (97) $$

In the following subsection the consistency of the above formula is discussed when a free surface $\partial\Omega_F$ is present in the domain boundary $\partial\Omega$.

### 7.1. Consistency of the continuity equation in presence of free surface

Because of the dynamic boundary condition the pressure is constant along $\partial\Omega_F$. Since the fluids are considered barotropic with polytropic laws, see section 2, the density is then also constant along $\partial\Omega_F$. As a consequence, the velocity divergence is zero along the free surface, therefore we are in the same situation as discussed in section 3.2.1, i.e. formula (97) is consistent even for a point belonging to the free surface (see [49] for more details):

$$ 
\begin{aligned}
\langle \nabla \cdot \mathbf{u} \rangle (\mathbf{r}) &= \nabla \cdot \mathbf{u}(\mathbf{r}) + \mathcal{O}(h^2), && \forall \mathbf{r} \in (\Omega \setminus \partial\Omega), \\
\langle \nabla \cdot \mathbf{u} \rangle (\mathbf{r}) &= 0 + \mathcal{O}(h), && \forall \mathbf{r} \in \partial\Omega.
\end{aligned} 
\qquad (98) 
$$

Note that equation (97) is Galilean invariant in all $\Omega$ including its boundary $\partial\Omega$. In section 5.3 it has been shown that this property is a crucial aspect when considering the SPH particle system.

In the next subsection it is shown that once the smoothed operator $\langle \nabla \cdot \mathbf{u} \rangle$ has been selected its counterpart $\langle \nabla p \rangle$ can be derived in such a way as to preserve the conservation of momenta and energy.

### 7.2. Energy and momenta conservation

Combining the equation for the specific internal energy $e$ and the momentum equation multiplied by the velocity $\mathbf{u}$ in the previous equations (95) for an inviscid flow with no heat flux, the variation of the specific total energy $e_{\mathrm{tot}}$ can be derived as:

$$ \rho \frac{d e_{\mathrm{tot}}}{dt} = -p \langle \nabla \cdot \mathbf{u} \rangle + \langle \nabla p \rangle \cdot \mathbf{u}, \quad e_{\mathrm{tot}} := \left( e + \frac{u^2}{2} + U \right), \tag{99} $$

where only conservative force fields $f$ are considered, i.e. $f = -\nabla U$ where $U$ is the specific energy potential which is assumed not to vary over time. Integrating equation (99) over the domain $\Omega$ and applying the Leibniz–Reynolds transport theorem (15) to the left-hand side, the variation of the total energy $\mathcal{E}$ is given by:

$$ \frac{d\mathcal{E}}{dt} = \frac{d}{dt} \int_{\Omega} \rho e_{\mathrm{tot}} \, d\mathbf{r} = \int_{\Omega} -p \langle \nabla \cdot \mathbf{u} \rangle \, d\mathbf{r} + \int_{\Omega} \langle \nabla p \rangle \cdot \mathbf{u} \, d\mathbf{r}. \tag{100} $$

If we consider the case where the fluid domain $\Omega$ is confined by a free surface and by fixed solid surfaces, the fluid system results to be isolated and the total energy of the system must be conserved. It follows that:

$$ \int_{\Omega} -p \langle \nabla \cdot \mathbf{u} \rangle \, d\mathbf{r} + \int_{\Omega} \langle \nabla p \rangle \cdot \mathbf{u} \, d\mathbf{r} = 0. \tag{101} $$

In [49] it is shown that, thanks to the kernel symmetry given in equation (23), when $\langle \nabla \cdot \mathbf{u} \rangle$ is evaluated by the interpolation formula (97), the above equation (101) yields:

$$ \langle \nabla p \rangle (\mathbf{r}) = \int_{\Omega} [p(\mathbf{r}') + p(\mathbf{r})] \nabla W d\mathbf{r}'. \tag{102} $$

This is a key point of the SPH method: to conserve the energy the smoothed operator for the pressure gradient cannot be defined with the same interpolation formula (37) as for the
velocity divergence; the minus sign needs to be replaced by a plus sign.

With that form (102) of $\langle \nabla p \rangle$ the conservation of the momenta is also guaranteed thanks to the equalities:

$$ \left\{ \begin{array}{l} \displaystyle \int_{\Omega} \int_{\Omega} [p(\mathbf{r}') + p(\mathbf{r})] \nabla W d\mathbf{r}' d\mathbf{r} = 0, \\[1.5em] \displaystyle \int_{\Omega} \int_{\Omega} (\mathbf{r} - \mathbf{r}_0) \times [[p(\mathbf{r}') + p(\mathbf{r})] \nabla W] d\mathbf{r}' d\mathbf{r} = 0. \end{array} \right. \qquad (103) $$

### 7.3. Viscous force modeling

In order to derive a smoothed version of the viscous stress tensor:

$$ \langle \mathbf{V} \rangle (\mathbf{r}) = \lambda \langle \nabla \cdot \mathbf{u} \rangle (\mathbf{r}) + 2\mu \langle \mathbf{D} \rangle (\mathbf{r}), \quad (104) $$

the rate-of-strain tensor $\mathbf{D}$ needs to be mollified. To this purpose the smoothed velocity gradient operator (96) cannot be used since the resulting $\langle \mathbf{D} \rangle$ does not vanish for rigid body motion when a free surface is present in $\partial \Omega$. Indeed, for a point for which $\partial \Omega_W$ intersects $\partial \Omega$ the property $I_4$ of equations (25) is no longer true. To overcome this problem the corrected kernel (38) can be used, see [29]:

$$ \begin{aligned} \langle \nabla \mathbf{u} \rangle^L (\mathbf{r}) &= \left[ \int_{\Omega} [\mathbf{u}(\mathbf{r}') - \mathbf{u}(\mathbf{r})] \otimes \nabla W d\mathbf{r}' \right] \cdot \mathbf{L}(\mathbf{r}) \\ &= \langle \nabla \mathbf{u} \rangle (\mathbf{r}) \cdot \mathbf{L}(\mathbf{r}), \end{aligned} \qquad (105) $$

where the identity:

$$ \mathbf{a} \otimes (\mathbf{A} \cdot \mathbf{b}) = (\mathbf{a} \otimes \mathbf{b}) \cdot \mathbf{A}^{\mathrm{T}} \qquad (106) $$

and the fact that the tensor $\mathbf{L}$ is symmetric were used.

It follows that the rate-of-strain tensor $\mathbf{D}$ and the spin tensor $\mathbf{W}$ can be mollified as:

$$ \begin{aligned} \langle \mathbf{D} \rangle &= \frac{1}{2} \left[ \langle \nabla \mathbf{u} \rangle \cdot \mathbf{L} + \mathbf{L} \cdot \langle \nabla \mathbf{u} \rangle^{\mathrm{T}} \right], \\ \langle \mathbf{W} \rangle &= \frac{1}{2} \left[ \langle \nabla \mathbf{u} \rangle \cdot \mathbf{L} - \mathbf{L} \cdot \langle \nabla \mathbf{u} \rangle^{\mathrm{T}} \right]. \end{aligned} \qquad (107) $$

It is worth noting that from the spin tensor $\langle \mathbf{W} \rangle$ the vorticity $\langle \nabla \times \mathbf{u} \rangle^L$ can be extracted, see, e.g. [11].

As discussed in [29], in an analogous way to the pressure forces, when evaluating the divergence of the viscous stress tensor the following formula needs to be used in order to preserve momenta and energy conservation:

$$ \langle \nabla \cdot \mathbf{V} \rangle (\mathbf{r}) = \int_{\Omega} [\langle \mathbf{V} \rangle (\mathbf{r}') + \langle \mathbf{V} \rangle (\mathbf{r})] \cdot \nabla W d\mathbf{r}'. \quad (108) $$

Note that, as explained in [50], the evaluation of $\langle \nabla \cdot \mathbf{V} \rangle$ through the above relation is CPU demanding since it requires first evaluating $\langle \mathbf{V} \rangle$. To avoid this double step a direct evaluation of $\langle \nabla \cdot \mathbf{V} \rangle$ is usually adopted by the SPH community:

$$ \langle \nabla \cdot \mathbf{V} \rangle (\mathbf{r}) = -\mu K \int_{\Omega} \pi(\mathbf{r}, \mathbf{r}') (\mathbf{r} - \mathbf{r}') F d\mathbf{r}', $$

$$ \pi(\mathbf{r}, \mathbf{r}') := \frac{[\mathbf{u}(\mathbf{r}') - \mathbf{u}(\mathbf{r})] \cdot (\mathbf{r}' - \mathbf{r})}{|\mathbf{r}' - \mathbf{r}|^2}. \qquad (109) $$

This relation is simply the limit of equation (88) for $\Delta \mathbf{r}/h \to 0$ with the viscosity coefficient $\mu$ considered constant in $\Omega$. In [81] it is demonstrated that:

$$ \begin{aligned} \langle \nabla \cdot \mathbf{V} \rangle (\mathbf{r}) &= \mu \nabla^2 \mathbf{u}(\mathbf{r}) + 2\mu \nabla (\nabla \cdot \mathbf{u})(\mathbf{r}) + \mathcal{O}(h^2), \\ \forall \mathbf{r} \in (\Omega \setminus \partial \Omega_{\mathrm{F}}). \end{aligned} \qquad (110) $$

As highlighted in [53], equation (109) implies a constraint on the secondary viscosity, $\lambda = \mu$, or in other words, the bulk viscosity $\kappa = \lambda + 2\mu/n$ is equal to $\kappa = 2\mu$ in 2D and $\kappa = 5/3\mu$ in 3D. Therefore it is not possible to prescribe a zero bulk viscosity. On the other hand, when considering fluids with limited compressibility, such as liquids, $\lambda/\mu = \mathcal{O}(1)$ and, in addition, the effect of the bulk viscosity on the flow evolution is always negligible, see, e.g. [58]. As shown in [17] it is possible to add the $\langle \nabla(\nabla \cdot \mathbf{u}) \rangle$ term in order to tune the right value of $\lambda$ for gas dynamics problems where the secondary viscosity plays a relevant role. However, this topic is beyond the scope of the present work. Even with that drawback of the constraint on the secondary viscosity $\lambda$, formula (109) has the great advantage of preserving momenta and energy, so that it is widely used within the SPH community.

If incompressible flow is considered, another way to evaluate the viscous forces is to directly smooth the velocity Laplacian in equation (110) using equation (33):

$$ \mu \langle \nabla^2 \mathbf{u} \rangle (\mathbf{r}) = 2\mu \int_{\Omega} [\mathbf{u}(\mathbf{r}') - \mathbf{u}(\mathbf{r})] F d\mathbf{r}'. \qquad (111) $$

However, the above formula does not preserve the angular momentum and, following the analysis of [50], when a free surface is present the induced viscous dissipation is different from zero for a pure rigid body motion.

The total energy variation given in equation (100) becomes, in presence of viscous forces:

$$ \begin{aligned} \frac{d\mathcal{E}}{dt} = & \int_{\Omega} -p \langle \nabla \cdot \mathbf{u} \rangle d\mathbf{r} + \int_{\Omega} \langle \nabla p \rangle \cdot \mathbf{u} d\mathbf{r} + \int_{\Omega} \langle \mathbf{V} : \mathbf{D} \rangle d\mathbf{r} \\ & + \int_{\Omega} \langle \nabla \cdot \mathbf{V} \rangle \cdot \mathbf{u} d\mathbf{r}. \end{aligned} \qquad (112) $$

In [50] it is demonstrated that using formula (109), even when a free surface is present, the following relation holds:

$$ \int_{\Omega} \langle \nabla \cdot \mathbf{V} \rangle \cdot \mathbf{u} d\mathbf{r} = - \int_{\Omega} \langle \mathbf{V} : \mathbf{D} \rangle d\mathbf{r} + \mathcal{O}(h). \qquad (113) $$

This is an important relation that demonstrates that equation (109) is consistent in the limit of $h \to 0$, giving the right viscous dissipation expected from the governing equations (95) even in presence of a free surface.

In [58] the conclusions drawn above are used to validate the SPH model in the prediction of the viscous attenuation of gravity waves.
### 7.4. Enforcement of the dynamic free-surface boundary condition

Still considering the case where the fluid domain $\Omega$ is confined by a free surface and by fixed solid surfaces, using the limit $h \to 0$ and the divergence theorem, equation (112) can be reshaped as:

$$ \frac{d\mathcal{E}}{dt} = \int_{\partial\Omega_F} [-p(\mathbf{u} \cdot \mathbf{n}) + (\nabla \cdot \mathbf{u}) \cdot \mathbf{n}] d\mathbf{r} = 0. \quad (114) $$

In the above equation the surface integral applies only to the free surface $\partial\Omega_F$ since on the fixed solid boundaries $\partial\Omega_B$ the impermeability condition $\mathbf{u} \cdot \mathbf{n} = 0$ holds.

Equation (114) implies that the pressure and viscous forces cannot do any mechanical work on $\partial\Omega_F$. In other words, this means that the dynamic free-surface boundary condition, $\mathbf{T} \cdot \mathbf{n} = 0$ on $\partial\Omega_F$ (without surface tension), is always satisfied in an integral sense. This important result justifies why in the SPH method there is no need to enforce the dynamic boundary condition on the free surface in an explicit way. This property is one of the most relevant that led to the continuously growing interest in the method. Indeed, when simulating complex 3D free-surface flow with large deformation of the free surface and complex fragmentation, it is not necessary to identify the particles that belong to the free surface in order to impose pressure values on them. It should be remembered that in a

Lagrangian framework the kinematic boundary condition is also intrinsically verified, at convergence, so that there is no need to impose it either.

### 7.5. Heat diffusion modeling

Using the Fourier law $\mathbf{q} = -k\nabla T$, the diffusion of the heat fluxes $\langle \nabla \cdot \mathbf{q} \rangle$ can be formulated as:

$$ \langle \nabla \cdot \mathbf{q} \rangle (\mathbf{r}) = k \langle \nabla^2 T \rangle (\mathbf{r}) = 2k \int_{\Omega} [T(\mathbf{r}') - T(\mathbf{r})] F d\mathbf{r}', \quad (115) $$

where the smoothed Laplacian operator (33) has been used. If one wants to prescribe that $\partial\Omega_F$ is also an isothermal surface, differently from what is discussed above for the dynamic boundary condition, it is then necessary to identify $\partial\Omega_F$ and explicitly enforce the Dirichlet boundary condition on it. In [253] an extension to thermo-capillary flows is proposed.

### 7.6. Strong formulation of the SPH equations

Using the quadrature formula presented in section 3.3 the smoothed differential operators presented in the previous subsections can be discretized at the fluid particle locations and a fifth SPH scheme is derived as:

$$ \left\{ \begin{aligned} \mathbf{S}_5 &\text{: SPH scheme for multiphase and free-surface flows} \\ \frac{d\rho_i}{dt} &= -\rho_i \sum_j (\mathbf{u}_{ij} \cdot \mathbf{r}_{ij}) F_{ij} V_j, \quad V_i(t) = \frac{m_i}{\rho_i(t)}, \quad p_i = p(\rho_i, e_i), \\ \rho_i \frac{d\mathbf{u}_i}{dt} &= \sum_j (p_i + p_j) \mathbf{r}_{ij} F_{ij} V_j - K \sum_j \mu_{ij} \pi_{ij} F_{ij} \mathbf{r}_{ij} V_j + \rho_i (\mathbf{f}_i + \mathbf{f}_i^\sigma + \mathbf{f}_i^S), \\ \rho_i \frac{de_i}{dt} &= -p_i \sum_j (\mathbf{u}_{ij} \cdot \mathbf{r}_{ij}) F_{ij} V_j + \frac{K}{2} \sum_j \mu_{ij} (\pi_{ij} \mathbf{r}_{ij})^2 F_{ij} V_j + 2 \sum_j k_{ij} (T_j - T_i) F_{ij} V_j, \\ \frac{d\mathbf{r}_i}{dt} &= \mathbf{u}_i, \end{aligned} \right. \quad (116) $$

where the term $\langle \mathbf{V} : \mathbf{D} \rangle$ in equations (95) has been discretized from equation (84) where $\Phi_D$ is the integral of $\langle \mathbf{V} : \mathbf{D} \rangle$ over the domain $\Omega$.

Without loss of generality, so as to be able to compare the different schemes, the forces $f_i^S$ related to the constraints on solid surfaces are added in the above system of equations, where the harmonic means of the viscosity and of the thermal conductivity are also inserted in the summations, as done in the previous sections.

### 7.7. SPH multiphase formulation

The scheme $\mathbf{S}_5$ was derived directly by discretizing the governing equations; as a result this scheme is applicable to multiphase and free-surface flows. This derivation allowed the discussion of the boundary conditions on the free-surface boundary.

This scheme $\mathbf{S}_5$ presents, as expected, strong similarities with $\mathbf{S}_4$, equations (89). The small differences between the two
schemes directly result from the different ways these two models are derived.

This scheme $S_5$ was derived in [54] following the previous work of [29] and was applied successfully to gas-liquid flows, underlining that using equation (102) to evaluate the pressure forces makes it possible, as with the scheme $S_4$, to overcome the problems linked to the schemes $S_1$ and $S_3$ when fluids with large density ratios are considered.

## 7.8. Surface tension modeling from a continuous description

As mentioned in section 6.6, the very large majority of authors see the surface tension as a local volumic force acting in the vicinity of the interface, i.e. equation (92) recalled here (in absence of the Marangoni effect):

$$ \rho f_i^\sigma = \sigma_{XY} \kappa_i n_i \delta_i^\sigma = -\sigma_{XY} (\nabla_i \cdot n_i) n_i \delta_i^\sigma. \quad (117) $$

Contrary to the works described in section 6.6, those reported here directly evaluate the normal and curvature in equation (117), and impose the physical value of $\sigma_{XY}$. These studies can be separated into two main branches. The first branch is the use of the continuum surface force (CSF) method introduced in [32] for mesh-based methods. The second branch is the use of a reconstruction of the interface. At the end of the present section two other ways are also mentioned where equation (117) is not used.

### 7.8.1. CSF method.
The majority of authors have worked with the CSF method. Interestingly, unlike the microscopic-inspired methods of section 6.6 which were first derived for free-surface single-phase flows and were much later extended to multiphase flows, the opposite occurred with CSF methods in SPH. The principle of CSF as proposed in [32] is to use a color function $c$ identifying each phase, which enables a direct calculation of the normal vector as $n_i = (\nabla c)_i / |\nabla c|_i$. The curvature can then be calculated as $\kappa_i = -\nabla \cdot n_i$. In the derivation proposed in [32], it is also shown that $\delta_i^\sigma = |\nabla c|_i / [c]$ where $[c]$ is the color jump between the two phases, often taken equal to 1 or to the jump of density [$\rho$] as also proposed by these authors. Static wall adhesion through CSF is also developed in [32], by imposing the interface normal at the wall.

The first work that adapted CSF to SPH was [192] where the author used, to improve accuracy, symmetric forms of the color SPH-gradient and normal SPH-divergence, like the ones in equation (39), a threshold on the particles contributing to the normal divergence calculation and a Shepard kernel as in equation (45). This author also compared the results with the original CSF to the ones using the variant proposed in [140] where:

$$ \rho f_i^\sigma = \sigma_{XY} \nabla \cdot [(\mathbb{I} - n_i \otimes n_i) \delta_i^\sigma], \quad (118) $$

is used rather than equation (117); he found comparable results with the two variants. Later on, in [120] this second formalism

based on equation (118) was adopted with minor modifications and applied also to cases with adhesion and more than two phases. The same group proposed in [1] another scheme based on the use of the density as the color function, as proposed in [32], and the use of a simpler way to calculate the divergence of the normal.

The latter work was extended in [34] to adhesion test cases with implementation of a restoring adhesion force imposed smoothly in a region close to the wall, tending to the static equilibrium angles. Finally, in [122] a more rigorous way of treating the contact line was derived, where a contact line force was proposed, in analogy to the derivation of the CSF for the surface tension, based on the physical derivation in [109]; accurate validation was provided. Without wall adhesion, an interesting set of applications of these two-phase SPH CSF models can be found in [275], including comparison to experiments. In addition, an interesting discussion on the need to use an additional artificial surface tension model to control the smoothness of interfaces, see section 6.6, can be found in [244].

Recently, the CSF method was extended to single-phase free-surface flows simultaneously in the works [114, 207], with application to droplet collisions in both cases. In the former work [114], a comparison is made between the weakly-compressible and the incompressible formulations, with quite similar results obtained and in rather good agreement to experiments. Lately, more accurate versions for single-phase flows were proposed in [156, 258] where both normal and curvature evaluation accuracy in equation (117) were improved by higher-order SPH interpolation. In the latter work [258], specific treatments were also included especially for thin jets and contact with a solid surface, and a large set of test cases was provided with reference solutions, including the Plateau-Rayleigh instability (see figure 7).

It can be mentioned that in [249] the authors derived a CSF model including a Marangoni force, i.e. the gradient of surface tension coefficient along the interface in equation (9), due to temperature variations in their case. This work on a thermocapillary model was further extended and applied by a number of authors, e.g. [118, 155, 174].

### 7.8.2. Reconstruction of a continuous interface from the particle arrangement.
As mentioned above, in a second branch of works authors opted for a reconstruction of the interface so as to be able to use geometric evaluations of the normal and curvature in equation (117). All these schemes need an accurate boundary detection algorithm from which they next reconstruct an interface with different algorithms. The surface tension forces $f_i^\sigma$ are directly applied on the particles belonging to the interfaces, which means that $\delta_i^\sigma$ is a delta Dirac function on the $i$th particle.

In [279] a polynomial fit was used in 2D and an MLS reconstruction in 3D. The authors of [274] opted for a different strategy where the interface is meshed initially by marching cubes and then this mesh is updated according to particle
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Plateau–Rayleigh instability (inviscid fluid): propagation of initial perturbations at the extremities of the fluid cylinder leading to the fragmentation into several drops; the free-surface deformation period just before fragmentation is then successfully compared to the analytical solution in the article. Reprinted from [258], Copyright (2023), with permission from Elsevier.</figcaption></figure>

velocities near the interface. In a simpler version [149] used Delaunay triangulations of each phase interface particles at each time step, and average the surface tension acting on these two proximal surfaces.

### 7.8.3. Surface tension through jump condition at the interfaces.
Another way is to impose directly the pressure jump induced by surface tension in the boundary condition (8). This is done in [142] through the use of a Riemann solver, see section 11.2, where this pressure jump is imposed directly in the contact discontinuity of the solution linked to the presence of the interface, for each pair of particles of different phases. Another possibility to impose this pressure jump is to use the ghost fluid method (GFM), popular in mesh-based methods, as in [280]. In the latter work the authors use a complex and probably computationally expensive strategy where the ghost values are obtained through a combination of: (i) LS construction to represent the interface and calculate its curvature, (ii) velocity reconstruction in both phases close to the interface to impose the correct stress balance at the interface, (iii) extrapolation to the particle of the other phase used as the ghost particles of the GFM.

### 7.8.4. Surface tension forces through minimization of the interface energy.
The surface tension acts such that the interface energy is minimal. Following this idea the authors of [110] minimized a discrete smoothed energy density calculated for each particle in the vicinity of the interface to get a local surface tension force in a single-phase free-surface context. Interestingly, this idea had already been developed in a more consistent way by J J Monaghan in an unpublished report [186] where he: (i) defined an SPH approximation of the free-surface area based on a volumic integration of the smoothed delta-function, (ii) used this free-surface area to define the surface

energy, (iii) inserted the surface energy into the Lagrangian of the system, (iv) derived a local conservative surface tension force in the momentum equation. Recent works [37, 132] are following that same path.

## 8. Derivation of SPH equations: from continuous conservation laws

A sixth SPH scheme can be derived following the pioneer work by Vila [259]. This scheme was written starting directly from a generic conservation law for a generic vector $\Phi$ of conservative variables. The conservation law reads:

$$ \mathcal{L}_u(\Phi) + \nabla \cdot \mathbf{F}(\Phi) = \mathbf{S}, \qquad \mathcal{L}_u(\Phi) := \frac{\partial \Phi}{\partial t} + \nabla \cdot (\mathbf{u} \otimes \Phi), \qquad (119) $$

where $\mathbf{F}$ is the flux tensor of the conservative variables, $\mathbf{S}$ are the source terms and $\mathcal{L}_u$ is the transport operator. Equation (119) is the one adopted to derive Finite Volume schemes in conventional mesh-based methods. Differently here, we want to discretize the above equation on a set of moving particles whose positions and volumes change under the action of a continuous velocity field $\mathbf{u}(\mathbf{r}, t)$:

$$ \frac{d\mathbf{r}_i}{dt} = \mathbf{u}(\mathbf{r}_i, t), \qquad \frac{dV_i}{dt} = V_i \nabla \cdot \mathbf{u}(\mathbf{r}_i, t), \qquad (120) $$

using again the volumetric strain-rate equation (4). In the simplest case the particles are initially distributed on a regular Cartesian lattice with spacing $\Delta r$ and $V_i(t=0) = (\Delta r)^n$.

In order to discretize the local conservation law (119) on the set of particles, we start by integrating it over the fluid domain $\Omega$:

$$ \int_{\Omega} \left[ \frac{\partial \Phi}{\partial t} + \nabla \cdot (\mathbf{u} \otimes \Phi) + \nabla \cdot \mathbf{F}(\Phi) \right] d\mathbf{r} = \int_{\Omega} \mathbf{S} d\mathbf{r}. \qquad (121) $$
Applying the Leibniz-Reynolds transport theorem to the first two terms we get:

$$ \frac{d}{dt} \int_{\Omega} \mathbf{\Phi} \, d\mathbf{r} + \int_{\Omega} \nabla \cdot \mathbf{F}(\mathbf{\Phi}) \, d\mathbf{r} = \int_{\Omega} S \, d\mathbf{r}, \quad (122) $$

where the second term on the left-hand side is nothing other than the flux of tensor $\mathbf{F}$ across $\partial\Omega$ once the divergence theorem is adopted. When this flux and the source terms are zero the global quantities $\int_{\Omega} \mathbf{\Phi} \, d\mathbf{r}$ must be conserved.

Using a simple quadrature formula, equation (122) can be discretized on the set of particles obtaining:

$$ \sum_{i} \left[ \frac{d(\mathbf{\Phi}_i V_i)}{dt} + \langle \nabla \cdot \mathbf{F}(\mathbf{\Phi}) \rangle_i V_i \right] = \sum_{i} S_i V_i, \quad (123) $$

where the divergent operator has been replaced with a particle approximation, i.e. using a smoothed version of the divergence operator. We should remark here that in the definition of Vila, the particle approximation of a function $f$ is its discrete convolution with the delta Dirac function over the set of particles:

$$ \Pi(f)(\mathbf{r}) = \sum_{i} V_i f_i \delta(\mathbf{r} - \mathbf{r}_i), \quad (124) $$

which is then regularized (or smoothed) to get the smoothed particle approximation $\langle f \rangle_i$, i.e. the other way round with respect to what was developed in section 3.

If we consider an isolated fluid domain delimited just by a free surface, i.e. $\partial\Omega = \partial\Omega_F$ the flux of the tensor $\mathbf{F}$ across the boundary has to be zero. In such a condition in [259] it is underlined that, in order to ensure a null flux across $\partial\Omega_F$ at the discrete level, the smoothed divergence operator needs to be in a symmetric form, i.e.:

$$ \sum_{i} \langle \nabla \cdot \mathbf{F}(\mathbf{\Phi}) \rangle_i V_i = \sum_{i} \sum_{j} [\mathbb{F}_i(\mathbf{\Phi}) + \mathbb{F}_j(\mathbf{\Phi})] \cdot \nabla W_{ij} V_j V_i = 0. \quad (125) $$

Indeed, thanks to the symmetry property (23) of the kernel the fluxes cancel each other out in pairs. It follows that for each particle the governing equations can be written as:

$$ \left\{ \begin{array}{l} \displaystyle \frac{d\mathbf{r}_i}{dt} = \mathbf{u}(\mathbf{r}_i, t), \qquad \frac{dV_i}{dt} = V_i \nabla \cdot \mathbf{u}(\mathbf{r}_i, t), \\[1.5em] \displaystyle \frac{d(\mathbf{\Phi}_i V_i)}{dt} + V_i \sum_j [\mathbb{F}_i(\mathbf{\Phi}) + \mathbb{F}_j(\mathbf{\Phi})] \cdot \nabla W_{ij} V_j = S_i V_i. \end{array} \right. \quad (126) $$

We can now apply the above conservation law to the Euler equations for a compressible inviscid fluid. In such a case the vector $\mathbf{\Phi}$ of conservative variables and the tensor

$\mathbb{F}(\mathbf{\Phi})$ of the related fluxes are:

$$ \begin{aligned} \mathbf{\Phi} : &= \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ \rho w \\ \rho e_{\text{tot}} \\ 0 \\ \rho f_x \\ \rho f_y \\ \rho f_z \\ 0 \end{pmatrix}, & \mathbb{F}(\mathbf{\Phi}) : &= \begin{pmatrix} 0 & 0 & 0 \\ p & 0 & 0 \\ 0 & p & 0 \\ 0 & 0 & p \\ pu & pv & pw \end{pmatrix}, \\ \mathbf{S} : &= \begin{pmatrix} \rho e_{\text{tot}} \\ 0 \\ \rho f_x \\ \rho f_y \\ \rho f_z \\ 0 \end{pmatrix}, & & \end{aligned} \quad (127) $$

where $u, v, w$ are the components of the velocity vector $\mathbf{u}$ and $f_x, f_y, f_z$ the components of the conservative force field $\mathbf{f}$. Introducing the above definition in equations (126) the discrete form of the equations for the generic $i$th particle is:

$$ \left\{ \begin{aligned} \frac{d\mathbf{r}_i}{dt} &= \mathbf{u}_i, & \frac{dV_i}{dt} &= V_i \langle \nabla \cdot \mathbf{u} \rangle_i, \\ \frac{d(\rho_i V_i)}{dt} &= 0, & p_i &= f(e_i, \rho_i), \\ \frac{d(\rho_i \mathbf{u}_i V_i)}{dt} &= V_i \sum_j (p_i + p_j) \mathbf{r}_{ij} F_{ij} V_j + m_i \mathbf{f}_i, \\ \frac{d(\rho_i e_{\text{tot}} V_i)}{dt} &= V_i \sum_j (p_i \mathbf{u}_i + p_j \mathbf{u}_j) \cdot \mathbf{r}_{ij} F_{ij} V_j, \end{aligned} \right. \quad (128) $$

where the divergence of the velocity field is now also approximated through a smoothed operator, and the force field $\mathbf{f}$ enters through the source terms vector $\mathbf{S}$. The derivation of such equations from the above conservative laws ensures the conservation of mass, momenta and energy of the particle system.

As done in the previous sections, the mass of $i$th particle $m_i = \rho_i V_i$ can be fixed at the initial time and remains constant during time evolution, satisfying the third equation of (128). Multiplying the fifth equation of (128) by $\mathbf{u}_i$ and subtracting it from the sixth one, the equation for the internal specific energy $e$ can be derived and the system (128) can be reshaped as:

$$ \left\{ \begin{aligned} \frac{d\mathbf{r}_i}{dt} &= \mathbf{u}_i, & \frac{dV_i}{dt} &= V_i \langle \nabla \cdot \mathbf{u} \rangle_i, \\ \rho_i(t) &= \frac{m_i}{V_i(t)}, & p_i &= f(e_i, \rho_i), \\ \rho_i \frac{d\mathbf{u}_i}{dt} &= \sum_j (p_i + p_j) \mathbf{r}_{ij} F_{ij} V_j + \rho_i \mathbf{f}_i, \\ \rho_i \frac{d\mathbf{e}_i}{dt} &= -\sum_j p_j (\mathbf{u}_{ij} \cdot \mathbf{r}_{ij}) F_{ij} V_j. \end{aligned} \right. \quad (129) $$
At this level the system (129) is not closed since we still have not derived an expression for the smoothed velocity divergence $\langle \nabla \cdot \mathbf{u} \rangle_i$.

To this end, in analogy with the procedure presented in section 5.3, using the pressure distribution instead of the density distribution, it is possible to show that the right-hand side of the energy equation can be written as:

$$ \rho_i \frac{\mathrm{d}e_i}{\mathrm{d}t} = -\langle p \nabla \cdot \mathbf{u} \rangle_i, \quad \langle p \nabla \cdot \mathbf{u} \rangle_i = \sum_j p_j (\mathbf{u}_{ij} \cdot \mathbf{r}_{ij}) F_{ij} V_j. \qquad (130) $$

Hence, the velocity divergence can be approximated again as:

$$ \langle \nabla \cdot \mathbf{u} \rangle_i = \sum_j (\mathbf{u}_{ij} \cdot \mathbf{r}_{ij}) F_{ij} V_j, \qquad (131) $$

which can be used consistently in the volumetric strain-rate equation in the system (129).

Note that, summarizing, the latter procedure described in [259] to preserve conservation of the continuous conservation laws at the discrete level prescribed: (i) the form of the pressure force term in the momentum equation and, (ii) the velocity divergence form in the energy equation. The latter was then used in the particle volume evolution equation. Conversely, in section 7 above, the derivation started from the continuity equation with a choice of a velocity divergence form which then determined the form of the pressure force term to ensure discrete conservation of the energy.

### 8.1. SPH multiphase formulation

From the previous subsection, adding the viscous, surface-tension and heat transfer models presented in the previous sections, the sixth and last SPH scheme presented in this review article results as:

$$ \left\{ \begin{aligned} S_6 &: \text{SPH schemes for multiphase and free-surface flows} \\ \frac{\mathrm{d}V_i}{\mathrm{d}t} &= V_i \sum_j (\mathbf{u}_{ij} \cdot \mathbf{r}_{ij}) F_{ij} V_j, \quad \rho_i(t) = \frac{m_i}{V_i(t)}, \quad p_i = p(\rho_i, e_i), \\ \rho_i \frac{\mathrm{d}\mathbf{u}_i}{\mathrm{d}t} &= \sum_j (p_i + p_j) \mathbf{r}_{ij} F_{ij} V_j - K \sum_j \mu_{ij} \pi_{ij} F_{ij} \mathbf{r}_{ij} V_j + \rho_i (\mathbf{f}_i + \mathbf{f}_i^\sigma + \mathbf{f}_i^S), \\ \rho_i \frac{\mathrm{d}e_i}{\mathrm{d}t} &= -\sum_j p_j (\mathbf{u}_{ij} \cdot \mathbf{r}_{ij}) F_{ij} V_j + \frac{K}{2} \sum_j \mu_{ij} (\pi_{ij} \mathbf{r}_{ij})^2 F_{ij} V_j + 2 \sum_j k_{ij} (T_j - T_i) F_{ij} V_j, \\ \frac{\mathrm{d}\mathbf{r}_i}{\mathrm{d}t} &= \mathbf{u}_i. \end{aligned} \right. \qquad (132) $$

This last scheme is the one adopted in [102, 103] and more recently in [106].

In these articles free-surface multiphase simulations with large density ratios are considered and the SPH models are validated against reference solutions and experimental data available in the literature.

In [106] the scheme $S_6$ was adopted to simulate the water entry of a corrugated panel. In this problem an air pocket is confined between the panel corrugation and the water surface and it is compressed during the water entry until it starts to escape laterally, see figure 8. It is worth noting that for this problem in [106] a multiphase flow is considered (left of the corrugation) together with a free surface one (right of the corrugation). The regularity of the pressure contours and the continuity of the pressure across the interface are noticeable.

In [86] the same problem was studied for different inclinations of the corrugated panel, see figure 9. Conversely to [106], in [86] the SPH scheme $S_5$ was adopted.

## 9. Differences between the derived SPH models in sections 6–8

As we have seen, the SPH models $S_4$, $S_5$ and $S_6$ result from different visions of how to derive an SPH system of discrete equations for solving multiphase NSEs. Nonetheless, all of these derivations strongly rely on preserving conservation of the discrete system.

In the end, the differences in these three schemes $S_4$, $S_5$ and $S_6$ are only reflected in some different positions and indices of the particle volumes and pressures in the systems
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Water entry of a corrugated panel with trapped air cavity and a free surface. Top: sketch of the problem. Bottom: evolution of the pressure field for four time instants. Only the right half of the domain is shown, $x^* = x/L > 0$. The black lines display the air-water interface position. The trapped air pocket lies in the top left part of each plot, under the corrugation. Reprinted from [106], Copyright (2020), with permission from Elsevier.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Water entry of a corrugated panel with trapped air cavity. Evolution of the trapped air pocket for two time instants and three different angles of impact. Reprinted from [86], Copyright (2022), with permission from Elsevier.</figcaption></figure>

of equations. Nevertheless, as the pressure field is continuous (in the absence of surface tension), the difference between two neighbor evaluations $p_i$ and $p_j$ reduces at convergence when $i$ and $j$ get closer and closer. In turn, this is also true for the density and the volumes $V_i$ which are linked to the pressure.

Thus, these three schemes become equivalent at convergence, so that they should not lead to significant differences in results for reasonably converged resolutions. E.g. in [106], it is shown that schemes $S_4$ and $S_6$ give practically the same outputs on the benchmarks traditionally adopted for the validation of these SPH models.
# 10. Particle arrangement within their Lagrangian evolution

## 10.1. Lagrangian trajectories, particle arrangement and interpolation accuracy

When using Lagrangian numerical methods the motion of the computation nodes can induce different problems linked both to the stability and the accuracy of the scheme. For example in the FEM the motion of the grid nodes leads to the so-called hourglass instability which is essentially a spurious deformation mode of a mesh, typically manifesting as a patchwork of zig-zag or hourglass-like element shapes where individual elements are severely deformed while the overall mesh section is undeformed.

In a meshless context as in SPH, following exactly Lagrangian trajectories induces anisotropic particle distributions, especially in regions where the flow converges, diverges or bifurcates, tending to aggregate particles in the flow direction and to take them away in the transverse direction, or the contrary (see example in figure 10).

This anisotropy negatively impacts the SPH interpolation quality. In fact, as stressed in section 3 the interpolation error is linked to the particle disorder which impacts the $I_i$-terms quadrature errors (41). A measure of these errors is given in [48, 215, 10]. In particular, in [10] an algorithm to measure the level of particle disorder is given and its connection with the $I_i$-terms is discussed. Furthermore, the above works show that if the level of particle disorder does not decrease with the spatial resolution, the SPH interpolation no longer converges.

A way to reduce the quadrature errors made in the initial condition was provided by [51] where a Packing algorithm was used to determine the particle positions in a fluid domain $\Omega$ confined by a boundary $\partial\Omega$ of arbitrary shape. The Packing algorithm uses the SPH kernel properties to minimize the $I_3$ error on the kernel gradient. In [10] it is proven in a heuristic way that the minimization of this $I_3$ quadrature error leads to a strong reduction of the particle disorder.

In [141] it is stressed that in fluid domain zones characterized by a positive pressure field the SPH pressure force formulas induce a self-reordering mechanism during the simulation which is proportional to the pressure intensity. Conversely, in regions of negative pressure this self-reordering becomes the source of the so-called 'tensile instability' which is addressed in section 11. For these two reasons, in SPH simulations of confined flow the addition of a background pressure generally improves the quality of the results and eliminates the occurrence of tensile instability. Thanks to this behavior the SPH solution then presents a convergence rate between 1 and 2, which would not be possible if the particles followed exactly their Lagrangian paths.

This dynamic self-reordering mechanism of SPH particles is a key feature of the method in its standard version. Although this mechanism itself originates from quadrature errors on the pressure gradient, it is such that particle disorder is lowered during the simulations, by locally breaking the pure Lagrangian paths; thus quadrature errors are lower and the scheme is, in turn, more accurate, convergent and stable.

A theoretical analysis of this process proves to be quite difficult to achieve since it involves the displacement of the particles during the time evolution, i.e. it involves the solution of the whole system of equations, and it is thus still an open problem.

However, in order to control and make more consistent the reduction of particle disorder and the consequent reduction of the $I_i$ quadrature errors during the time evolution of the particle system, specific algorithms have been derived in recent years. This topic is addressed in the following subsections.

## 10.2. Remeshing

A possible way to improve the accuracy of the SPH scheme was proposed in [41], involving a periodic remeshing of the particle set into a regular lattice. The latter drastically reduces the $I_i$ quadrature errors. Such a remeshing requires interpolating the quantities associated to the particles at their new regular positions. As stressed in [41] it is possible to choose interpolation functions which ensure conservation of mass and momenta, and special treatments are needed for the fluid particles close to the solid boundaries $\Omega_B$. This apparently advantageous remeshing procedure has however two main drawbacks. The first is that the required periodic interpolation of the particle data induces a level of numerical diffusion. The second drawback is that the remeshing procedure can significantly alter the interfaces between fluids, and free-surface configurations, as shown in [141]. While the first drawback can be alleviated through a proper choice of interpolation functions, the second drawback is more complicated to deal with. As a result, this technique is seldom used in practice.

## 10.3. Particle Shifting Technique (PST)

A second strategy for controlling the disorder of the particle distribution is the PST. The latter has become popular in the SPH community in the last decade. It consists in adding a dynamic modification $\hat{u}$ to the transport velocity $u$, or equivalently modifying directly the particle displacements, in order to continuously rearrange particles in uniform configurations, resulting in a great improvement in the accuracy and robustness of the SPH scheme.

The first PST was introduced by Monaghan in [182], which he afterward named XSPH in [183]. Despite good conservation and consistency properties, see [180], this XSPH shifting velocity has the major defect of not intrinsically pointing towards 'void spaces' between particles. Twenty years later, the current class of PST emerged, first proposed by [200] in the finite volume particle method (FVPM) context, relying on locally using $\sum_j \omega_{ij} r_{ij}$ as shifting velocity, which actually points towards the least dense region around particle $i$. $\omega_{ij}$ stands here for a given weight and in SPH variants most authors relate it to $F_{ij}$, so that they use rather forms as $\sum_j (\hat{\omega}_{ij} - h \nabla W_{ij} V_j)$ where $\hat{\omega}_{ij}$ is a given weight depending on authors, whose dimension is a velocity. This PST by [200] was then first extended to the Incompressible SPH method by [151, 269] and to the Weakly-Compressible SPH method by [3, 227].
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Viscous flow past a square cylinder at Reynolds number 50. Snapshots of passive Lagrangian tracers evolution. Reprinted from [205], Copyright (2016), with permission from Elsevier.</figcaption></figure>

Note that in the moving particle semi-implicit (MPS) context, authors worked rather on a pressure gradient correction based on the same shape, which in turn yields a particle shifting acceleration rather than velocity or displacement, with similar effects however, see, e.g. [250].

In [236], where the PST was adapted to simulations with different particle resolution zones, leading to the so-called $\delta^+$-SPH model, it is stressed that the PST introduced in a weakly-compressible scheme leads to non-conservation of momenta and energy. In [235] this problem of non-conservation of the $\delta^+$-SPH scheme is analyzed more in depth. In particular, following the work by [205] it is shown that a consistent PST can be derived by introducing a shifting velocity $\hat{u}$ instead of a shifting displacement of the particles.

Different laws are possible for the shifting velocity $\hat{u}$, see [180]. The choice is on how to define $\hat{\omega}_{ij}$ and on how to limit its magnitude in circumstances where it would result in too large shifting displacements within one time step. For example in [170], where violent sloshing flow is considered, the limiter involves the global maximum velocity of the flow and the local maximum velocity difference over the neighbors.

As remarked in [130], the shifting velocity close to interfaces has to be modified to be consistent with the kinematic boundary condition on the interface. In particular, the component of $\hat{u}$ normal to the interface has to be nullified while the tangential component is maintained unaltered. More details are also given in [235].

Recently in [180], a review and theoretical study of existing PSTs was presented, especially analyzing their consistency and conservation properties. A new law was also proposed based on a weight $\hat{\omega}_{ij}$ involving $\max_j(\pi_{ij})$ (see equation (77)). This shifting velocity law presents several benefits with respect to previous laws. In particular, it is consistent both in shifting displacement and velocity, it is Galilean and local rotation invariant, and it does not depend on global flow quantities. Furthermore, special attention is also given to correct the shifting velocity in the proximity of the free surface and of solid boundaries.

Finally, in [219, 220] the PST is performed in an implicit way. During the time evolution the algorithm controls the particle spatial distribution through an implicit minimization problem, related to the particle concentration gradient and therefore, to the particles' anisotropy (see figure 11). The latter procedure provides a shifting displacement field from which a shifting velocity $\hat{u}$ is evaluated and introduced in an ALE formalism as explained in the next subsection. Within this procedure, it is no longer necessary to introduce a limitation to the shifting velocity magnitude.

### 10.4. Consistent particle shifting in the ALE framework

As presented in section 2.3 the NSEs can be written in an Arbitrary-Lagrangian-Eulerian framework [259]. In those equations the fluid elements are advected with a modified
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Jet impinging on a flat plate. Left and middle plots: particle distributions without (left) and with PST (middle). Right plot: reduction of the maximum particle shifting velocity $\hat{u}$ in the domain when increasing the spatial resolution. (Left and Middle) Reproduced from [219]. CC BY 4.0. (Right) Reprinted from [180], Copyright (2022), with permission from Elsevier.</figcaption></figure>

velocity ($u + \hat{u}$), with $\hat{u}$ an arbitrary velocity field. Thanks to this framework the PST can be introduced in a consistent way at the discrete level, independently of the choice of $\hat{u}$ or of the fact that it tends to zero when refining since it preserves the conservation properties of the scheme, as discussed in [205]. The latter conservation benefit is also stressed in [16] where the ALE framework was used to introduce the PST for the $\delta$-SPH model (see section 11.2). Then, in the context where the ALE framework is used as a quasi-Lagrangian framework (i.e. where the PST introduces a $\hat{u}$ tending to zero at convergence), this direct use of ALE equations, where both particle mass and volume evolve, can be compared to using the same equations where the following quasi-Lagrangian assumption is made with: (i) either enforcing constant particle masses as is usually done in Lagrangian schemes, or (ii) enforcing constant particle volumes. This comparison among those schemes is performed in [16] where it is also shown that directly using the ALE equations imposes adding a mass diffusion operator to stabilize the scheme, see section 11.2. In [205] it is noted that the use of a Riemann solver for evaluating the fluxes of the conservative variable $F(\Phi)$ (see section 8) automatically introduces enough numerical diffusion to obtain a stable scheme in this ALE framework, see section 11.2.

## 10.5. Volume conservation issue

As underlined in section 5, geometrical volume conservation is not ensured within the SPH models, namely the actual volume occupied by a particle in physical space is not linked to its transported volume $V_i$. For the SPH schemes $S_1$ and $S_2$ the particle volumes are directly linked to the particle positions, hence large volume errors are not expected if the ratio $h/\Delta x$ remains large enough during the time evolution. The latter condition can be compromised when dealing with high compressibility effects in the gaseous phases. As already remarked, in that flow condition a variable-$h$ scheme would be needed. In addition to that solution, in section 13.2 an effective correction is discussed which allows for maintaining a limited variation of the smoothing length.

Regarding the schemes $S_{3,4,5,6}$ the main source of error on the volume is simply due to the time integration of the continuity equations. However, with the small time steps needed for

a stable time integration, see section 11, error accumulation is expected to be very limited for these four SPH schemes as shown in, e.g. [213].

A different situation holds when corrections are added to SPH schemes, such as the PST presented in the above subsections, or the numerical diffusion terms discussed in section 11. As emphasized in [131, 235] these corrections, beside their benefits, introduce secondary terms where volume errors accumulate with the same sign. As a consequence volume errors can become non negligible when simulating violent free-surface flows for long-time evolution. In this kind of problem the fragmentation of the interfaces is the main source of volume errors. It is worth noting that these errors are not strictly related to the weakly-compressible assumption since in [105, 198, 282] the authors show that even with incompressible SPH models the geometrical volume conservation is an important issue, and this can even be generally extended to other particle methods [233].

Possible solutions to this issue have recently been proposed in the literature. For example, it is shown in [96, 125, 160] that PSTs can be designed in order to control volume error. Conversely, in [213, 273] the ALE SPH model exhibits some advantages for volume conservation due to the possibility of controlling both particle mass and volume evolutions.

## 10.6. Detection of free-surface particles

As discussed in the previous subsections the PST implies the detection of the particles belonging to the free surface. One of the most widely used algorithms is the one proposed in [171]. This method uses a two-stage procedure.

In the first stage only the intrinsic properties of the SPH interpolation are used, making it easy to embed in any SPH solver. As proposed in [75], a subset of the particles which are close to the free surface can be efficiently identified through the values of the minimum eigenvalue $\lambda_i$ of the renormalization tensor $L_i^{-1}$, see equation (47).

In a second stage a more precise and reliable control is performed on this subset in order to complete the free-surface detection. This second step requires the evaluation of the normal vector field to the free surface which is performed using again the kernel properties:
$$ \mathbf{n}_i = \frac{\nabla c_i}{|\nabla c_i|}, \quad c_i := \mathbb{L}_i \sum_j \nabla_i W_{ij} V_j. \quad (133) $$

Here, compared to the normal vectors evaluation of equation (93), the renormalization tensor $\mathbb{L}$ is used. The latter helps to get a normal vector field that is less dependent on the particle disorder. In [236] the properties of the renormalization tensor are further inspected, leading to the following improvement of the normal vector field evaluation:

$$ \mathbf{n}_i = \frac{\langle \nabla \lambda_i \rangle}{|\langle \nabla \lambda_i \rangle|}, \quad \langle \nabla \lambda_i \rangle := \sum_j (\lambda_j - \lambda_i) \nabla_i W_{ij} V_j. \quad (134) $$

This form was then used in [258] for the evaluation of the surface tension for single-phase flow, see section 7.8.

This free-surface particle detection algorithm of [171] applies to both 2D and 3D frameworks and does not need any complex operation linked to the geometrical distribution and intersections of the particles modeled as spheres, conversely to in [72, 107]. Furthermore, this algorithm can also be used for post-processing purposes and for assigning a LS function on underlying Cartesian meshes.

There are later articles where different detection algorithms are presented, for example in [149, 222]. The latter is based on a ‘Hidden Point Removal’ operator which determines the visible points in a point cloud as viewed from a given viewpoint, while the former uses a Delaunay triangulation. In [240] following the ‘arc method’ of [72] another detection algorithm called ‘central angle method’ is proposed which is simpler and computationally more advantageous than the ‘arc method’.

# 11. Stability and convergence of the SPH model

## 11.1. Time integration of the SPH schemes

The SPH systems of equations $S_{1..6}$ are usually integrated in time using an explicit scheme. In the literature different time integration schemes were discussed, one of the most widely used being the Modified Euler Scheme. Symplectic integration algorithms have also been adopted, see, e.g. [60, 191]. The latter have been developed within the context of molecular dynamics [144] and have the advantage of being time reversible in the absence of friction or viscous effects. However, for the SPH systems of equations there is no evident convenience of this class of integration schemes with respect to other classic schemes, see, e.g. [13].

Conversely, high-order Runge-Kutta schemes appear to be more computationally efficient, allowing larger time steps with respect to other integration schemes [54, 203]. As discussed in [12, 235] another advantage of high-order Runge-Kutta schemes is the possibility, classically used with other methods, to partly freeze the update of the time-evolved quantities, e.g. by simple linear extrapolation, which limits computational costs, improving also the stability and the robustness of the SPH scheme.

Besides the choice of integration scheme its explicit nature imposes a stability criterion on the time step, $\Delta t$, obtained as

the minimum of the limitations imposed by physical information propagation times:

$$ \left\{ \begin{array}{ll} \Delta t_c = K_c \left( \frac{h}{c_0} \right), & \Delta t_v = K_v \left( \frac{h^2}{\nu} \right), \\ \Delta t_\sigma = K_\sigma \sqrt{\frac{(\rho_X + \rho_Y)/2}{2\pi\sigma_{XY}}} h^3, & \Delta t_a = K_a \min_i \sqrt{\frac{h}{\|\mathbf{a}_i\|}}, \\ \Delta t = \min(\Delta t_c, \Delta t_v, \Delta t_\sigma, \Delta t_a), & \end{array} \right. \quad (135) $$

where $\|\mathbf{a}_i\|$ is the particle acceleration, and the Courant-Friedrichs-Lewy constants $K_c$, $K_v$, $K_\sigma$ and $K_a$ depend on the specific time integration scheme adopted. The first constraint $\Delta t_c$ is linked to the weak-compressibility of the medium and to the related speed of sound $c_0$, and it usually represents the most restrictive condition, as the constraint of equation (13) needs to be satisfied for the liquid phases. The second and the third constraints $\Delta t_v$ and $\Delta t_\sigma$ are linked to viscous diffusion and to the surface tension modeling [103, 193]. When working with liquids such as water or oil the low magnitude of the kinematic viscosity $\nu$ and surface tension $\sigma_{XY}$ coefficients make these constraints the least restrictive in most of the free-surface or multiphase problems addressed in the literature. Note also that when using equation (13) to define $c_0$, the surface tension constraint is already embedded within $\Delta t_c$. A theoretical analysis of the first two constraints of equations (135) was provided in [263] where it was also shown that the choice of the kernel shape used in the SPH model has not a great influence on the parameter $K_c$, when the kernel standard deviation is considered instead of the usual smoothing length $h$. As stressed in [57], a limit on the time step for the advection of the particles should be considered since the SPH method is sensitive to the particle spatial distribution, therefore the last constraint $\Delta t_a$ needs to be added. Note that for incompressible SPH schemes the acoustic time constraint does not exist, see section 12.

## 11.2. Stability of the scheme adding numerical dissipation

As commented in the previous section the SPH schemes $S_{1..6}$ are integrated in time explicitly and the spatial discretization can be seen as centered since the support of the kernel is spherical and is centered at the particle $i$ location. Therefore, seen as a generalization of, e.g. a Finite Difference scheme, these SPH schemes would be classified as unstable and can be stabilized through the introduction of numerical dissipation, without which the pressure field becomes extremely noisy.

Two traditional ways to achieve stabilization then arise from the literature of the numerical solution of the compressible Euler equations, dating from the 1950s: adding an artificial viscosity [265], or using a Riemann solver to evolve the solution in time [99].

In the SPH framework the simplest way, namely adding an artificial viscosity, was first followed where an additional term appears in the momentum equation. Different forms of this term exist, the most frequently used is the one presented
in [189] which was reshaped in [141] as:

$$ m_i f_i^{\text{av}} = -\alpha c_0 \rho_0 h \sum_j \min[\pi_{ij}, 0] F_{ij} r_{ij} V_i V_j, \\ \pi_{ij} := \left( \frac{\mathbf{u}_{ij} \cdot \mathbf{r}_{ij}}{r_{ij}^2} \right). \qquad (136) $$

The intensity of the $\alpha$ parameter tunes the level of the artificial viscosity. The latter is proportional to the smoothing length, so that, when the spatial resolution increases (that is, when $h$ decreases), the action of the artificial viscosity becomes smaller and smaller, and the consistency with the Euler equations is recovered. Generally, the switch through the sign of the $\pi_{ij}$ operator is useful in gas dynamics to improve the SPH stability in the neighborhood of shock waves. In the context of the present review this switch is usually avoided in order to get an artificial viscosity which is consistent with the real one of equation (88), see, e.g. [11].

Following the path of using Riemann solvers instead, Vila and Ben Moussa proposed [259] and studied [24] SPH schemes based on the Godunov method and the use of Riemann solvers for approximating the fluxes of conservative variables $F(\Phi)$, see section 8, discarding the use of an explicit artificial viscosity term such as the one of equation (136). They used weak forms of the Euler equations and gave proof of convergence with the use of these Riemann solvers in [24], extracting in particular the convergence conditions. Such a Riemann-based SPH scheme, also called Godunov-SPH scheme, relies on solving a Riemann problem between each pair of interacting particles, with a left and a right state reconstructed at a fictitious contact discontinuity located at the midpoint of the two particles. As mentioned by Vila [259] different Riemann solvers can be used, either exact or approximate, see, e.g. [214]. It must be noted that an important consequence of this way of introducing Riemann solvers in the SPH scheme is that the mass of the particles varies in time, contrary to in the classic SPH schemes. This was studied especially in [205].

Alternatively, at the same period of time Parshikov et al [209] proposed a slightly different way to introduce Riemann solvers in which they directly replaced the particle interaction, e.g. $P_i + P_j$, by the solution of the approximate Riemann solver taken again at the midpoint of the two particles. As a consequence, the mass of the particles does not vary in time in that case, as is the case in the classic SPH schemes. Note that a previous attempt to introduce Riemann solvers was proposed by Inutsuka [123] where only one global 1D-Riemann problem was solved per particle $i$, but this was not conservative and was subsequently abandoned by the author who later proposed a Riemann-SPH scheme in [124] quite similar to that of [209].

Another advantage of using Riemann solvers is the capability of increasing the global order of convergence of the SPH method through the reconstruction order of the variables used in the Riemann problem, namely through how the particle variables are extrapolated to reconstruct the right and left states of the Riemann problem at the midpoint of the two considered

particles. As proposed in [123, 259], a first-order extrapolation can be used, through a MUSCL (Monotonic Upstream-centered Scheme for Conservation Laws) scheme [254], so as to gain one order of convergence, see [259]. This principle can then be extended to get high-order schemes, see, e.g. [18, 175, 257]. However, even though these high-order schemes prove to be effective and efficient in one-dimensional problems, their efficient application to three-dimensional problems remains an open question, see [257]. Another complexity when increasing the order lies in the fact that the extrapolation used to reconstruct the right and left states of the Riemann problem needs to be limited within physical bounds. The choice and use of such limiters can prove to be critical, see, e.g. [101, 195].

Through the use of Riemann solvers numerical diffusion is added in both continuity and momentum equations. This has the great advantage with respect to the artificial viscosity term (136) of better stabilizing the pressure field by filtering out spurious numerical high-frequency noise. In the last two decades the successful use of Riemann-based SPH solvers has largely increased in the community, providing accurate results on a variety of problems, see, e.g. [178].

Exploiting the theory of Riemann solvers, a simplified numerical diffusive term based on a Rusanov flux was then built in [89]. This term was added only into the continuity equation. A similar approach was also proposed by [181]. However, the latter two techniques present the problem of not conserving the total volume of the fluid system since the hydrostatic component is filtered improperly, making them inapplicable for long-time simulations. To overcome this drawback Antuono et al [12] proposed the $\delta$-SPH scheme based on the following correction:

$$ \left\{ \begin{aligned} \mathcal{D}_i^\rho &:= \delta c_0 h \sum_j \psi_{ji} F_{ij} V_j, \\ \psi_{ji} &:= 2 \left[ (\rho_j - \rho_i) - \frac{1}{2} \left( \langle \nabla \rho \rangle_i^L + \langle \nabla \rho \rangle_j^L \right) \cdot \mathbf{r}_{ji} \right], \end{aligned} \right. \qquad (137) $$

where the term $\mathcal{D}_i^\rho$ is the numerical density diffusion term added in the right-hand side of the SPH continuity equation, e.g. the first equation of equations (116). The superscript $L$ in equations (137) indicates that the gradient is evaluated through the renormalized gradient, see equation (47). The dimensionless parameter $\delta$ is related to the magnitude of the numerical diffusive contribution added for stability purposes. In [13] it is shown that the numerical density diffusion $\mathcal{D}^\rho$ acts as a pure energy dissipation in the particle system. During the last decade the $\delta$-SPH model has been successfully applied to several hydrodynamics problems, proving to be accurate and robust, and it has become one of the most popular SPH variants, see, e.g. [101].

An example of the benefits of using the density diffusion on removing pressure instabilities is given in figure 12, from [12]. The plots show the pressure fields inside an inclined jet impinging on a flat plate without (left) and with (right) adding density diffusion (137).

In [106] the extension of this model to the multiphase context is provided and compared against a Riemann-SPH model.
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Pressure field inside an inclined jet flow impinging on a flat plate. Result without (left) and with (right) density diffusion. Note that in this figure the density diffusion coefficient is noted $\xi$ and not $\delta$. Reprinted from [12], Copyright (2010), with permission from Elsevier.</figcaption></figure>

This article also provides the formulation of the $\delta$-SPH scheme for the scheme $S_6$ (132) where the continuity equation is written for particle volumes instead of particle densities.

In [177] the benefits of using a combination of equations (136) and (137) for stabilizing SPH schemes are discussed. The two corrections act in a complex and non-linear way; they are not independent and it is possible to optimize their actions by changing the $\alpha$ and $\delta$ parameters for each particle using simple sub-scale models. The new scheme derived stabilizes the SPH scheme while drastically reducing the level of numerical dissipation, making it possible to simulate high Reynolds number flow using Large Eddy Simulation models. That new SPH model was named $\delta$-LES-SPH, see [14, 177].

### 11.3. Tensile instability

In the two studies [71, 72] different SPH formulations were investigated, considering the connections between consistency, stability, conservation properties and enforcement of the boundary conditions. In particular it is stressed that conservative formulations (in terms of momenta and energy), such as schemes $S_{1..6}$ discussed in this review, lead to schemes that are robust schemes but that have a low accuracy. Furthermore such schemes are unstable in specific pressure/tension conditions as demonstrated by [242], where this specific instability is called ‘Tensile Instability’. Conversely non-conservative formulations based for example on the use of MLS interpolation (see section 3) can present higher accuracy, free from tensile instability, but the drawback is their lack of robustness especially when the fluid/solid domain is subjected to large deformation. Recently, this robustness issue was solved in [179].

In the context of solid deformation different techniques were developed for removing tensile instability. For example in [28] it is demonstrated that this instability can be removed if the particle interactions are performed while keeping the

same neighboring particles, that is using the so-called ‘Total Lagrangian Formalism’. Another possible solution is the use of so-called ‘stress points’, proposed in [77], in which the stresses are calculated at points other than the SPH centroids in order to remove the instability.

These kinds of approaches are usable when considering deformation of solids but for fluid dynamics problems the deformations are generally very large, making them inapplicable. In the fluid dynamics context it was shown that the use of PST, aside the accuracy improvement of the solution, can also alleviate tensile instability but not completely remove it, see [205, 236]. For example in [14, 234] it is shown that, especially for high Reynolds number regimes, tensile instability can induce numerical cavitation in the cores of eddies, even using the PST, due to the pressure drop inside these flow regions.

One way to eliminate the instability is to introduce a switch in the pressure force:

$$ \left\{ \begin{array}{l} f_i^p = - \sum_i (p_j + p_i) \nabla_i W_{ij} V_j + S_i \sum_j \nabla_i W_{ij} V_j, \\ S_i = \begin{cases} 0 & \text{for } p_i \ge 0 \text{ or } i \in S_F, \\ 2p_i & \text{for } p_i < 0, i \notin S_F, \end{cases} \end{array} \right. \qquad (138) $$

where $S_F$ denotes the region of the fluid domain close to the free surface, that is the free-surface particles and their neighbor particles. The free surface particles are detected through the algorithm described in [171], see section 10.6. Equation (138) corresponds to a switch from the ‘plus’ formulation (namely, $p_j + p_i$) to the ‘minus’ formulation (that is, $p_j - p_i$) when the pressure $p_i$ is negative inside the fluid region. It is worth stressing that equation (138) can be used only in combination with the PST; indeed, if this switch is introduced alone in an SPH scheme it drastically reduces its robustness as remarked in [14, 141].

The left plots of figure 13 depict the development of numerical cavitation due to tensile instability for the flow past a
<figure><img src="image_17.png" /><figcaption>**Figure 13.** Flow past a NACA0010 airfoil with an angle of attack of 30° at Reynolds number 10 000. Snapshots of the SPH solutions without (left) and with (right) Tensile Instability Control (equation (138)). Reprinted from [234], Copyright (2018), with permission from Elsevier.</figcaption></figure>

NACA0010 profile with an angle of attack of 30° at Reynolds number 10 000. The right plots of the same figure show the flow evolution using the correction (138).

### 11.4. Pairing instability

Beside tensile instability, in [67] a second instability related to the SPH scheme is discussed, called ‘pairing instability’. This instability is caused by the negativity of the multi-dimensional Fourier transform of the kernel which induces, for low $\Delta r/h$ ratios, i.e. for large numbers of neighbors, the clumping in pairs of the particles. This clumping can strongly reduce the convergence speed of the SPH scheme. The analysis performed in [67] shows that, so far, the Wendland functions [267] are the only kernels which are stable against pairing for all the $\Delta r/h$ ratios.

### 11.5. Convergence of the SPH schemes

In previous sections the consistency and the stability of the SPH schemes were addressed. The consistency was discussed in terms of the SPH integral interpolation and in terms of the quadrature formulas adopted to discretize the convolution integral with the kernel function, see section 3. In section 7,

the consistency of the SPH equations to the NSEs was demonstrated. However, the convergence of the SPH schemes $S_{1..6}$ is still missing.

Different attempts to address the convergence of SPH schemes are reported in the literature, at least for the Euler equations (inviscid flow) written in a weak form. The most relevant studies are [23, 24, 68, 173, 259]. The main result of these studies is that the only way to ensure convergence is to decrease $h$ much more slowly than $\Delta r$ when refining, due to the double approximation made in the SPH interpolation, first mollification and then quadrature, see section 3. This consequently implies that the number of neighbor particles grows rapidly when refining. Conversely, for obvious calculation time reasons the $\Delta r/h$ ratio is traditionally kept constant in SPH simulations. The impact of this choice on the discrete SPH integral interpolation error is discussed in section 3.3. A heuristic study of this choice of a constant number of neighbor particles is discussed in [283] with different simulated test cases, and in [141] it is underlined that a convergence rate between 1 and 2 can be attained in a range of spatial resolutions of practical use, if a small enough constant $\Delta r/h$ ratio is adopted.

Furthermore, the convergence is linked to the regularity of the particle distribution and, as highlighted in [10], adding
numerical dissipation to SPH schemes, beside improving their stability, helps in the reduction of particle disorder and consequently increases the interval where the scheme is convergent. The use of the PST, see section 10, which has become current practice in recent years, is thus crucial in this respect as well.

Nonetheless, as pointed out in [24, 71, 259], a way to guarantee the convergence of the SPH scheme with a fixed $\Delta r/h$ ratio is through the use of renormalized kernel gradients such as equations (46) and (47). However, the use of these higher-order schemes leads to better accuracy, thus a better respect of the Lagrangian trajectories and, in turn, less regular particle distributions, see section 10. As already discussed, this has conversely a negative impact on the accuracy and besides, the renormalization tensor becomes ill-conditioned, leading to stability issues, see [141]. Again, the use of the PST is thus crucial in this context of higher-order variants.

Finally, recent theoretical studies on the improvement of SPH convergence results are also given in [84, 91, 92].

## 12. Modeling liquid phase: weakly-compressible vs. incompressible approach

As stated in section 2, liquid phases are treated as barotropic fluids with an EOS that directly links the density to the pressure field. Within the constraint of equation (13), since the Mach number is sufficiently low to guarantee negligible compressible effects, the schemes $S_{1..6}$ are generally named ‘weakly-compressible SPH’ (wcSPH) schemes. Note that since equation (13) is an initial estimation, it is necessary to monitor the density variations during the simulation in order to check that the weakly-compressible constraint $\Delta\rho/\rho_{0x} < 0.01$ is not exceeded.

However, the liquid phase can also be treated by making the incompressible assumption, which is the usual choice made in mesh-based methods. The enforcement of zero velocity divergence can be made in the SPH method using the classical Projection Method by Chorin [46]. A first attempt at using a Projection Method for particle methods for free-surface incompressible viscous flow was made by Koshizuka et al [136], resulting in a method named moving particle semi-implicit method (MPS). It uses a semi-implicit time-marching scheme which is based on a fractional step algorithm subdivided as follows:

(i) The particles are advected only by viscous and external forces without considering the pressure force.

(ii) A PPE is solved in order to obtain a pressure force that corrects the velocity at the new time step by making it divergence-free.

(iii) The particles are advected at the new time step by averaging the velocity at previous time step with the new one.

As shown in [229] the MPS method can be seen as an SPH method with specific kernel choices, therefore [136] can be seen as the first work in which an incompressible SPH (iSPH) scheme is presented. During the last decades different researchers developed enhanced iSPH versions, see, e.g.

[63, 129, 157, 228, 269], while an extension to multiphase flow is provided in [121]. It is worth noting that, as remarked in section 10.3, the PST is a crucial tool to obtain robust and accurate solutions with iSPH models.

### 12.1. Explicit vs. semi-implicit solution

From a computational point of view, the main difference between wcSPH and iSPH models is that the acoustic time-step constraint in equations (135) is no longer needed in iSPH models and larger time step can then be adopted. On the other hand:

(i) At each time iteration iSPH needs the solution of a PPE, i.e. the solution of an algebraic system with a sparse matrix.

(ii) Explicit detection of the free surface is required in order to assign the dynamic boundary condition. Consequently, the detection algorithm (see section 10.6) may affect the solution of the Poisson equation when considering violent free-surface flows with very complex interface configurations.

(iii) The other time-step constraints of equations (135) still limit the time step which can thus only be increased by a maximum factor of 10 with respect to the wcSPH one.

The above three points make the CPU costs of iSPH models of the same order of magnitude as wcSPH ones. Furthermore, iSPH algorithms are more difficult to efficiently parallelize. This is the reason why 3D simulations with billions of particles have already been performed with wcSPH (see, e.g. [33, 60, 62, 204]), while, to the authors’ knowledge, to date there has been no application of iSPH involving such a large amount of particles, or solving problems with such geometrically highly-complex interfaces.

A possible compromise between iSPH and wcSPH is explicit incompressible schemes (eiSPH), see, e.g. [20, 64, 65, 119, 196, 216, 252]. These schemes adopt an approximated solution of the PPE. In this way the computational costs are reduced with respect to iSPH, while at the same time maintaining an explicit time integration. eiSPH schemes are able to greatly reduce the velocity divergence within the domain, though they cannot enforce strict zero divergence, which has consequences on the regularity of the predicted pressure field.

### 12.2. Differences in the physical modeling

In the previous section it is highlighted that the liquid phases can be treated as incompressible or as weakly-compressible media. The two different approaches imply differences in the solutions obtained.

Following [194, 224], in the limit of Mach number going to zero, the pressure $p$ given by the solution of the compressible NSEs can be decomposed as $p = p_I + p_c$, where $p_I$ is the pressure solution of the incompressible NSEs and $p_c$ is the acoustic pressure obtained by the d’Alembert wave equation, related to the non-null divergence component of the velocity field. It
should be underlined that the above decomposition is phys-
ical and related to the governing equations and not to the spe-
cific numerical schemes at hand. As remarked in [59], in free-
surface problems the displacement of the air–water interface is
practically unaffected by the component $p_c$. This is something
that is well known within the SPH community, where the free-
surface time evolutions obtained by iSPH or wcSPH are very
similar while the pressure fields can present large differences.

Since the speeds of sound used for the liquid phases are not
the physical ones in wcSPH simulations, the acoustic pressure
$p_c$ is an unwanted output of the simulations, and it can be gen-
erally removed from time series outputs of SPH simulations
through filtering procedures, as discussed, e.g. in [176].

There are different techniques to partially remove $p_c$ dur-
ing the time evolution through approximated incompressible
flow solvers. One of these is through the use of eiSPH as
remarked in section 12.1. A second technique was recently
explored in [131, 239] using a numerical secondary viscosity
in order to quickly damp acoustic waves. This technique is also
used in some Finite Volume Methods where an incompress-
ible flow constraint is enforced using a pseudo-compressibility
approach. Within this kind of schemes the numerical sec-
ondary viscosity is used to speed up the convergence
towards the incompressible NSEs in the pseudo-time, see,
e.g. [35, 97, 147].

12.3. Liquid impacts

From a physical point of view, the air phase plays a crucial
role during liquid impacts, especially when the liquid–liquid
or liquid–solid collision at stake occurs with a flat or close-to-
flat angle, that is the air–liquid interface is parallel to the other
air–liquid or air–solid interface. There, the cushioning effect
of the air trapped between the two interfaces prevents a singu-
lar behavior of the impact. In this context, it is also important
to recall the analysis performed by Peregrine [211] regarding
the effect of the presence of dispersed air bubbles in the bulk of
water during an impact event: even a small volume fraction of
air in water greatly increases its compressibility. This induces
in particular a drastic variation in the speed of sound. At atmo-
spheric pressure, just 1% of dispersed air decreases the speed
of sound of the mixture to 120 m s⁻¹, i.e. more than a factor 10
compared to pure water, and it drops to 30 m s⁻¹ for 20% of
dispersed air. It is worth noting that the speeds of sound typic-
ally adopted in wcSPH models are often in the range of these
values, even if the speed of sound of an air–water mixture is,
in reality, not constant in space and time, depending on the air
bubble concentration and migration, see also [163, 176, 212].

However, in SPH simulations liquid impacts are often
modeled by single-phase free-surface models, i.e. neglecting
the presence of air, transformed into void. This has several
implications since the nature of the impact becomes singular,
especially when the impact angle is flat. Under this assump-
tion of neglecting the presence of air, due to its compressib-
ility part of the mechanical energy of the liquid is converted
into elastic energy during the impact. That elastic energy is
generally dissipated by viscous effects after several oscillation
cycles during which the generated traveling acoustic waves in

the fluid domain transform the elastic energy into mechanical
one and vice-versa. If an incompressible assumption is made,
on the contrary, instantaneous mechanical energy losses occur
as pointed out in [245].

Thus, when making the weakly-compressible assumption
as is usual in single-phase free-surface SPH models, specific
attention has to be paid to these impact situations. First, the
acoustic waves originating during liquid impacts can be char-
acterized by large intensities and, in some cases, cannot be
simply considered as a perturbation of the incompressible
component since the compressible pressure component, $p_c$,
can be of the same order of magnitude as the incompress-
able one, $p_I$, or even much larger. It follows that the filter-
ing procedure needs to be performed carefully since remov-
ing the acoustic component is not always possible, depend-
ing especially on the continuous or discontinuous behavior
of the mechanical energy of the incompressible NSEs solu-
tion at the liquid impact, see [169]. For example, for a flat
liquid impact the maximum pressure predicted by an incom-
pressible model is proportional to $I/\Delta t$, I being the pres-
sure impulse, whereas it is proportional to 1/Ma under the
weakly-compressible assumption. Therefore, in the limit of
the incompressible solution Ma → 0, or equivalently in the
limit of $\Delta t \to 0$ for an incompressible model, the impact pres-
sure becomes singular, clearly indicating the inability of these
single-phase free-surface models to deal with this kind of phe-
nomenon, see [134, 169]. Finally, as underlined in [169], dur-
ing these impacts the pressure scales with $\rho c_0 U$, U being the
flow velocity scale, rather than with $\rho U^2$ in a non-impact situ-
ation, making the condition (13) more severe, drastically reduc-
ing the time step $\Delta t_c$ in equation (135).

13. SPH simulation of liquid–gas flows

Thanks to the Lagrangian nature of the SPH method, phases
always remain strictly separated, i.e. (i) no mass flux occurs
between phases, and (ii) interfaces are always sharp and
simply defined as separating two neighbor particles of differ-
ent phases. These features are advantageous with respect to
other numerical techniques such as the VoF or the LS tech-
niques used in conventional mesh-based methods.

However, even with these sharp interfaces defined between
phases a direct enforcement of complex interface conditions
(see section 2.1) is generally not performed in SPH models:

(i) Density discontinuities can be handled through the SPH
operators as discussed in section 5.

(ii) The SPH pressure force formulas of the schemes $S_{2,4,5,6}$ are based on a direct interaction between particles belonging to different phases. Therefore, similarly to what is discussed in [49, 50] for the free surface, an explicit enforcement of the dynamic boundary condition (8) is not performed either for interfaces between two phases.

(iii) The SPH particle viscous and thermal stresses are simply given through harmonic means of the particle viscosities and thermal conductivities, see section 6.3.1.