# A study on the energy consistency in SPH surface tension modelling

Salvatore Marrone, Matteo Antuono, Antonio Agresta, Andrea Colagrossi

► To cite this version:

Salvatore Marrone, Matteo Antuono, Antonio Agresta, Andrea Colagrossi. A study on the energy con- sistency in SPH surface tension modelling. Computer Methods in Applied Mechanics and Engineering, 2024, 433, pp.117473. 10.1016/j.cma.2024.117473. hal-05173678

HAL Id: hal-05173678

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

Submitted on 21 Jul 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>
<figure><img src="image_2.png" /></figure>

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

# A study on the energy consistency in SPH surface tension modelling

S. Marrone<sup>a</sup>, M. Antuono<sup>a</sup>, A. Agresta<sup>b,*</sup>, A. Colagrossi<sup>a,c</sup>

## ARTICLE INFO

**Keywords:**

Smoothed Particle Hydrodynamics

Free-surface flows

Surface tension

Energy balance

Droplet oscillations

Sloshing flows

## ABSTRACT

In the present work the evolution of viscous drops oscillating under the action of surface tension is tackled. Thanks to its structure, the SPH scheme allows for an analysis of the energy balance that is rarely addressed to in the general Computational Fluid Dynamics literature for this kind of flows. A procedure for checking the consistency between the energy of the surface-tension force and the free-surface evolution is proposed. Such a procedure relies on well-known analytical relations for the surface tension and on the evaluation of the free-surface area through a level-set function. Several test cases, in both two and three dimensional frameworks, are considered for validation. The study is performed by selecting a specific SPH scheme with a specific single-phase surface tension model. In any case, the procedure proposed is general and extendable to other SPH surface tension models and SPH schemes.

## 1. Introduction

The definition of an appropriate surface tension models in SPH has been an active and prolific research subject in the last decade (see e.g. [1–7]). Generally, the most widespread approach for the inclusion of surface tension in SPH relies on the so-called Continuum Surface Force (CSF) method, which dates back to the pioneering work of [8]. Succinctly, the main idea is to represent the force along the interface in the form of a volumetric contribution inside the fluid. This can be done by modelling the surface tension through the divergence of a proper tensor (see e.g. [9]) or through the definition of a volume force distributed along the interface and vanishing inside the fluid bulk (see e.g. the pioneering work by [10] within the SPH framework). The latter approach (which is the one selected in our study) is somehow more complex since it involves the explicit computation of the surface curvature, but it has proven to be rather robust and accurate. In particular, in the present work we consider the model of Vergnaud et al. [11]. This choice is motivated by the fact that such a model has been validated on many benchmarks, showing to accurately model the physics of the surface tension in comparison to other simpler surface-tension models (see, for example, [10,12–14]).

The idea at the basis of the present work is to analyse how the energy balance of the SPH evolves when a surface tension model is implemented. Due to the complexity of this problem, the analysis is limited to single-phase problems, while the study of multiphase cases is postponed to future works. In the last decade, the study of the energy balance in SPH has attracted an increasing interest because of two main factors: (i) the simplicity in identifying the different contributions of the energy budget, (ii) the possibility of inspecting in detail the properties of energy conservation of the scheme. This latter aspect is becoming more and more compelling as the SPH models evolve towards non-conservative formulations (namely, those formulations that do not allow for an exact conservation of momenta and energy). In this latter case the study of the energy balance permits to check the extent of the deviations from the theoretical conserved quantity and to inspect the presence of spurious numerical dissipation. We

<footnote><sup>a</sup> CNR-INM, Institute of Marine Engineering, Via di Vallerano 139, 00128 Roma, Italy</footnote>

<footnote><sup>b</sup> ENEA - Centro Ricerche Casaccia, Via Anguillarese 301, 00123, Rome, Italy</footnote>

<footnote><sup>c</sup> École Centrale Nantes, LHEEA Lab. (ECN/CNRS), Nantes, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: antonio.agresta@enea.it (A. Agresta).</footnote>
highlight that the study of the energy balance described in the present work belongs to a line of research developed by the authors in over a decade. This starts with the analysis of the viscous effects in SPH [15–17] and goes on with the investigation of the energy balance for its diffusive variants [18,19]. The same kind of analysis is also applied in several different works, including wave impact problems and violent sloshing dynamics in 2D and 3D (see [20–22]).

In the present work the surface tension model of Vergnaud et al. [11] is implemented inside the Regularized High-Order Diffusive Smoothed Particle Hydrodynamics scheme (RHOD-SPH) described in Michel et al. [23]. This scheme belongs to the class of non-conservative SPH solvers and proved to be reliable and very accurate thanks to an effective action of the shifting technique, to the use of approximate Riemann solvers and to the control of the tensile instability (see [23] for more details). Incidentally, we observe that the model of Vergnaud et al. [11] is non-conservative as well, since it does not guarantee the perfect conservation of the linear momentum. In the presence of surface tensions, the analysis of the energy balance benefits from some theoretical results which allow us to check the accuracy of both the RHOD-SPH and the Vergnaud’s model. In particular, it is possible to represent the work done by the surface tension in terms of a geometrical measure of the free-surface length. The latter one is computed by interpolating the particles’ distance on a discrete grid and by using a level-set function. Such an approach allows us to extract the drop surface and to compute its area, checking the consistency between the geometrical evolution of the free surface and the work done along it. Furthermore, the analytical findings described in [24] are used to verify the correctness of the viscous damping during the drop oscillations and the corresponding energy contributions during the evolution.

The paper is structured as follows: in Sections 2 and 3 we briefly introduce the governing equation, the main theoretical results and the RHOD-SPH, then in Section 4 we describe the numerical applications. These contain 2D and 3D problems with singularities along the free surface (i.e. squared drops) and the evolution of 2D circular drops subjected to large deformations.

## 2. Governing equation

In the present work we consider the evolution of a single-phase viscous fluid in the domain $\Omega$ subjected to surface tension along the frontier $\partial\Omega$ and without external forces. Phase changes and miscibility phenomena are not considered. The fluid is assumed to be weakly-compressible (more details are given below) and barotropic, so that it is governed by the weakly-compressible Navier–Stokes Equations (NSEs):

$$ 
\left\{ 
\begin{aligned}
\frac{D\rho}{Dt} &= -\rho \nabla \cdot \mathbf{u}, \\
\frac{D\mathbf{u}}{Dt} &= -\frac{\nabla\rho}{\rho} + \frac{\nabla \cdot \mathbf{V}}{\rho} + f^{\sigma}, \\
\frac{D\mathbf{r}}{Dt} &= \mathbf{u}, \quad p = c_{0}^{2}(\rho - \rho_{0}).
\end{aligned}
\right. 
\qquad (1) $$

where $D/Dt$ indicates the Lagrangian derivative with respect to the fluid velocity $\mathbf{u}, \mathbf{r}$ is the position of the material points, $p$ is the fluid pressure, $\rho$ the fluid density and $\mathbf{V}$ denotes the viscous stress tensor. Finally, $\rho f^{\sigma}$ represents the surface tension which, in the present work, is modelled as a concentrated volumetric force acting on the interfaces $\partial\Omega$ as described below:

$$ \rho f^{\sigma}(\mathbf{r}) = \sigma \kappa(\mathbf{r}) n(\mathbf{r}) \delta_{S}, \quad \kappa(\mathbf{r}) = -\nabla \cdot \mathbf{n}(\mathbf{r}), \quad \forall \mathbf{r} \in \Omega \qquad (2) $$

where $\sigma$ is the surface tension coefficient, $\delta_S$ is a Dirac distribution along the interface $\partial\Omega$ and $\mathbf{n}$ is the normal unit vector of $\partial\Omega$ pointing outward $\Omega$. Accordingly, the dynamic boundary condition at the free surface $\partial\Omega$ including the surface tension action reads as:

$$ \mathbf{V} \cdot \mathbf{n} = (\rho + \sigma \kappa) \mathbf{n}, \quad \forall \mathbf{r} \in \Omega \qquad (3) $$

As indicated at the beginning of this section, the liquid phase is considered a weakly-compressible Newtonian fluid, that means that the density variations with respect to the reference value $\rho_0$ (i.e. the density at the free surface) have to be always smaller than 1%. Under this hypothesis, it is possible to use a linearized state equation (like the one adopted in the system (1)), where the speeds of sound $c_0$ can be modified for computational convenience. In order to ensure the attainment of the weakly-compressible regime, however, the choice of the speed of sound is subjected to the following constraints:

$$ c_0 \ge 10 \max \left( \Delta U_{max}, \sqrt{\frac{\Delta P_{max}}{\rho_0}}, \sqrt{\frac{\sigma |\kappa|_{max}}{\rho_0}} \right), \qquad (4) $$

where $|\kappa|_{max}$, $\Delta U_{max}$, $\Delta P_{max}$ are the maximum interface curvature (in absolute value), velocity variation and pressure variation expected during the flow evolution. The dimensionless numbers characterizing the motions described by the system (1) are the Reynolds and Weber numbers. These are respectively defined as:

$$ \text{Re} = \frac{U_{\text{ref}} L}{\nu}, \quad \text{We} = \frac{\rho_0 U_{\text{ref}}^2 L}{\sigma} $$

where $\nu$ is the kinematic viscosity and $U_{\text{ref}}$ and $L$ are the reference velocity and length of the problem under consideration.

Following the approach defined in [17] and adding the action of surface tension, the energy balance of system (1) reads:

$$ 
\left\{ 
\begin{aligned}
\dot{\mathcal{E}}_K + \dot{\mathcal{E}}_\sigma + \dot{\mathcal{E}}_C = \mathcal{P}_v, & \quad \mathcal{E}_K := \frac{1}{2} \int_{\Omega} (\mathbf{u} \cdot \mathbf{u}) \rho \, dV, & \quad \mathcal{E}_\sigma := \sigma \int_{\partial\Omega} dS, \\
\mathcal{P}_v := 2\mu \int_{\Omega} (\mathbf{D} : \mathbf{D}) \, dV, & \quad \mathcal{E}_C := c_0^2 \int_{\Omega} \left[ \log\left(\frac{\rho}{\rho_0}\right) + \frac{\rho_0}{\rho} - 1 \right] \rho \, dV,
\end{aligned}
\right. 
\qquad (5) $$
where $\mathcal{E}_K$ is the kinetic energy, $\mathcal{E}_\sigma$ the surface tension energy, $\mathcal{E}_C$ is the potential energy associated to the fluid compressibility and the dot indicates time derivation. The power dissipated by viscosity is denoted by $P_v$, where $\mathbb{D}$ is the rate-of-strain tensor. During the numerical simulations, we also check the energy dissipated by viscous effects, defined as:

$$Q_v = \int_{t_0}^{t} P_v d\tau. \qquad (6)$$

The terms $\mathcal{E}_\sigma$ and $\mathcal{E}_C$ represent pure reversible energies. The latter quantity, however, is always negligible because of the weak-compressibility assumption and, hence, it is not considered hereinafter. The initial energy components at $t = t_0$ will be denoted as $\mathcal{E}_{K0}$ and $\mathcal{E}_{\sigma0}$.

## 3. Numerical model

### 3.1. Brief description of the RHOD-SPH model

In this section the Regularized High-Order Diffusive Smoothed Particle Hydrodynamics scheme (RHOD-SPH) derived in Michel et al. [23] is briefly recalled. In particular the governing Eq. 2 are discretized within a Quasi-Lagrangian SPH model as below:

$$\left\{
\begin{aligned}
\frac{d\rho_i}{dt} &= -\rho_i \nabla \cdot (\mathbf{u}_i + \delta\mathbf{u}_i) + \nabla \cdot (\rho_i \delta\mathbf{u}_i) + \Theta_{i,Rie}^p \\
\frac{d\mathbf{u}_i}{dt} &= F_i^p + F_i^\mu + F_i^{ad} + F_i^\sigma + \nabla \cdot (\rho_i \mathbf{u}_i \otimes \delta\mathbf{u}_i) + \Theta_{i,Rie}^u \\
\frac{d\mathbf{r}_i}{dt} &= \mathbf{u}_i + \delta\mathbf{u}_i, \quad V_i(t) = m_i / \rho_i(t), \quad p_i = c_0^2 (\rho_i - \rho_0),
\end{aligned}
\right.
\qquad (7)$$

where the subscript i refers to the generic ith particle while $F_i^p$, $F_i^\mu$, $F_i^{ad}$ and $F_i^\sigma$ are respectively the pressure, the viscous, the acoustic damper defined in [25] and the surface tension forces acting on the particles.

The time derivative $d/dt$ in (7) indicates a quasi-Lagrangian derivative i.e.:

$$\frac{d(\bullet)}{dt} := \frac{\partial(\bullet)}{\partial t} + \nabla(\bullet) \cdot (\mathbf{u} + \delta\mathbf{u}),$$

since the particles move with the modified velocity $(\mathbf{u} + \delta\mathbf{u})$.

#### 3.1.1. Particle Shifting Technique (PST)

The vector $\delta\mathbf{u}$ is the so-called Particle Shifting velocity and it is used to regularize the particle distribution during the motion. The law chosen for $\delta\mathbf{u}$ is the same reported in [23]. Accordingly, the continuity and the momentum equations contain terms depending on $\delta\mathbf{u}$. These are discretized as follows:

$$\left\{
\begin{aligned}
\nabla \cdot (\mathbf{u}_i + \delta\mathbf{u}_i) &= \sum_j (\mathbf{u}_j - \mathbf{u}_i + \delta\mathbf{u}_j - \delta\mathbf{u}_i) \cdot \nabla_i W_{ij} V_j \\
\nabla \cdot (\rho_i \delta\mathbf{u}_i) &= \sum_j (\rho_i \delta\mathbf{u}_i + \rho_j \delta\mathbf{u}_j) \cdot \nabla_i W_{ij} V_j \\
\nabla \cdot (\rho_i \mathbf{u}_i \otimes \delta\mathbf{u}_i) &= \sum_j (\rho_i \mathbf{u}_i \otimes \delta\mathbf{u}_i + \rho_j \mathbf{u}_j \otimes \delta\mathbf{u}_j) \nabla_i W_{ij} V_j
\end{aligned}
\right.
\qquad (8)$$

where the subscript j indicates the neighbouring particles to the ith particle, while $W_{ij} = W(\|\mathbf{r}_j - \mathbf{r}_i\|)$ is the kernel function (hereinafter a $C^2$-Wendland kernel is adopted [26]).

The particle masses $m_i$ remain constant in time and they are set using the initial particle volumes, $V_{i0}$, and the initial particle densities $\rho_{i0}$, i.e. $m_i = \rho_{i0} V_{i0}$. The particles are initially set on a lattice with an almost homogeneous spacing $\Delta x$, and, consequently, the initial particles' volumes $V_{i0}$ are evaluated as $\Delta x^n$ where n is the number of spatial dimensions (see [27]).

Finally, the radius of the support of the kernel function is set to $2h = 4\Delta x$ in two dimensions and $2h = 3\Delta x$ in three dimensions, where h is the smoothing length.

It is worth noting that the ratio $h/\Delta x = 2$ implies 50 neighbour particles in the kernel support in two dimensions and 266 in three dimensions. In order to lower the latter number and, as a consequence, limit the CPU costs, it is a common practice in SPH to reduce the ratio in 3D simulations.

#### 3.1.2. Numerical diffusion and tensile instability prevention

The symbols $\Theta_{i,Rie}^p$ and $\Theta_{i,Rie}^u$ in Eq. (7) denote the numerical diffusive based on Riemann solvers and added to the scheme for stability purposes. These are expressed following the approach described in Parshikov et al. [28] and, for the sake of brevity, they are not reported in this work. The interested reader can find the details in [23].

Within the RHOD-SPH scheme the pressure force $F_i^p$ is modelled as:

$$F_i^p = \begin{cases} -\rho_i \mathbb{L}_i \sum_j (p_j - p_i) \nabla_i W_{ij} V_j & \text{if } i \in I \\ -\rho_i \sum_j (p_j + p_i) \nabla_i W_{ij} V_j & \text{otherwise} \end{cases} \qquad (9)$$
where the subset I refers to the inner fluid particles whose kernel support does not intersect the free surface (for more details see [23]) and L_i is the renormalization matrix (see [29]) which ensures the recovering of the exact gradient of a linear pressure field. The switch in Eq. (9) changes from a non-conservative formulation on the inner particles to a conservative formulation for the particle belonging to the thin region close to the free-surface. Even though the momenta conservation is lost, the proposed formulation is invariant by change of reference pressure (for internal flows), free from the Tensile Instability, and more accurate in the evaluation of the pressure forces with respect to the conservative scheme.

It is important to remark that during the droplet oscillations the surface tension can induce large negative pressure values (with respect to the zero level of the free-surface), potentially leading to the onset of the tensile instability. If the analysis is limited to high viscous regimes, the viscous forces can counteract this phenomenon. Anyway this is not the case when large Reynolds number regimes are considered, as in the present work (see e.g. Section 4.3).

### 3.1.3. Viscous forces and damping of the acoustic waves

The viscous forces are modelled following Monaghan & Gingold [30]:

$$F_i^\mu = 2\mu(n+2) \sum_j \pi_{ij} \nabla_i W_{ij} V_j, \quad \pi_{ij} = \frac{(u_j - u_i) \cdot (r_j - r_i)}{\|r_j - r_i\|^2}. \qquad (10)$$

This classical viscous model is selected following the work by [15] where it is demonstrated that formula (10) is consistent for free-surface flows. On the contrary, the formulation proposed in Morris et al. [31] and e.g. [32] presents some inconsistencies. For example, it predicts a non-null viscous dissipation for a droplet that moves following a pure rigid rotation. Since our analysis is mainly focused on the energy balance of the droplet dynamics, the formulation of Monaghan & Gingold is preferred in order to avoid possible inconsistencies. Incidentally, we highlight that the above viscous model was widely validated in the context of free-surface flows (see e.g. [16,33,34]).

When the Reynolds number related to the droplet dynamics is larger than ten, acoustic waves may be generated inside the fluid domain as a consequence of the weak-compressibility assumption. This issue is particularly evident when the initial configuration of the free surface includes singularities (like corners or edges). Under these circumstances the use of the acoustic damping term described in [25] proves to be advantageous. This term reads:

$$F_i^{ad} = \alpha_2 \rho_0 c_0 h \sum_j (\nabla \cdot \mathbf{u}_i + \nabla \cdot \mathbf{u}_j) \nabla_i W_{ij} V_j, \qquad (11)$$

where $\alpha_2$ is a dimensionless parameter. As discussed in Section 3.2, the choice of $\alpha_2$ is important, since for $\alpha_2 > 1$ a further constraint on the time step has to be considered for the time integration of the numerical scheme (see Section 3.2).

### 3.1.4. Continuum Surface Force (CSF) model

The specific model for surface tensions adopted in the RHOD-SPH is the Continuum Surface Force (CSF) model described in [11]. This is reported below:

$$\left\{ \begin{aligned} F_i^\sigma &= \sigma \kappa_i \mathbf{n}_i \delta_{Fi}, & \delta_{Fi} &= 2 \left\| \sum_j \nabla_i W_{ij} V_j \right\| \\ \nabla \lambda_i &= \mathbb{L}_i \sum_j (\lambda_j - \lambda_i) \nabla_i W_{ij} V_j, & \mathbf{n}_i &= - \frac{\nabla \lambda_i}{\| \nabla \lambda_i \|}, \\ \kappa_i &= \mathbb{L}_i \sum_j (\mathbf{n}_i - \mathbf{n}_j) \cdot \nabla_i W_{ij} V_j, & \mathbb{L}_i := & \left[ \sum_j (\mathbf{r}_{ji} \otimes \nabla_i W_{ij}) V_j \right]^{-1}. \end{aligned} \right. \qquad (12)$$

where the improvement in the evaluation of the normal vector through the minimum eigenvalue $\lambda$ of the tensor $\mathbb{L}^{-1}$ is described in [35]. Following [11], four critical corrections for thin jets or small drops are also added in order to increase accuracy and robustness of the scheme.

## 3.2. Time integration schemes

A fourth-order Runge-Kutta scheme is adopted for integrating in time the system (7), since this is computationally efficient and allows for larger time steps with respect to other classical integration schemes. As discussed in [36,37], another advantage related to the use of high-order Runge-Kutta schemes is the possibility, classically adopted in other numerical methods, to partly freeze the update of the time-evolved quantities.

The explicit nature of the time integrator imposes a stability criterion on the time step, $\Delta t$, which is obtained as the minimum among all propagation times for the physical information. In particular, this leads to:

$$\left\{ \begin{aligned} \Delta t_a &= 0.25 \min_i \sqrt{\frac{h}{\|\mathbf{a}_i\|}}, & \Delta t_v &= 0.125 \left(\frac{h^2}{v}\right), & \Delta t_\sigma &= 0.1 \max_i \sqrt{\frac{\rho_0 h^2}{\sigma |\kappa_i|}}, \\ \Delta t_c &= 1.0 \left(\frac{h}{c_0}\right), & \Delta t_{ad} &= \frac{\Delta t_c}{\alpha_2}, & \Delta t &= \min(\Delta t_a, \Delta t_v, \Delta t_\sigma, \Delta t_c, \Delta t_{ad}), \end{aligned} \right. \qquad (13)$$
where $||a_i||$ is the particle acceleration and $\kappa_i$ the local free-surface curvature as evaluated in Eq. (12). The constraint on the maximum particle acceleration is relevant in problems involving liquid impacts, whereas in the present work it can be always neglected. The second and the third constraints, namely $\Delta t_v$ and $\Delta t_\sigma$, are linked to the viscous diffusion and to the surface tension modelling [31,38].

The time step $\Delta t_c$ is related to the weak compressibility of the medium through the speed of sound $c_0$, whereas $\Delta t_{ad}$ is the time step associated to the use of the acoustic damping term [25]. The latter one can be, however, discarded when $\alpha_2 < 1$. It is also important to highlight that:

(i) if the speed of sound, $c_0$, is evaluated using all the constraints in (4), the formula for $\Delta t_c$ already embeds that for $\Delta t_\sigma$;

(ii) for computational convenience, one may choose to evaluate $c_0$ by using only the constraints in (4) that depend on $\Delta U_{max}$ and $\Delta P_{max}$. This might lead to the occurrence of local compressible effects in the fluid regions close to high curvatures of the free-surface. Such an approach is, however, allowed as long as these effects play a minor role (see the discussion in Section 4.1). Otherwise, it is necessary to consider the time-step constraint $\Delta t_\sigma$ to avoid numerical instabilities.

### 3.3. SPH energy balance in the presence of surface tension

The energy balance of the RHOD-SPH with surface tension is obtained by following the approach described in [18] and, later, in [23]. It reads:

$$ 
\left\{ 
\begin{aligned}
\dot{\mathcal{E}}_K &= P_{diss} + P_\sigma, && \mathcal{E}_K = \frac{1}{2} \sum_i m_i (\mathbf{u}_i \cdot \mathbf{u}_i), \\
P_\sigma &= \sum_i (F_i^\sigma \cdot \mathbf{u}_i), && P_{diss} := P_v + P_N. \\
P_v &= \sum_i (F_i^v \cdot \mathbf{u}_i),
\end{aligned}
\right. 
\qquad (14) 
$$

where $P_\sigma$ is the power of the surface tension forces while $P_{diss}$ is the overall power dissipated by the scheme. This includes the power dissipated by the viscous forces, $P_v$, and the power dispelled by the numerical dissipation, $P_N$. The latter is linked to the adopted scheme and is expected to converge to zero in the limit of $\Delta x \to 0$. The work exerted by the surface tension and the overall work dissipated by the scheme are defined as:

$$ W_\sigma = \int_{t_0}^{t} P_\sigma d\tau, \quad Q_{diss} = \int_{t_0}^{t} P_{diss} d\tau \qquad (15) $$

It is worth noting that, even if $W_\sigma$ should be theoretically equal to $(\mathcal{E}_{\sigma 0} - \mathcal{E}_\sigma)$ according to Eq. (5), this cannot be ensured numerically. In particular the reversibility of $W_\sigma$ cannot be guaranteed in the adopted SPH model. As a matter of fact, the work $W_\sigma$ is intrinsically related to the CSF model (see Section 3.1.4) which is based on the evaluation of local quantities (like the normal and the curvature of the free-surface shape), without any direct connection with global quantities (like the surface energy associated to the deformation of $\partial \Omega_F$).

### 3.4. Numerical procedure for verifying the surface tension energy consistency

In order to check the consistency of the SPH scheme in the evaluation of $W_\sigma$, we compute $\mathcal{E}_\sigma$ through a Level-Set function $\phi$ and monitor the deviations of the numerical outputs from the identity $W_\sigma = (\mathcal{E}_{\sigma 0} - \mathcal{E}_\sigma)$.

The level-set function is obtained using the algorithm described in [39] based on interpolation of the particles set on a Cartesian mesh. Fig. 1 depicts an example of the level-set function for a 2D liquid square droplet (see Section 4.1). After the contour level $\phi = 0$ is extracted, the length $L_F$ of the free-surface is evaluated, and the surface energy is computed as $\mathcal{E}_\sigma = \sigma L_F$. It is important to point out that the level-set function is evaluated in a post-processing stage and for a limited number of time instants. Because of the CPU costs, it is not calculated at runtime for each time iteration.

At the end of each simulation we check whether the work $W_\sigma$ is consistent with $\mathcal{E}_\sigma$ by measuring the error:

$$ E_{rr\sigma} = \frac{1}{(t_{end} - t_0)} \int_{t_0}^{t_{end}} \left| W_{\sigma}(t) - \left[ \mathcal{E}_{\sigma 0} - \mathcal{E}_{\sigma}(t) \right] \right| dt \qquad (16) $$

The error $E_{rr\sigma}$ is evaluated for different spatial resolutions in order to check if it converges toward zero when increasing $N = L/\Delta x$.

Note that, the surface energy $\mathcal{E}_\sigma$ is always evaluated geometrically by the particles positions through the Level-Set function $\phi$ whereas the work of the surface tension forces $W_\sigma$ is evaluated from the CSF model which is embedded in the SPH scheme.

## 4. Numerical test cases

### 4.1. Oscillation of a 2D liquid square droplet

The benchmarks presented in this section and in the next one are among the most used for testing surface tension models (see e.g. [40–51]). At the initial instant the fluid bulk is a square with side $L$, with null pressure and velocity fields. During the time evolution the initial surface energy, i.e. $\mathcal{E}_{\sigma 0} = 4L\sigma$, is converted into kinetic energy resulting in an oscillatory deformation of the
<figure><img src="image_4.png" /><figcaption>Fig. 1. Top: Level-set function $\phi$ contour calculated during a RHOD-SPH simulation of 2D liquid square droplet (see Section 4.1). Bottom: enlarged view showing the Cartesian mesh adopted for the evaluation of $\phi$ and the normal vectors $n_i$ on $\partial\Omega_F$ calculated through Eq. (12). The red dash-dotted curve represents the $\phi = 0$ contour level.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. Oscillation of a 2D liquid square droplet: Pressure field in the initial and final instants of the simulations. Dash-dotted circle is the final drop configuration expected at the end of the simulation.</figcaption></figure>

drop. This motion is gradually damped by the viscous forces and by the numerical dissipation. At $t \to \infty$ the drop becomes circular with a radius $R = L/\sqrt{\pi}$ with a constant pressure field equal to the Laplace pressure:

$$P_0 = \sigma/R = \sqrt{\pi}\sigma/L$$

Fig. 2 depicts the initial and final states of the simulation. The energy gap:

$$\Delta\mathcal{E} = \mathcal{E}_{\sigma 0} - \mathcal{E}_{\sigma}(\infty) = 2\sigma L (2 - \sqrt{\pi}), \qquad (17)$$

drives the flow evolution whereas the magnitude of the viscous forces influences the damping of the oscillation cycles.

If we assume that the energy gap $\Delta\mathcal{E}$ is converted in kinetic energy according to the relation below:

$$\Delta\mathcal{E} = \frac{1}{2} \rho_0 L^2 U_{ref}^2,$$

the reference velocity can be represented as:

$$U_{ref} = \sqrt{\frac{2\Delta\mathcal{E}}{\rho_0 L^2}} = K \sqrt{\frac{\sigma}{L \rho_0}}, \quad K = 2\sqrt{(2 - \sqrt{\pi})}$$

and, being $K \approx 1$, we approximate it as $U_{ref} = \sqrt{\sigma/(L\rho_0)}$ without any loss of generality. As a consequence of this choice, the Reynolds and Weber numbers become:
<caption>**Table 1**</caption>

Oscillation of a 2D liquid square droplet: numerical evaluation of the maximum curvature $|κ|_{max}$ at initial configuration. Calculation are performed for different spatial resolutions, namely $N = L/Δr = 25, 50, 100, 200$, and for three different ratios $h/Δr$.

<table><thead><tr><th rowspan="2">N</th><th colspan="3">h/Δr = 2.82</th><th colspan="3">h/Δr = 2.00</th><th colspan="3">h/Δr = 1.49</th></tr><tr><th>|κ|<sub>max</sub> L</th><th>L/h</th><th>|κ|<sub>max</sub> h</th><th>|κ|<sub>max</sub> L</th><th>L/h</th><th>|κ|<sub>max</sub> h</th><th>|κ|<sub>max</sub> L</th><th>L/h</th><th>|κ|<sub>max</sub> h</th></tr></thead><tbody><tr><td>25</td><td>9.33</td><td>8.86</td><td>1.05</td><td>13.9</td><td>12.5</td><td>1.11</td><td>19.1</td><td>16.8</td><td>1.14</td></tr><tr><td>50</td><td>18.7</td><td>17.7</td><td>1.05</td><td>27.9</td><td>25.1</td><td>1.11</td><td>38.2</td><td>33.5</td><td>1.14</td></tr><tr><td>100</td><td>37.3</td><td>35.5</td><td>1.05</td><td>55.7</td><td>50.1</td><td>1.11</td><td>76.5</td><td>67.0</td><td>1.14</td></tr><tr><td>200</td><td>74.6</td><td>79.9</td><td>1.05</td><td>111</td><td>100</td><td>1.11</td><td>153</td><td>134</td><td>1.14</td></tr><tr><td>400</td><td>149</td><td>142</td><td>1.05</td><td>223</td><td>201</td><td>1.11</td><td>306</td><td>268</td><td>1.14</td></tr></tbody></table>

$$ \text{Re} = \frac{\rho_0 U_{ref} L}{\mu} = \frac{\sqrt{\rho_0 \sigma L}}{\mu}, \quad \text{We} = \frac{\rho_0 U_{ref}^2 L}{\sigma} = 1, $$

hence, only the Reynolds number remains as a free dimensionless parameter. Note that this depends on the ratio $\sqrt{\sigma}/\mu$ which is the only non geometrical parameter. In the following simulations, five different Reynolds numbers, namely $Re = 3, 10, 50, 100$ and $500$, are considered.

Since at the initial time the droplet has a squared shape, the fourth oscillation mode is the mainly excited mode during the time evolution. It follows that for the time scale, the period of oscillation is expected to be close to the analytical estimation given in [11] for $n = 4$:

$$ T_{(n,\sigma)} = \frac{2\pi}{\sqrt{n(n+1)(n-1)}} \sqrt{\frac{\rho_0 R^3}{\sigma}} \Rightarrow T_{4,\sigma} = \sqrt{\frac{\sqrt{\pi}}{15}} \frac{L}{U_{ref}}, \quad (18) $$

The above relation is valid for oscillations characterized by small amplitudes, while for large amplitudes, as in the present case, the recorded period is expected to be slightly larger during the early stages of the evolution.

In particular, for the simulations with $Re = 100$ the period of the first oscillation is $T \approx 1.05 T_{4,\sigma}$.

The initial time of the simulation is arbitrarily set equal to $t_0 = -T/4$, so that the maximum of the pressure and kinetic energy at the drop centred is attained at about $t = 0$.

### 4.1.1. Choice of the speed of sound

For the present benchmark the speed of sound should be, in principle, set through the constraints of Eq. (4). This would, however, imply that an infinite speed of sound should be required at $t = t_0$, as a consequence of the singularities at the square corners. In any case, being consistent with the smoothing procedure at the foundation of SPH, the maximum allowable/modelling curvature in the numerical simulation can be assumed equal to $|κ|_{max} \approx 1/h$. A heuristic proof for the latter statement is given in Table 1, where the numerical evaluation of $|κ|_{max}$ in the initial squared drop configuration is performed varying the spatial resolution $N = L/Δr$ for three different ratios $h/Δr$.

Substituting $|κ|_{max} \approx 1/h$ in Eq. (4), we obtain:

$$ c_0 = 10 \sqrt{\frac{\sigma |κ|_{max}}{\rho_0}} = 10 U_{ref} \sqrt{|κ|_{max} L} \approx 10 U_{ref} \sqrt{N/2}, \quad (19) $$

Since the maximum resolution is $N = 400$, it follows that $c_0 \approx 140 U_{ref}$, that is more than ten times larger than $10 U_{ref}$. The related weakly-compressible pressure limit is $\pm 111 P_0$ and these large pressure values are only reached during the initial iterations.

Fig. 3 depicts the initial evolution of the pressure field using two different speeds of sound, namely $c_0 = 140 U_{ref}$ and $c_0 = 50 U_{ref}$. As a consequence of the change of $c_0$, the following observations can be made:

(i) only small differences in the highest pressure regions are detectable (see the enlarged views on the bottom plots of Fig. 3);

(ii) a negligible effect on the free-surface evolution is observed;

(iii) at the indicated time instants, the pressure levels have already dropped in the range of $\pm 10 P_0$;

(iv) for both values of the speed of sound, the singular initial condition leads to the generation of non negligible acoustic-wave systems;

(v) in the initial stages of the evolution both simulations have the same time step, namely $\Delta t_σ$, since this is the most restrictive constraint and does not depend on $c_0$. An example is shown in Fig. 3 where for both the sound speeds the evolution at $t = t_0 + 0.028T$ correspond to 30 time iterations. For longer times, however, $\Delta t_c$ becomes dominant and, consequently, the smaller sound speed implies a reduction of the CPU costs.

The above results show that a smaller speed of sound can be adopted as long as the compressible effects are limited in time and space.

To the authors' knowledge, no reference in the state-of-the-art literature has addressed the potential issue related to the speed of sound and to the fulfilment of the weak-compressibility assumption in presence of free-surface geometrical singularities, as in this benchmark.
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Oscillation of a 2D liquid square droplet: Top: Pressure field in the initial stage of the evolution for the case at Re = 100 using two different speeds of sound, $c_0 = 50 U_{ref}$ (left) and $c_0 = 141 U_{ref}$ (right) Bottom: enlarged views of the pressure field in the right top corner.</figcaption></figure>

### 4.1.2. Damping the acoustic wave components

As discussed in the previous subsection, an acoustic waves system is generated during the early stages of the evolution of the present test-case. The free-surface singularities induce intense free-surface forces which, in turn, generate high pressure regions and lead to the formation of acoustic waves.

Fig. 4 shows the time histories of the pressure recorded at the centre of the droplet for $c_0 = 140 U_{ref}$ (top panel) and $c_0 = 50 U_{ref}$ (middle panel) without the acoustic damping term $F_{ad}$ described in Eq. (11). For the largest value of the sound speed, a high frequency mode is well visible, while for $c_0 = 50 U_{ref}$ lower and more irregular frequency components are detected. The bottom plot of Fig. 4 shows the time histories of the pressure recorded at the centre of the droplet for both values of the speed of sound when the acoustic damping term is used. For $t > 0.2T$ the pressure signals are practically the same, even though for $c_0 = 140 U_{ref}$ the removal of the acoustic component is more effective (see [52] for more details).

It is worth noting that large values of $\alpha_2$ (in this case $\alpha_2 = 6$) are necessary to damp out the acoustic components completely. Consequently, for these simulations the more restrictive time step is $\Delta t_{ad}$. Incidentally, we highlight that in all the four simulations the free-surface configurations (not shown here for the sake of brevity) are practically unchanged, confirming that the acoustic wave system does not alter the droplet dynamics.

Finally, Fig. 5 shows the time history of the kinetic energy without using the acoustic damper (i.e. $\alpha_2 = 0$) and with the acoustic damper with $\alpha_2 = 6$. The two curves are practically superimposed. This behaviour confirms the observations in Sun et al. [25] where it is highlighted that the acoustic damper does not act on the “incompressible flow component” but only on the acoustic component. As a consequence $\alpha_2$ has negligible effect on the kinetic energy of the present test-case.

### 4.1.3. Droplet time evolution at $re = 100$: convergence study

The case Re = 100 discussed in this section corresponds to a water droplet of about $L = 0.1$ mm. As shown in the left plot of Fig. 2, at $t = t_0$ the drop has a squared shape and, consequently, the surface tension is singular at the four vertices. As displayed in the top-left plot of Fig. 6, the vertices are suddenly rounded inducing an high pressure region in their neighbourhood. In the same
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Oscillation of a 2D liquid square droplet, time histories of the pressure recorded at the droplet centre. Top: $c_0 = 140 U_{ref}$ without acoustic damper. Middle: $c_0 = 50 U_{ref}$ without acoustic damper. Bottom: with acoustic damper $\alpha_2 = 6$, comparison for both the speeds of sound $c_0 = 50 U_{ref}$ and $c_0 = 140 U_{ref}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Oscillation of a 2D liquid square droplet, time histories of the kinetic energy without and with the acoustic damper $\alpha_2 = 6$.</figcaption></figure>

picture it is possible to appreciate the absence of acoustic waves which are suppressed thanks to the acoustic damping term (11) (see Section 4.1.2). At $t = t_0 + T/4$ (top-right panel of Fig. 6) the drop has an almost circular shape with the pressure core which reaches its maximum of about $2 P_0$, while at $t = t_0 + T/2$ the droplet elongation reaches a second maximum (bottom-left plot). In the subsequent time steps the drop keeps oscillating with rounded shapes. Finally, at $t = t_0 + 10T$ (right plot of Fig. 2) the 99% of the kinetic energy has been damped and the pressure field is uniform at the Laplace pressure level $P_0$.

The plot of Fig. 7 displays the time histories of the kinetic energy for the five different spatial resolutions adopted. The convergence of the solution is well visible; within the first oscillation a rate of convergence close to 2 is measured, as expected with the RHOD-SPH scheme. For longer times the signals get out of phase, thus reducing the convergence rate. Fig. 8 shows the time histories of the viscous and numerical dissipation, namely $Q_V$ and $Q_N$. Increasing the spatial resolution, both quantities converge monotonously in all the time range of the simulations. A convergence rate close to one is measured for both the dissipation quantities. It is worth highlighting that only for $N > 200$ the numerical dissipation decreases below 10% of the energy gap $\Delta E$. This suggests that a good energy balance is only possible for high spatial resolutions.

The left plot of Fig. 9 shows the time histories of the pressure recorded at the centre of the drop for $Re = 100$ and for five different spatial resolutions, namely $N = 25, 50, 100, 200, 400$. For this case more than 12 periods are required to damp the pressure field oscillations. The left plot of Fig. 9 clearly shows that the solution converges towards the Laplace pressure $P_0$ at the end of the simulation. The convergence rate is between one and two. The maximum pressure recorded is about $2.4 P_0$, that is rather below the weakly-compressible pressure limit. Doing a comparison with the approximate solution described in [11], this case indicatively corresponds to a motion amplitude of about $a \approx 0.02R$ for the four-mode oscillating drop.

The right plot of Fig. 9 depicts the time histories of the work $W_\sigma$ performed by the surface tension $f_\sigma$ for the case $Re = 100$, using the five different spatial resolutions indicated before. Again, the convergence rate is between one and two. It is evident that at the end of the simulations the curves converge towards the correct analytical value, namely $\Delta E$.

Finally, the left plot of Fig. 10 shows the time histories of the surface energy $E_\sigma$ measured via a Level-Set function as described in Section 3.4. Similarly to what shown for the work of the surface tension force, the surface energy also converges. However, the rate of convergence is close to one, that is lower than the one measured for the work $W_\sigma$.

To verify the coherence between the work $W_\sigma$ and the surface energy $E_\sigma$, the error $E_{err\sigma}$ defined in the Eq. (16) is evaluated for the five spatial resolutions considered (see the right plot of Fig. 10). The error $E_{err\sigma}$ converges with a rate between one and two.
<figure><img src="image_9.png" /><figcaption>Fig. 6. Oscillation of a 2D liquid square droplet: Pressure field evolution for the case at Re = 100. Dash-dotted circle is the final drop configuration expected at the end of the simulation. The video of the simulation is available at [Link Video N°1](https://www.youtube.com/watch?v=dQw4w9WgXcQ).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. Oscillation of a 2D liquid square droplet. Time histories of the kinetic energy $E_K$. Five different spatial resolutions are adopted, i.e. $N = 25, 50, 100, 200, 400$.</figcaption></figure>

### 4.1.4. Inconsistency and convergence issues when using more classical SPH schemes

It is relevant to underline that the RHOD-SPH scheme is expected to be a high-order SPH scheme with convergence rate close or above second order (see [23]). Here we showed that, when the CSF surface tension model is included, the order of the RHOD-SPH is reduced. This is due to the intrinsic nature of the CSF model: the surface tension forces are evaluated only using the particles in the free-surface region, a region where the particle spatial arrangement is difficult to keep regular (the PST can only act along the tangential directions in this region). In this subsection we show that SPH schemes that do not implement the PST provide outputs that are even worse.

The simulations for the test-case at Re = 100 at five different spatial resolutions were repeated using three different SPH schemes: the Standard SPH, the $\delta$-SPH and the Riemann-SPH. All these three models can be obtained by the RHOD-SPH by switching off some of the corrections described in Section 3.1. All these SPH schemes do not use the PST (i.e. $\delta u = 0$) and the tensile instability is not prevented. On the other hand, conversely to the RHOD-SPH, those schemes are conservative in terms of linear and angular momenta.

In the $\delta$-SPH a numerical diffusion is adopted and inserted in the continuity equation, instead of using Riemann solvers in both continuity and momentum equations (see [36]). In the Standard SPH the numerical diffusion corrections are not adopted at all.
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Oscillation of a 2D liquid square droplet. Top: time histories of the viscous dissipation $Q_v$. Bottom: time histories of the numerical dissipation $Q_N$. Five different spatial resolutions are adopted, i.e. $N = 25, 50, 100, 200, 400$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Oscillation of a 2D liquid square droplet using five different spatial resolutions, namely $N = 25, 50, 100, 200, 400$. Left: time histories of the pressure recorded at the centre of the drop. Right: Time histories of the work performed by the surface tension forces, $W_σ$.</figcaption></figure>

As well documented in literature, the consequence is that in this scheme the pressure field is generally affected by non-physical high-frequency noise.

Fig. 11 shows the error $E_{rr\sigma}$ as a function of the spatial resolution $N$ for the Standard SPH, the $\delta$-SPH and the Riemann-SPH. For all these three models the convergence is not attained. The reason is related to the particles spatial distribution which becomes locally irregular when increasing the spatial resolution, introducing large errors in the evaluation of the normal vector field and, consequently, in the surface tension forces. Fig. 12 displays a snapshot obtained with the Standard SPH scheme where pressure instabilities develops in some specific regions where particle spatial distribution becomes too irregular. Similar behaviours are observed with the $\delta$-SPH and Riemann-SPH schemes but they are not reported here for the sake of brevity.

To the authors' knowledge, the above convergence issues related to the use of the CFS model in the SPH schemes are not well documented in the literature. This is likely due to the fact that the present test-case is usually solved only for very low Reynolds numbers, leading to a weaker dynamics of the free-surface.
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Oscillation of a 2D liquid square droplet using five different spatial resolutions, namely $N = 25, 50, 100, 200, 400$. Left: time histories of the surface energy $E_σ$ measured through the level-set function (see Section 3.4). Right: the error $E_{rr,σ}$ (see Eq. (16)) as a function of the spatial resolution $N$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Oscillation of a 2D liquid square droplet using five different spatial resolutions, namely $N = 25, 50, 100, 200, 400$. The error $E_{rr,σ}$ (see Eq. (16)) as a function of the spatial resolution $N$ for the schemes: SPH Standard, δ-SPH and Riemann-SPH.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Oscillation of a 2D liquid square droplet at $Re = 100$. Left: pressure field at time instant $t = 0.471T$ evaluated with the Standard SPH using the resolution $N = 400$. Right: enlarged view showing the particle spatial distribution and the free-surface normal vectors.</figcaption></figure>

### 4.1.5. Droplet time evolution for five different Reynolds numbers

Fig. 13 contains the time histories of the energy balance for $Re = 3, 10, 50, 100$ and $500$ at the highest resolution, i.e. $N = 400$. In all the cases a good agreement between $(E_{σ0} - E_σ)$ and $W_σ$ is attained. For the lowest Reynolds numbers, i.e. $Re = 3$ and $Re =$
<figure><img src="image_16.png" /><figcaption>Fig. 13. Oscillation of a 2D liquid square droplet. Time histories of the kinetic $\mathcal{E}_K$, surface energy $\mathcal{E}_\sigma$, dissipated energy $Q_{diss}$ and the work of the surface tension forces $W_\sigma$ (in the legend $W_\sigma$ is denoted as $W_{st}$), for the Reynolds numbers $Re = 100, 50, 10$ and 3.</figcaption></figure>

10, the kinetic energy, $\mathcal{E}_K$, is almost completely damped within the first period of oscillation and the maximum occurs before the time $t = 0$.

The highest Reynolds number considered for this benchmark is $Re = 500$. The top plot of Fig. 14 displays the time histories of the energy balance at the highest resolution, namely $N = 400$. On the bottom plot of the same figure the dashed-dotted line represents the damping predicted by the linear theory for a viscous drop in zero-gravity oscillating under surface tension with mode $n = 4$ (see [24]). Since this test is initially highly nonlinear, the linear solution is not valid during the first oscillations, while for $t > 2T$ a good matching between the RHOD-SPH and the analytical solution is obtained.

It is worth noting that, the complex behaviour of the energy components displayed in Fig. 14 in comparison to the time series discussed for the lower Reynolds numbers is an indication on how the droplet dynamics becomes more complicated when the Reynolds number is increased. This is further confirmed by comparing the evolution of the pressure field shown in Fig. 15 to that in Fig. 6 for a lower Reynolds number. Generally, for $Re \ge 100$ a suitable numerical scheme is required in order to preserve the coherence between the surface energy $\mathcal{E}_\sigma$ and the work of the surface forces $W_\sigma$ given by the specific CFS model adopted. For this reason, the test-cases provided in the present work can be useful for testing other numerical schemes.

## 4.2. Oscillation of a 3D liquid cubic droplet

The analysis performed in a 2D framework in the previous section is replicated here for a cubic droplet.

In the initial condition the liquid drop has a cubic shape with side length $L$, with null pressure and velocity fields. The initial energy of the surface tension is $\mathcal{E}_\sigma(t_0) = 6L^2\sigma$ and, during the evolution, it is converted in kinetic energy and damped by the viscous forces and the numerical dissipation.
<figure><img src="image_17.png" /><figcaption>Fig. 14. Oscillation of a 2D liquid square droplet for Reynolds numbers Re = 500. Top: time histories of the kinetic $\epsilon_K$, surface energy $\epsilon_\sigma$, dissipated energy $Q_{diss}$ and the work of the surface tension forces $W_\sigma$ (in the legend $W_\sigma$ is denoted as $W_{st}$). Bottom: time histories of the kinetic $\epsilon_K$ compared with the energy decay predicted by the linear theory [24].</figcaption></figure>

For $t \to \infty$ the drop becomes a sphere with radius $R = L\sqrt{3/(4\pi)}$ and constant pressure field equal to the Laplace pressure $P_0 = 2\sigma/R$. For this test-case the energy gap is:

$$ \Delta\mathcal{E} = \mathcal{E}_{\sigma}(t_0) - \mathcal{E}_{\sigma}(\infty) = 2\sigma L^2 \left[ 3 - 2\pi \left( \frac{3}{4\pi} \right)^{2/3} \right] \quad (20) $$

It follows that the reference velocity for the problem is:

$$ U_{ref} = \sqrt{\frac{2\Delta\mathcal{E}}{\rho_0 L^3}} = K \sqrt{\frac{\sigma}{L \rho_0}}, \quad K = 2\sqrt{3 - 2\pi \left(\frac{3}{4\pi}\right)^{2/3}}, $$

Since $K \approx 1.5$, the reference velocity $U_{ref} = \sqrt{\sigma/(L\rho_0)}$ can still be used, hence, the definitions of the Reynolds and Weber numbers are the same of the 2D case. In the following analysis only Re = 10 is considered in order to limit the time range of the simulations and reduce the computational costs. In the numerical simulations, the measured oscillation period of the drop is $T = 1.22 T_{4,\sigma}$ and, again, the initial time is set equal to $t_0 = -T/4$.

In Fig. 16 the pressure field for six different time instants is shown. Analogously to the 2D test-case, at $t = t_0$, the surface tension is singular at the eight vertices of the cube, which are suddenly rounded inducing high pressure zones in the neighbourhood. At $t = t_0 + T/4$ the kinetic energy is close to its maximum value, while in the subsequent instants the droplet shape becomes more and more spherical. For $t > T$ the pressure field is almost uniform and equal to the Laplace pressure level $P_0$.

Different spatial resolutions, namely $N = 12, 25, 50$ and $200$, are examined. As a reference, the finest spatial resolution corresponds to 8 millions particles. As in the 2D test-case, the speed of sound is $c_0 = 140 U_{ref}$, and at the maximum resolution $\Delta t_\sigma \approx \Delta t_c$ with $\Delta t_c = 2 \cdot 10^{-4} T$. Hence, about 11,000 time steps were required.

The top plot of Fig. 17, reports the time histories of the pressure recorded at the drop centre. This plot highlights the convergence towards the Laplace pressure $P_0$ when $t > T$. The energy balance in time is depicted in the bottom plot of Fig. 17 for the highest spatial resolution, i.e. $N = 200$. As in the 2D test-case a good agreement between the energy $\epsilon_\sigma$ (evaluated through the level-set function) and the work of the surface tension force $W_\sigma$ (computed by the SPH scheme) is attained.

The left side of Fig. 18, shows precisely the convergence of $W_\sigma$ for five different resolutions. The right panel of the same figure shows the time histories of the viscous dissipation at different resolutions. A good convergence toward the analytical energy gap $\Delta\mathcal{E}$ is observed.

Fig. 19 displays the relative errors for the pressure at the droplet centre, the work $W_\sigma$ and the viscous dissipation evaluated at the end of the simulations with respect to the analytical values.
<figure><img src="image_18.png" /><figcaption>**Fig. 15.** Oscillation of a 2D liquid square droplet: Pressure field evolution for the case at Re = 500. Dash-dotted circle is the final drop configuration expected at the end of the simulation. The video of the simulation is available at Link Video N°2.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 16.** Oscillation of a 3D liquid square droplet: Pressure field evolution for the case at Re = 10. Only one half of the domain is plotted to show the pressure field inside the drop. The spatial resolution is N = 100.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Fig. 17.** Oscillation of a 3D liquid cubic droplet. Top: time histories of the pressure recorded at the drop centre. Five different spatial resolutions are used, namely *N* = 12, 25, 50, 100, 200. Bottom: time histories of the kinetic energy *E*<sub>k</sub>, surface energy *E*<sub>s</sub>, dissipated energy *Q*<sub>diss</sub> and the work of the surface tension forces *W*<sub>s</sub> (in the legend *W*<sub>s</sub> is denoted as *W*<sub>st</sub>).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 18.** Oscillation of a liquid cubic droplet. Left: time histories of surface tension work. Right: time histories of the viscous dissipation and of the work of the surface tension force. Five different spatial resolutions are used, namely *N* = 12, 25, 50, 100, 200.</figcaption></figure>

Finally, Fig. 20 depicts the pressure field evolution on the free-surface in the initial stage for the case with $Re = 50$. A capillary wave system originates from the twelve edges of the cube, propagates and focuses close to the centres of the cube faces. This result shows that the increase of the Reynolds number leads the free-surface configurations to evolve in a rather complex way. Consequently, high spatial resolutions are needed to capture the propagation of the capillary wave system. For this reason, much
<figure><img src="image_22.png" /><figcaption>Fig. 19. Oscillation of a liquid cubic droplet. Relative errors at the end of the simulations with respect to the analytical values for the: pressure measured at the centre of the droplet (left), work of the surface tension forces (middle), viscous dissipated energy (right). Five different spatial resolutions are used, namely N = 12, 25, 50, 100, 200.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 20. Oscillation of a 3D liquid square droplet: Pressure field associated to the capillary waves generated by the 12 edges. The spatial resolution is N = 200, Reynolds number Re = 50.</figcaption></figure>

higher spatial resolutions are needed to resolve the Reynolds number range investigated in 2D with a considerable increase of the CPU costs. This topic is, however, left open for future works.

For the cases studied in this section, due to the low Reynolds numbers considered, $\alpha_2 = 1$ can be adopted for the acoustic damper (11) with no restrictions on the time steps.

### 4.3. Large amplitude oscillation of a 2D liquid drop

In this section we simulate the evolution of a viscous drop in zero-gravity oscillating under surface tension with mode $n = 2$. In the numerical simulations we assume that the drop is initially circular with a drop radius equal to R. Following the approach by Lamb [53], the fluid velocity is approximated through the velocity potential $\phi$. Specifically following [24]), for $\sigma > 8 \rho_0 v^2/(3R)$, the solution exhibits a damped oscillatory behaviour with a solution in polar coordinates given by:

$$ \phi(r, \theta, t) = \frac{a \omega R}{2} \left(\frac{r}{R}\right)^2 \cos(2\theta) \cos(\omega t + \psi) \exp(-\beta t), \quad (21) $$

where a is the amplitude of the motion and $\psi$ is a proper phase shift that is defined later. The above analytical solution is obtained through a perturbation expansion based on the assumption that the parameter $\epsilon = a/R$ is small (i.e. $\epsilon \ll 1$). The damping coefficient $\beta$ and the angular frequency $\omega$ are given by:

$$ \beta = 4 \frac{v}{R^2}, \quad \omega = \sqrt{\frac{2}{R^3} \left( \frac{3\sigma}{\rho_0} - \frac{8v^2}{R} \right)}, \quad T_{2,\sigma} = \frac{2\pi}{\omega}, \quad (22) $$
where $T_{2,σ}$ is the oscillation period. The velocity components for the damped-oscillatory regime are:

$$ egin{cases} u_r(r, \theta, t) = a \omega \left(\frac{r}{R}\right) \cos(2\theta) \cos(\omega t + \psi) \exp(-\beta t) \\ u_\theta(r, \theta, t) = -a \omega \left(\frac{r}{R}\right) \sin(2\theta) \cos(\omega t + \psi) \exp(-\beta t). \end{cases} \quad (23) $$

The solution for the free surface is:

$$ r_{FS} = R + \frac{a\omega R}{(\omega^2 + \beta^2)} \cos(2\theta) [\omega \sin(\omega t + \psi) - \beta \cos(\omega t + \psi)] \exp(-\beta t). \quad (24) $$

Imposing the drop to be initially circular, we find $\psi = \arctan(\beta/\omega)$.

The expression for the pressure field is not provided in [24] and is derived here by substituting the approximate solution for the velocity field in the Bernoulli equation. According to Eq. (3), the leading order of the pressure field along the free surface is imposed to be equal to the Laplace pressure, i.e. $P_0 = \sigma/R$:

$$ \frac{p(r, \theta, t)}{\rho_0} = \frac{P_0}{\rho_0} + \frac{a \omega R^2}{2} \left(\frac{r}{R}\right)^2 \cos(2\theta) \left[ \beta \cos(\omega t) + \omega \sin(\omega t) \right] \exp(-\beta t) + \frac{a^2 \omega}{2 R^2} \left[ 1 - \left(\frac{r}{R}\right)^2 \right] \cos^2(\omega t + \psi) \exp(-2\beta t). \quad (25) $$

Note that the pressure field along the circumference at the initial time is not identically equal to $P_0$ because of the presence of the viscous stress at the free surface. This latter term, however, is not solved exactly, since the procedure described in [53] just provides an approximation of the dynamic boundary condition (3).

The above analytical solution is valid only for small amplitudes, i.e. $\epsilon = a/R \ll 1$. In any case, in the present section we consider a large amplitude condition, namely $\epsilon = a/R = 1$, so that a highly non-linear behaviour is expected during the early stages of the evolution. The comparison of the energy dissipation ratio to the analytical value $\beta$ is still possible on the last oscillation periods where a linear behaviour is attained.

Eqs. (23) and (25) evaluated at $t = 0$ are used to get the initial conditions for the SPH simulation. To set the simulation, a drop with radius $R = 10.5$ mm, density $\rho_0 = 1000$ kg/m³, kinematic viscosity $\nu = 1.37 \cdot 10^{-5}$ m²/s and surface tension $\sigma = 72$ mN/m is considered. Since the viscosity is low, Eq. (26) can be approximated as follows:

$$ \omega \approx \sqrt{\frac{6\sigma}{R^3 \rho_0}}, \quad T_{2,\sigma} \approx 2\pi \sqrt{\frac{\rho_0 R^3}{6\sigma}} \approx 0.3253 \text{ s} \quad (26) $$

Specifically, for the present problem the Weber and Reynolds numbers are defined as:

$$ \text{We} = \frac{\rho_0 a^2 \omega^2 R}{\sigma}, \quad \text{Re} = \frac{a \omega R}{\nu}. $$

and, for the chosen parameters, we obtain $We = 6$ and $Re = 155$.

### 4.3.1. Choice of the speed of sound

Differently to the previous test-cases, the present one has no singularities at the free surface. Hence, the estimation $|k|_{max} h \approx 1$ can be relaxed in the constraint (4) and the speed of sound can be set by using the maximum pressure variation estimated by the linear theory, that is about $10 P_0$. It follows that $c_0 = 13 \omega R$ is a good compromise for all the spatial resolutions considered here.

We verify a posteriori that the weakly compressible effects maintain small by monitoring that the variation of the energy associated to compressibility, namely $\Delta E_C$, remains below the 0.2% of the initial kinetic energy $E_{k0}$. Furthermore, the absence of singularities at the free surface (as for the test-cases studied in Sections 4.1 and 4.2) reduces the generation of acoustic waves; it follows that $\alpha_2 = 1$ is sufficient and, therefore, no restrictions on the time step are required (see Section 3.2).

### 4.3.2. Results

Fig. 21 shows the pressure field at six different time instants, where the drop firstly elongates along the x-axis, reaches the maximum stretching around $t = 0.24T$ and goes back.

At $t = 0.49T$ a “pinch-off” configuration occurs and the surface drop gets the maximum shrinkage at $x = 0$ with a thickness of $0.026R$. For the maximum spatial resolution, namely $N = R/\Delta r = 400$, about 13 particles are present in such a small layer of fluid. It is important to underline that in this time instant, during the pinch-off configuration, a pressure drop of about $-4.1P_0$ is reached. Such a negative value can excite tensile instability leading to a non-physical fragmentation. In this case the pressure gradient correction proposed in Section 3.1.2 becomes crucial for this test-case.

After this time instants the drop re-expand in the vertical axis, reaching the maximum elongation on the y-axis at $t = 0.78T$. After one period the drop comes back to a more rounded configuration, even if this deviates from the initial circular shape because of the high non-linearity of the problem.

Fig. 22 sketches the time evolution of the kinetic energy up to $t = 1.25T$. The different curves refers to different spatial resolutions, i.e. $N = 50, 100, 200$ and $400$. With the coarsest resolution it is not possible to resolve the drop shrinkage, since only one particle is present in the small free-surface gap and, consequently, the drop numerically breaks-up because of the tensile instability which cannot be prevented in this condition. For $N = 100, 200, 400$ this non-physical event does not occur and the kinetic energy converges with a rate of about 1.1 evaluated with a $L_1$ norm.
<figure><img src="image_24.png" /><figcaption>**Fig. 21.** Pressure evolution for a large amplitude oscillation of a 2D liquid drop. The dash-dotted circle is the drop configuration in the beginning of the simulation. The spatial resolution is $N = R/\Delta r = 400$. Link Video N°3.</figcaption></figure>

It is worth noting that, if the convergence is evaluated in the time range $t \in [0, 0.45T]$, thus excluding the pinch-off configuration, the rate of convergence rises to 1.5. This is something expected since the pinch off configuration is an almost singular configuration where few particles presented in the liquid thin region influence the subsequent droplet dynamics.

At $t = T$ the second peak of the kinetic energy occurs. The time evolution of the kinetic energy clearly reflects the non-linear behaviour of the drop evolution in its firsts oscillations. This is also confirmed by the difference observed in the oscillation period between the numerical outputs and the linear theory. In fact at the beginning of the simulation $T = 1.77 T_{2,σ}$. Conversely, Fig. 23 reports the time history of the kinetic energy for the last 5 periods of oscillations at the finest spatial resolution (N = 400). The time axis is shifted with the time $t_1$ which corresponds to 2.25 drop oscillations. The time axis is made dimensionless with the oscillation period $T_{2,σ}$ predicted by the linear theory. In this time range the signal has an almost sinusoidal behaviour for with a good matching
<figure><img src="image_25.png" /><figcaption>Fig. 22. Large amplitude oscillation of a 2D liquid drop: initial stages of the evolution of the kinetic energy for different spatial resolutions.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 23. Large amplitude oscillation of a 2D liquid drop: long-time evolution of the kinetic energy and of the dissipated energy $Q_{diss}$. The spatial resolution is N = 400.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 24. Large amplitude oscillation of a 2D liquid drop: long-time evolution of the main energy components for the finest spatial resolution (namely, N = 400). Note that in the legend $W_σ$ is denoted as $W_{st}$.</figcaption></figure>

with the damping predicted by the linear theory, i.e. $\exp(-2\beta t)$, is obtained. A further confirmation of the good matching between the simulation and the linear theory is given by the period registered for the last oscillation which is $T = 1.03 T_{2,σ}$.

Finally, Fig. 24 displays the long-time evolution of the main energy components for the finest spatial resolution (namely, N = 400) and the comparison between the energy of the surface tension evaluated by the SPH scheme and by the Level-Set function. An overall good match is observed all over the considered time window.
## 5. Conclusions

The present study describes an approach to the modelling of the surface tension through a single-phase method embedded in the framework of the SPH scheme. The analysis is performed by studying the evolution of viscous drops oscillating under the action of the surface tension only and inspecting the consistency of the overall energy balance with respect to the contribution of the surface energy. The surface energy is obtained in two different ways: (i) in terms of the work done by the surface tension forces $W_σ$ evaluated with the SPH model presented in [11], (ii) through the computation of the liquid surface S with the use of a level-set function and evaluating the surface energy directly as $E_σ = σS$. Tests in both 2D and 3D are presented, proving for all the cases the accuracy of the selected SPH model.

In particular, the evolution of initially squared and cubic drops are considered, starting with null pressure and velocity fields and without external volume forces. Under these conditions, the droplets start their oscillations and exchange the contribution of surface energy with the kinetic energy, which is gradually damped by both the viscosity and the numerical dissipation. Different Reynolds number are tested in the 2D framework (i.e. $Re = 3, 10, 50, 100$ and $500$) using different spatial resolutions, while in 3D the analysis is limited to $Re = 10$, because of the higher CPU costs.

A detailed convergence analysis is performed for the Reynolds number 100 (2D test-case). Convergence rates for the different energy components and for the local pressure measured at the centre of the droplet are provided. As an important outcome, it is shown that the use of a CFS model reduces the order of convergence of the RHOD-SPH scheme from 2 to an order between 1 and 2. More importantly, convergence cannot be attained when using classical SPH schemes as the Standard, the $δ$- and the Riemann-SPH models. The above findings suggest that, for an accurate and robust modelling of the surface tension in the SPH framework, a number of numerical tools are necessary. In particular, these include the use of the PST, the renormalized pressure gradient and the diffusive terms.

As for the future perspective, the application and validation of the RHOD-SPH model needs to be extended to problems where drop coalescence and breakup phenomena occur. A further extension of the model concerns regimes with low Weber numbers ($We < 1$) or regimes with high Reynolds numbers ($Re > 500$) where numerical issues can arise and enhanced solutions may be required.

## CRediT authorship contribution statement

**S. Marrone:** Writing – review & editing, Writing – original draft, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **M. Antuono:** Writing – review & editing, Writing – original draft, Validation, Supervision, Methodology, Investigation, Formal analysis. **A. Agresta:** Writing – review & editing, Writing – original draft, Validation, Resources, Investigation, Formal analysis. **A. Colagrossi:** Writing – original draft, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Data curation, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

The work was supported by the project "Next Generation SPH schemes for complex multiphase flows" NEOGEO (CUP B83C23003850006) in collaboration with Ecole Centrale de Nantes in the framework of their Chair programme funded by Siemens Digital Industries Software.

The work was also partially supported by the Project BIO-EMBRACE, PRIN PNRR P202298P25 - CUP B53D23026860001, financed by the European Union – Next Generation EU.

## Data availability

Data will be made available on request.

## References

[1] A. Zhang, P. Sun, F. Ming, An sph modeling of bubble rising and coalescing in three dimensions, Comput. Methods Appl. Mech. Engrg. 294 (2015) 189–209.

[2] A.M. Tartakovsky, A. Panchenko, Pairwise force smoothed particle hydrodynamics model for multiphase flow: surface tension and contact line dynamics, J. Comput. Phys. 305 (2016) 1119–1146.

[3] M. Huber, F. Keller, W. Säckel, M. Hirschler, P. Kunz, S.M. Hassanizadeh, U. Nieken, On the physically based modeling of surface tension and moving contact lines with dynamic contact angles on the continuum scale, J. Comput. Phys. 310 (2016) 459–477.

[4] M. Hirschler, G. Oger, U. Nieken, D. Le Touzé, Modeling of droplet collisions using smoothed particle hydrodynamics, Int. J. Multiph. Flow 95 (2017) 175–187.

[5] S. Braun, L. Wieth, S. Holz, T.F. Dauch, M.C. Keller, G. Chaussonnet, S. Gepperth, R. Koch, H.J. Bauer, Numerical prediction of air-assisted primary atomization using smoothed particle hydrodynamics, Int. J. Multiph. Flow 114 (2019) 303–315.
[6] L. Yang, M. Rakhsha, W. Hu, D. Negrut, A consistent multiphase flow model with a generalized particle shifting scheme resolved via incompressible sph, J. Comput. Phys. 458 (2022) 111079.

[7] Y. Xu, G. Yang, S. Liu, D. Hu, Improvement of surface tension discrete model in the isph-fvm coupling method, Int. J. Multiph. Flow 160 (2023) 104347.

[8] J. Brackbill, D. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (1992) 335-354.

[9] B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, G. Zanetti, Modelling merging and fragmentation in multiphase flows with surfer, J. Comput. Phys. 113 (1994) 134-147.

[10] J.P. Morris, Simulating surface tension with smoothed particle hydrodynamics, Int. J. Numer. Methods Fluids 33 (2000) 333-353.

[11] A. Vergnaud, G. Oger, D. Le Touzé, M. De Leffe, L. Chiron, C-CSF: Accurate, robust and efficient surface tension and contact angle models for single-phase flows using SPH, Comput. Methods Appl. Mech. Engrg. 389 (2022) 114292.

[12] A.M. Tartakovský, P. Meakin, A smoothed particle hydrodynamics model for miscible flow in three-dimensional fractures and the two-dimensional rayleigh-taylor instability, J. Comput. Phys. 207 (2005) 610-624.

[13] M. Becker, M. Teschner, Weakly compressible sph for free surface flows, in: Proceedings of the 2007 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, 2007, pp. 209-217.

[14] N. Akinci, G. Akinci, M. Teschner, Versatile surface tension and adhesion for sph fluids, ACM Trans. Graph. 32 (2013) 1-8.

[15] A. Colagrossi, M. Antuono, A. Souto-Iglesias, D. Le Touzé, Theoretical analysis and numerical verification of the consistency of viscous smoothed-particle-hydrodynamics formulations in simulating free-surface flows, Phys. Rev. E 84 (2011) 026705.

[16] A. Colagrossi, A. Souto-Iglesias, M. Antuono, S. Marrone, Smoothed-particle-hydrodynamics modeling of dissipation mechanisms in gravity waves, Phys. Rev. E 87 (2013) 023302.

[17] A. Colagrossi, B. Bouscasse, S. Marrone, Energy-decomposition analysis for viscous free-surface flows, Phys. Rev. E 92 (2015) 053003.

[18] M. Antuono, S. Marrone, A. Colagrossi, B. Bouscasse, Energy balance in the δ-SPH scheme, Comput. Methods Appl. Mech. Engrg. 289 (2015).

[19] J. Michel, M. Antuono, G. Oger, S. Marrone, Energy balance in quasi-lagrangian riemann-based sph schemes, Comput. Methods Appl. Mech. Engrg. 410 (2023a) 116015.

[20] S. Marrone, A. Colagrossi, A. Di Mascio, D. Le Touzé, Prediction of energy losses in water impacts using incompressible and weakly compressible models, J. Fluids Struct. 54 (2015) 802-822.

[21] S. Marrone, A. Colagrossi, F. Gamboli, L. González-Gutiérrez, Numerical study on the dissipation mechanisms in sloshing flows induced by violent and high-frequency accelerations. i. theoretical formulation and numerical investigation, Phys. Rev. Fluids 6 (2021) 114801.

[22] J. Michel, D. Durante, A. Colagrossi, S. Marrone, Energy dissipation in violent three-dimensional sloshing flows induced by high-frequency vertical accelerations, Phys. Fluids 34 (2020) 102114.

[23] J. Michel, A. Colagrossi, M. Antuono, S. Marrone, A regularized high-order diffusive smoothed particle hydrodynamics scheme without tensile instability, Phys. Fluids 35 (2023b).

[24] A. Aalilija, C.A. Gandin, E. Hachem, On the analytical and numerical simulation of an oscillating drop in zero-gravity, Comput. & Fluids 197 (2020) 104362.

[25] P. Sun, C. Pilloton, M. Antuono, A. Colagrossi, Inclusion of an acoustic damper term in weakly-compressible SPH models, J. Comput. Phys. 483 (2023) 112056.

[26] H. Wendland, Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree, Adv. Comput. Math. 4 (1995) 389-396.

[27] A. Colagrossi, B. Bouscasse, M. Antuono, S. Marrone, Particle packing algorithm for SPH schemes, Comput. Phys. Comm. 183 (2012) 1641-1683.

[28] A. Parshikov, S. Medin, Smoothed particle hydrodynamics using interparticle contact algorithms, J. Comput. Phys. 180 (2002) 358-382.

[29] P. Randles, L. Libersky, Normalized SPH with stress points, Internat. J. Numer. Methods Engrg. 48 (2000) 1445-1462.

[30] J. Monaghan, R. Gingold, Shock simulation by the particle method SPH, J. Comput. Phys. 52 (1983) 374-389.

[31] J. Morris, P. Fox, Y. Zhu, Modeling low Reynolds Number Incompressible Flows Using SPH, J. Comput. Phys. 136 (1997) 214-226.

[32] S. Shao, E.Y. Lo, Incompressible sph method for simulating newtonian and non-newtonian flows with a free surface, Adv. Water Resour. 26 (2003) 787-800.

[33] M. Antuono, A. Colagrossi, S. Marrone, C. Lugni, Propagation of gravity waves through an SPH scheme with numerical diffusive terms, Comput. Phys. Comm. 182 (2011) 866-877.

[34] M. Antuono, A. Colagrossi, The damping of viscous gravity waves, Wave Motion 50 (2013) 197-209.

[35] P. Sun, A. Colagrossi, S. Marrone, A. Zhang, The δplus-SPH model: Simple procedures for a further improvement of the SPH scheme, Comput. Methods Appl. Mech. Engrg. 315 (2017) 25-49.

[36] M. Antuono, A. Colagrossi, S. Marrone, Numerical diffusive terms in weakly-compressible SPH schemes, Comput. Phys. Comm. 183 (2012) 2570-2580.

[37] P. Sun, A. Colagrossi, S. Marrone, M. Antuono, A. Zhang, A consistent approach to particle shifting in the δ-Plus-SPH model, Comput. Methods Appl. Mech. Engrg. 348 (2019) 912-934.

[38] N. Grenier, D.L. Touzé, A. Colagrossi, M. Antuono, G. Colicchio, Viscous bubbly flows simulation with an interface SPH model, Ocean Eng. 69 (2013) 88-102.

[39] S. Marrone, A. Colagrossi, D. Le Touzé, G. Graziani, Fast free-surface detection and level-set function definition in sph solvers, J. Comput. Phys. 229 (2010) 3652-3663.

[40] G. Zhou, W. Ge, J. Li, A revised surface tension model for macro-scale particle methods, Powder Technol. 183 (2008) 21-26.

[41] M. Zhang, Simulation of surface tension in 2d and 3d with smoothed particle hydrodynamics method, J. Comput. Phys. 229 (2010) 7238-7259.

[42] S. Adami, X. Hu, N. Adams, A new surface-tension formulation for multi-phase SPH using a reproducing divergence approximation, J. Comput. Phys. 229 (2010) 5011-5021.

[43] M. Zhang, S. Zhang, H. Zhang, L. Zheng, Simulation of surface-tension-driven interfacial flow with smoothed particle hydrodynamics method, Comput. & Fluids 59 (2012) 61-71.

[44] M. Huber, S. Reinhardt, D. Weiskopf, B. Eberhardt, Evaluation of surface tension models for SPH-based fluid animations using a benchmark test, in: VRIPHYS, 2015, pp. 41-50.

[45] M. Moghimi, N. Quinlan, A model for surface tension in the meshless finite volume particle method without spurious velocity, Comput. & Fluids 179 (2019) 521-532.

[46] W. Liu, D. Ma, M. Zhang, A. He, N. Liu, P. Wang, A new surface tension formulation in smoothed particle hydrodynamics for free-surface flows, J. Comput. Phys. 439 (2021) 110203.

[47] B. Zheng, L. Sun, P. Yu, A novel interface method for two-dimensional multiphase SPH: Interface detection and surface tension formulation, J. Comput. Phys. 431 (2021) 110119.

[48] F. Zhang, P. Yang, M. Liu, An improved continuum surface tension model in SPH for simulating free-surface flows and heat transfer problems, J. Comput. Phys. 490 (2023) 112322.

[49] D.S. Morikawa, M. Asai, Surface tension simulations with corrected ale-isph and density-based shifting technique, Comput. Part. Mech. (2023) 1-12.

[50] Y. Xu, G. Yang, D. Hu, Comparison of surface tension models for the simulation of two-phase flow in an ISPH-FVM coupling method, Eur. J. Mech. B Fluids 105 (2024) 57-96.

[51] C. Cen, G. Fourtakas, S. Lind, B.D. Rogers, A single-phase gpu-accelerated surface tension model using sph, Comput. Phys. Comm. 295 (2024) 109012.

[52] D.D. Meringolo, A. Colagrossi, S. Marrone, F. Aristodemo, On the filtering of acoustic components in weakly-compressible SPH simulations, J. Fluids Struct. 70 (2017) 1-23.

[53] H. Lamb, Hydrodynamics, in: Cambridge Mathematical Library, The University Press, New York, USA, 1932.