A parameter-free total Lagrangian smooth particle hydrodynamics algorithm applied to problems with free surfaces

Abstract

This paper presents a new Smooth Particle Hydrodynamics computational framework for the solution of inviscid free surface flow problems. The formulation is based on the Total Lagrangian description of a system of first-order conservation laws written in terms of the linear momentum and the Jacobian of the deformation. One of the aims of this paper is to explore the use of Total Lagrangian description in the case of large deformations but without topological changes. In this case, the evaluation of spatial integrals is carried out with respect to the initial undeformed configuration, yielding an extremely efficient formulation where the need for continuous particle neighbouring search is completely circumvented. To guarantee stability from the SPH discretisation point of view, consistently derived Riemann-based numerical dissipation is suitably introduced where global numerical entropy production is demonstrated via a novel technique in terms of the time rate of the Hamiltonian of the system. Since the kernel derivatives presented in this work are fixed in the reference configuration, the non-physical clumping mechanism is completely removed. To fulfil conservation of the global angular momentum, a posteriori (least-squares) projection procedure is introduced. Finally, a wide spectrum of dedicated prototype problems is thoroughly examined. Through these tests, the SPH methodology overcomes by construction a number of persistent numerical drawbacks (e.g. hour-glassing, pressure instability, global conservation and/or completeness issues) commonly found in SPH literature, without resorting to the use of any ad-hoc user-defined artificial stabilisation parameters. Crucially, the overall SPH algorithm yields equal second order of convergence for both velocities and pressure.

1 Introduction

Modelling interfacial flows characterised by complex free surface patterns is regarded as one of the most challenging topics in the field of Computational Fluid Dynamics [1,2,3]. Efforts in developing a Lagrangian meshfree method to accurately capture the complex evolution of the flow interfaces still remains a challenging problem [4, 5].

The Smoothed Particle Hydrodynamics (SPH) method, originally developed by Lucy [6] and Gingold and Monaghan [7] for astrophysical and cosmological applications, has since been explored in a wide range of engineering applications, namely shock hydrodynamics [8], free surface flows [9,10,11,12], multi-phase flows [13,14,15,16,17,18], rheology [19,20,21], magneto-hydrodynamics [22] and solid mechanics [23,24,25,26,27,28,29,30,31]. One attractive feature of SPH is its mesh free nature, thus not requiring the use of an underlying (or background) mesh. The absence of mesh and the calculation of the interactions among particles based exclusively on their separation allow ease of computation for large deformation problems. For its low computational cost, reasonable accuracy and stability, as well as its ability to handle extremely large distortions [32,33,34,35,36,37], the SPH methodology has been shown to be very competitive [4].

However, as it is well established, the traditional SPH method presents a number of numerical drawbacks, namely: (1) non-physical clumping of particles due to tensile instability [38,39,40,41] when using a non-Lagrangian description of the problem, (2) numerical errors near the boundaries due to reduced compact support [42, 43], (3) zero-energy modes due to rank-deficiency inherent to the use of a nodal integration approach [44], (4) pressure instabilities when considering weakly compressible fluids [45, 46] and (5) reduced order of convergence for volumetric strain and pressure [47,48,49,50].

In the last two decades, substantial effort has been devoted to address above deficiencies and thus enhance the robustness of the standard SPH method. Specifically, corrections have been applied in order to ensure reproducibility of complete polynomials in finite domains, as well as to pass the patch test. Without being exhaustive, Johnson et al. [25, 51] proposed a normalised smoothing technique; Chen et al. [52, 53] introduced corrected kernel approximations on the basis of a Taylor series expansion; Bonet et al. [54,55,56] introduced corrections in the kernel functions and in their derivatives; Liu et al. [27] presented an in-depth discussion about reproducibility properties of the SPH method. Other notable modifications to the SPH method include the moving least-square particle hydrodynamics [57, 58] and reproducing kernel particle method [59, 60]. Problems associated with the stability and convergence of meshfree methods [24, 43, 61] were also reported in [61].

However, above enhanced (also known as modified or corrected) SPH methodologies still suffer from one or more of the following spurious mechanisms, namely zero-energy modes, pressure oscillations and the possible development of long-term instabilities. These numerical artefacts can be alleviated through the use of numerical dissipation mechanisms as follows. First, for the treatment of zero-energy modes, one popular option is the incorporation of ad-hoc artificial viscosity stabilisation, as proposed by Monaghan [8, 9, 34], to the linear momentum equation. Even though the viscosity term requires several user defined (problem dependent) stabilisation parameters, this approach is still being widely accepted and implemented in many open source codes [10, 62, 63]. Morris and Monaghan [64] further improved the accuracy of the stabilisation by introducing an appropriate limiter. Second, for the treatment of pressure instabilities, Monaghan [9, 65] introduced a smoothing procedure for the density update, named “X-SPH”, by penalising the difference between the particle velocity and the average velocity in its neighbourhood. Another a posteriori approach [10, 14, 62, 63] is the re-initialisation of the density field via a corrected (or re-normalised) kernel function after a certain number of time steps. Unfortunately, as reported in [11], both of these techniques may potentially lead to excessive dissipation, which is clearly undesirable in practice. In parallel, Molteni and Colagrossi [45] introduced another ad-hoc approach, named $\delta$-SPH method [66, 67], by incorporating a (Laplacian based) diffusive term to the continuity equation. Further improvement has been made by Sun et al. [68] whereby a particle shifting procedure is combined. On another front, several interesting attempts have also been reported at aiming to solve the instability issues via the use of a physically based Riemann solver [69,70,71,72,73,74, 74,75,76,77]. Despite the tremendous development in the field, the ab initio stability of SPH schemes is still an open problem.

In the current paper, a parameter-free Total Lagrangian SPH computational framework is presented with particular attention paid to the analysis of fluid flow problems without significant topological changes on surface (to be explored in a subsequent publication). A first-order system of hyperbolic equations¹ is presented in terms of the linear momentum $\varvec{p}$ and the Jacobian J of the deformation. The mixed $\{\varvec{p},J\}$ system ensures the pressure to converge at the same spatial rate than the velocities and displacements, advantageous in those scenarios (especially when considering inviscid fluid flow problems) where the pressure plays a dominant role. One contribution of this paper is the reformulation of the Updated Lagrangian version of the inviscid fluid equations typically used in SPH in a Total Lagrangian description. In this case, the evaluation of spatial integrals is carried out with respect to the initial undeformed configuration, yielding an extremely efficient formulation where the need for continuous particle neighbouring search is avoided altogether. One of the objectives of this paper is to lay the foundation for seeking a new SPH approach with the final aim to handle in the near future violent free surface flows typically accompanied with large topological changes.

A crucial aspect that requires special attention is that of the stability of the SPH formulation. Specifically, a Riemann based (upwinding) approach is pursued where a consistently derived numerical stabilisation is carefully introduced ensuring global production of numerical entropy. The latter is demonstrated by the monitoring of the so-called Hamiltonian of the system [88]. For completeness, a global a posteriori angular momentum projection procedure is also presented with the aim to preserve the angular momentum of the system. On the shown test problems, the proposed SPH algorithm is shown to effectively eliminate the appearance of spurious hourglass-like modes and pressure instabilities.

The paper is organised as follows. Section 2 summarises the system of first-order conservation laws $\{\varvec{p}, J\}$ for inviscid fluid flow problems used in this paper. Section 3 presents the weak variational statements of the problem as well the second law of thermodynamics written in terms of the so-called Hamiltonian [86, 89, 90]. Section 4 presents the Smooth Particle Hydrodynamics discretisation numerical scheme where special attention is paid to the Riemann based (upwinding) numerical dissipation employed. A proof of entropy production is included as a necessary condition for stability at the semi-discrete level. Section 5 describes the explicit Total Variation Diminishing Runge-Kutta time integrator used for temporal discretisation. Section 6 summarises the flowchart of the proposed SPH methodology. In Sect. 7, an extensive set of challenging numerical examples is examined to assess the performance of the proposed algorithm. Finally, Sect.  8 presents some concluding remarks and current directions of research.

2 Total Lagrangian first-order conservation laws

Fig. 1

Motion of a deformable continuum domain

2.1 Kinematics and conservation of volume map

Consider the three dimensional deformation of a continuum moving from its initial (material) configuration occupying a volume $\Omega _0$, of boundary $\partial \Omega _0$ with outward unit normal $\varvec{N}$, into a current configuration at time t occupying a volume $\Omega (t)$ of boundary $\partial \Omega (t)$ with outward unit normal $\varvec{n}$, see Fig. 1. The standard notation and definitions for the deformation gradient tensor $\varvec{F}$ (used to map differential fibre elements $d\varvec{x}=\varvec{F}d\varvec{X}$) and its determinant J (used to map differential volume elements $d\Omega = J d \Omega _0$) are introduced

$$\begin{aligned} \varvec{F} = \frac{\partial \varvec{x}}{\partial \varvec{X}} = \varvec{\nabla }_0 \varvec{x}; \qquad J = \det \varvec{F}, \end{aligned}$$

(1)

where $\varvec{x}$ represents the current position of a particle originally at $\varvec{X}$ and $\varvec{\nabla }_0$ denotes the gradient with respect to the material configuration. Similarly, material differential area vectors $d\varvec{A}$ (co-linear with $\varvec{N}$) are mapped to current differential area vectors $d \varvec{a}$ (co-linear with $\varvec{n}$) through the co-factor $\varvec{H}$ of the deformation, which is related to the deformation gradient tensor $\varvec{F}$ via the Nanson’s rule [86]

$$\begin{aligned} \varvec{H} = J \varvec{F}^{-T}. \end{aligned}$$

(2)

An alternative way to Eq. (1b) to compute the Jacobian of the deformation is possible via an integral conservation law [29,30,31, 80, 84,85,86] as follows,

$$\begin{aligned} \frac{{\mathrm{d}}}{{\mathrm{d}}t} \int _{\Omega _0} J \, {\mathrm{d}}\Omega _0 = \int _{\partial \Omega _0} \left( \varvec{H}^{\mathrm{T}} \varvec{v} \right) \cdot \varvec{N} {\mathrm{d}}A, \end{aligned}$$

(3)

where $\varvec{v}=\frac{{\mathrm{d}}\varvec{x}}{{\mathrm{d}}t}$ is the velocity field. The equivalent local conservation law and associated jump condition across a discontinuity surface with normal $\varvec{N}$, propagating with speed U in the reference space [78, 83, 84], are

$$\begin{aligned} \frac{\partial J}{\partial t} - \text {DIV} \left( \varvec{H}^{\mathrm{T}} \varvec{v} \right) = 0; \qquad U \llbracket J \rrbracket = - \llbracket \varvec{v} \rrbracket \cdot \left( \varvec{H}^{\text {Ave}} \varvec{N} \right) . \end{aligned}$$

(4)

In this expression, DIV represents the material divergence operator, $\varvec{H}^{\text {Ave}} = \frac{1}{2} \left( \varvec{H}^- + \varvec{H}^+ \right)$ represents an average state of the co-factor of the deformation between the left and right states across the discontinuity surface and $\llbracket \cdot \rrbracket = \left[ \cdot \right] ^+ - \left[ \cdot \right] ^-$ represents the jump operator.

2.2 Conservation of linear momentum

The conservation of linear momentum per unit undeformed volume $\varvec{p} =\rho _0 \varvec{v}$ (with $\rho _0$ the material density of the continuum) [29,30,31, 78,79,80,81,82,83,84,85,86,87] is established for any arbitrary Lagrangian material volume $\Omega _0$ by

$$\begin{aligned} \frac{{\mathrm{d}}}{{\mathrm{d}}t} \int _{\Omega _0} \varvec{p} \, {\mathrm{d}}\Omega _0 = \int _{\partial \Omega _0} \varvec{t} \, {\mathrm{d}}A + \int _{\Omega _0} \varvec{f}_0 \, {\mathrm{d}}\Omega _0, \end{aligned}$$

(5)

where $\varvec{f}_0$ is a body force per unit of undeformed volume and $\varvec{t} = \varvec{P} \varvec{N}$ is the traction vector associated with the material outward unit normal surface vector $\varvec{N}$ with $\varvec{P}$ being the first Piola-Kirchhoff stress tensor. The equivalent local equilibrium equation and the corresponding jump condition across a discontinuity [78, 82,83,84] can be written as

$$\begin{aligned} \frac{\partial \varvec{p}}{\partial t} - \text {DIV} \varvec{P} = \varvec{f}_0; \qquad U \rho _0 \llbracket \varvec{v} \rrbracket = - \llbracket \varvec{P} \rrbracket \varvec{N}. \end{aligned}$$

(6)

2.3 Combined equations

Combining expressions (3) and (5), above global conservation laws can now be summarised in a concise manner as

$$\begin{aligned} \frac{{\mathrm{d}}}{{\mathrm{d}}t} \int _{\Omega _0} \varvec{\mathcal {U}} \, {\mathrm{d}} \Omega _0 = - \int _{\partial \Omega _0} \varvec{\mathcal {F}}_{\varvec{N}} \, {\mathrm{d}}A + \int _{\Omega _0} \varvec{\mathcal {S}} \, {\mathrm{d}} \Omega _0, \end{aligned}$$

(7)

representing a system of Total Lagrangian conservation laws where $\varvec{\mathcal {U}}$ denotes the set of conservation variables, $\varvec{\mathcal {S}}$ the source term and $\varvec{\mathcal {F}}_{\varvec{N}}$ the surface flux vector, described as follows

$$\begin{aligned} \varvec{\mathcal {U}} = \left[ \begin{array}{c} \varvec{p} \\ J \\ \end{array} \right] ; \qquad \varvec{\mathcal {F}}_{\varvec{N}} = - \left[ \begin{array}{c} \varvec{t} \\ \frac{1}{\rho _0} \varvec{p} \cdot \left( \varvec{H} \varvec{N} \right) \\ \end{array} \right] ; \qquad \varvec{\mathcal {S}} = \left[ \begin{array}{c} \varvec{f}_0 \\ 0 \\ \end{array} \right] , \end{aligned}$$

(8)

where $\varvec{H}$ is the cofactor of the deformation previously defined in (2). The computation of $\varvec{H}$ is carried out from the use of J in (4) and the computation of the deformation gradient tensor (1) obtained from the time integration of the velocity field $\varvec{v}$ deduced from (6).

For closure, above system (8) requires the introduction of an appropriate “elastic fluid” type of constitutive law (also known as polytropic equation of state [34, 54]), relating the stress tensor $\varvec{P}$ with the Jacobian of the deformation J. In this specific case, a strain energy $\Psi$ per unit undeformed volume is introduced as

$$\begin{aligned} \Psi (J) = \left\{ \begin{array}{ll} \kappa \left( J - 1 - \frac{J^{1 - \gamma } - 1}{1 - \gamma } \right) , &{} \gamma > 1 \\ \kappa (J - 1 - \ln J), &{} \gamma = 1 \end{array} \right. \end{aligned}$$

(9)

with $\{ \gamma , \kappa \} > 0$ the positive material constants.² Thus, the first Piola-Kirchhoff stress tensor $\varvec{P}$ can be expressed as

$$\begin{aligned} \varvec{P} = - p \varvec{H}; \qquad p = - \frac{{\mathrm{d}} \Psi }{{\mathrm{d}}J} = \kappa \left[ \left( \frac{1}{J} \right) ^{\gamma } - 1 \right] , \end{aligned}$$

(10)

where a positive value of the pressure p in above equation indicates compression. So far only an elastic volumetric contribution (10) is considered. General materials however will also exhibit additional inelastic effects $\varvec{P}^{v}$, where the combined stress becomes

$$\begin{aligned} \varvec{P} = - p \varvec{H} + \varvec{P}^{v}. \end{aligned}$$

(11)

In the case of fluids this is due to viscosity whereas for solids more complex constitutive equations possibly involving plasticity [31] can be required. In the current paper, we consider the case of perfect fluid, that is $\varvec{P}^{v} = \varvec{0}$. Physically, this implies that no shear/deviatoric stresses, viscosity or heat conduction are taken into account.

Importantly, with the aid of (9), hyperbolicity of system (8) can be demonstrated (see “Appendix”) ensuring the existence of real wave speeds at any state of deformation. Finally, for the complete definition of the initial boundary value problem, initial and boundary (essential and natural) conditions must also be specified as appropriate.

Remark 1

An alternative non-conservative form for J (4a) is achieved by introducing the Piola identity, namely $\text {DIV} \varvec{H} = \varvec{0}$, into Eq. (4a) to give

$$\begin{aligned} \frac{\partial J}{\partial t} - \varvec{H} : \varvec{\nabla }_0 \varvec{v} = 0. \end{aligned}$$

(12)

3 Weak form statements

In general, a standard weak variational statement for the system (8) is established by multiplying the local differential equations $\varvec{\mathcal {U}}$ (8) with their appropriate work conjugate virtual fields $\delta \varvec{\mathcal {V}}$, and integrating over the reference domain $\Omega _0$ of the body, to give

$$\begin{aligned}&\int _{\Omega _0} \delta \varvec{\mathcal {V}} \bullet \frac{\partial \varvec{\mathcal {U}}}{\partial t} \, {\mathrm{d}}{\Omega _0} = \int _{\Omega _0} \delta \varvec{\mathcal {V}} \bullet \varvec{\mathcal {S}} \, {\mathrm{d}}{\Omega _0} + \int _{\Omega _0} \varvec{\mathcal {F}}_I \bullet \frac{\partial \delta \varvec{\mathcal {V}}}{\partial X_I} \, {\mathrm{d}}{\Omega _0} \nonumber \\&\qquad \quad - \int _{\partial \Omega _0} \delta \varvec{\mathcal {V}} \bullet \varvec{\mathcal {F}_N} \, {\mathrm{d}}A, \end{aligned}$$

(13)

where the normal fluxes $\varvec{\mathcal {F}}_{\varvec{N}} = \varvec{\mathcal {F}}_I N_I$ with $N_I$ being the material outward normal in I-th material direction and the symbol $\bullet$ is used to denote the inner (dual) product of work conjugate pairs.

In order to give a proper physical meaning to the conjugate virtual fields $\delta \varvec{\mathcal {V}}$ and pave the way for the proof of numerical entropy production presented in a subsequent section, it is useful to introduce the notion of Hamiltonian $\mathcal {H}$ = $\mathcal {H} (\varvec{X}, t)$ [86, 88]. This indeed can be understood as a generalised convex entropy function of the system of conservation laws (8). Specifically, the Hamiltonian $\mathcal {H}$ is defined by

$$\begin{aligned} \mathcal {H} (\varvec{X}, t) = \hat{\mathcal {H}} (\varvec{\mathcal {U}}) = \frac{1}{2 \rho _0} \varvec{p} \cdot \varvec{p} + \Psi (J), \end{aligned}$$

(14)

which represents the summation of the kinetic energy (i.e. the first term on the right hand side of (14)) and Helmholtz free energy $\Psi$ per unit of undeformed volume. Note here that $\mathcal {H} (\varvec{X}, t)$ and $\hat{\mathcal {H}} (\varvec{\mathcal {U}})$ represent alternative functional representations of the same quantity. With expression (14) at hand, the work conjugates $\varvec{\mathcal {V}}$ can then be obtained by taking simple differentiation to give

$$\begin{aligned} \varvec{\mathcal {V}} = \frac{\partial \hat{\mathcal {H}}}{\partial \varvec{\mathcal {U}}} = \left[ \begin{array}{c} \frac{\partial \hat{\mathcal {H}}}{\partial \varvec{p}} \\ \frac{\partial \hat{\mathcal {H}}}{\partial J} \end{array} \right] = \left[ \begin{array}{c} \varvec{v} \\ - p \end{array} \right] . \end{aligned}$$

(15)

Finally, upon substitution of (15) into (13), individual variational statement for the linear momentum $\varvec{p}$ and the Jacobian of the deformation J are

$$\begin{aligned} \begin{aligned} \int _{\Omega _0} \delta \varvec{v} \cdot \frac{\partial \varvec{p}}{\partial t} \, {\mathrm{d}}{\Omega _0}&= - \int _{\Omega _0} \varvec{P} : \varvec{\nabla }_0 \delta \varvec{v} \, {\mathrm{d}}{\Omega _0} \\&\quad +\,\int _{\Omega _0} \delta \varvec{v} \cdot \varvec{f}_0 \, {\mathrm{d}}{\Omega _0} + \int _{\partial \Omega _0} \delta \varvec{v} \cdot \varvec{t}_B \,{\mathrm{d}}A; \\ \int _{\Omega _0} \delta p\frac{\partial J}{\partial t} \, {\mathrm{d}} \Omega _0&= \int _{\Omega _0} \delta p \left( \varvec{H} : \varvec{\nabla }_0 \varvec{v} \right) {\mathrm{d}}{\Omega _0}. \end{aligned} \end{aligned}$$

(16)

Here, both $\delta \varvec{v}$ and $\delta p$ represent the virtual velocity and pressure fields, respectively.

3.1 Second law of thermodynamics

It is instructive to revisit the second law of thermodynamics when written in terms of the Hamiltonian. The time rate of the volume integral of the Hamiltonian is obtained as follows

$$\begin{aligned} \begin{aligned} \frac{{\mathrm{d}}}{{\mathrm{d}}t}\int _{\Omega _0} \mathcal {H} \, {\mathrm{d}}\Omega _0&= \int _{\Omega _0} \varvec{\mathcal V} \bullet \frac{\partial \varvec{\mathcal U}}{\partial t}\, {\mathrm{d}}{\Omega _0}\\&= \int _{\Omega _0} \left( \varvec{v} \cdot \frac{\partial \varvec{p}}{\partial t} - p \frac{\partial J}{\partial t} \right) \, {\mathrm{d}}{\Omega _0} \\&= \int _{\Omega _0} \left( \varvec{v} \cdot \frac{\partial \varvec{p}}{\partial t} - p \varvec{H} : \varvec{\nabla }_0 \varvec{v} \right) \, {\mathrm{d}}{\Omega _0} \\&= \int _{\Omega _0} \left( \varvec{v} \cdot \frac{\partial \varvec{p}}{\partial t} + (\varvec{P} - \varvec{P}^{v}) : \varvec{\nabla }_0 \varvec{v} \right) \, {\mathrm{d}}{\Omega _0} \\ \end{aligned} \end{aligned}$$

(17)

where, in the first line of (17), use has been made of the conjugacy of the fields (15). In addition, Eqs. (12) and (11) have been substituted in the second and third lines of (17), respectively. Subsequently, we can substitute the linear momentum Eqs. (6) into (17) to give

$$\begin{aligned}&\frac{{\mathrm{d}}}{{\mathrm{d}}t}\int _{\Omega _0} \mathcal {H}\, {\mathrm{d }}\Omega _0 =\int _{\Omega _0} \left( \varvec{v} \cdot \varvec{f}_0 + \varvec{v} \cdot \text {DIV} \varvec{P} + \varvec{P} : \varvec{\nabla }_0 \varvec{v} \right) \, {\mathrm{d}}{\Omega _0} \nonumber \\&\qquad \qquad - \int _{\Omega _0} \varvec{P}^{v} : \varvec{\nabla }_0 \varvec{v} \, {\mathrm{d}}{\Omega _0}. \end{aligned}$$

(18)

Recalling that $\varvec{v} \cdot \text {DIV} \varvec{P} + \varvec{P} : \varvec{\nabla }_0 \varvec{v} = \text {DIV} (\varvec{P}^{\mathrm{T}} \varvec{v})$, above equation reduces to

$$\begin{aligned}&\frac{{\mathrm{d}}}{{\mathrm{d}}t}\int _{\Omega _0} \mathcal {H} \, {\mathrm{d}} \Omega _0 = \int _{\Omega _0} \left( \varvec{v} \cdot \varvec{f}_0 + \text {DIV} \left( \varvec{P}^{\mathrm{T}} \varvec{v}\right) \right) \,{\mathrm{d}} \Omega _0 \nonumber \\&\qquad \qquad - \int _{\Omega _0} \varvec{P}^{v} : \varvec{\nabla }_0 \varvec{v} \,d \Omega _0. \end{aligned}$$

(19)

By performing integration by parts of the $\text {DIV}$ term in Eq. (19) and after some re-arrangement, it yields

$$\begin{aligned} \frac{{\mathrm{d}}}{{\mathrm{d}}t}\int _{\Omega _0} \mathcal {H} \, {\mathrm{d}} \Omega _0 - \dot{\Pi }_{\text {ext}} = - \int _{\Omega _0} \varvec{P}^{v} : \varvec{\nabla }_0 \varvec{v} \,{\mathrm{d}} \Omega _0, \end{aligned}$$

(20)

where $\dot{\Pi }_{\text {ext}}$ denotes the power introduced by external forces, defined as

$$\begin{aligned} \dot{\Pi }_{\text {ext}} = \int _{\Omega _0} \varvec{v} \cdot \varvec{f}_0 \,{\mathrm{d}}\Omega _0 + \int _{\partial {\Omega _0}} \varvec{v} \cdot \varvec{t}_B\,{\mathrm{d}}A. \end{aligned}$$

(21)

In expression (20) , the term on the right hand side represents possible physical dissipation and is always non-positive (namely, $-\varvec{P}^v : \varvec{\nabla }_0 \varvec{v} \le 0$) and, consequently, Eq. (20) can be transformed into the following inequality

$$\begin{aligned} \frac{{\mathrm{d}}}{{\mathrm{d}}t}\int _{\Omega _0} \mathcal {H} \, {\mathrm{d}}\Omega _0- \dot{\Pi }_{\text {ext}} \le 0, \end{aligned}$$

(22)

which represents a valid expression for the second law of thermodynamics [91]. Satisfaction of inequality (22) is a necessary ab initio condition to ensure stability, otherwise referred to as the Coleman–Noll procedure [91].

In the specific case of an isolated system (i.e. $\dot{\Pi }_{\text {ext}}=0$), inequality (22) reduces to

$$\begin{aligned} \frac{{\mathrm{d}}}{{\mathrm{d}}t}\int _{\Omega _0} \mathcal {H} \, {\mathrm{d}}\Omega _0 \le 0. \end{aligned}$$

(23)

This implies that for an isolated system, the decrease in the Hamiltonian is intrinsically related to the dissipation introduced by any inelastic (i.e. viscous) effects. This key concept will be further exploited in Sect. 4.2 at a semi-discrete level.

4 Spatial discretisation

For the case of corrected SPH methods [54], both the problem variables $\varvec{\mathcal {U}}$ and the conjugate pairs $\delta \varvec{\mathcal {V}}$ are in general interpolated at any given position via corrected SPH shape functions $\tilde{\mathrm {N}}$ (or kernel functions $\tilde{W}$) with a given compact support of radius 2h around every particle (see Fig. 2). Specifically, for a given position $\varvec{X}_a$, both $\varvec{\mathcal {U}}$ and $\delta \varvec{\mathcal {V}}$ can be approximated as

$$\begin{aligned}&\varvec{\mathcal {U}}_a (t) \approx \sum _{b \in \Lambda ^b_a} \tilde{\mathrm {N}}_b (\varvec{X}_a) \varvec{\mathcal {U}}_b (t); \qquad \delta \varvec{\mathcal {V}}_a \approx \sum _{b \in \Lambda ^b_a} \tilde{\mathrm {N}}_b (\varvec{X}_a) \delta \varvec{\mathcal {V}}_b; \nonumber \\&\quad \tilde{\mathrm {N}}_b (\varvec{X}_a) = \Omega _0^b \tilde{W}_b (\varvec{X}_a). \end{aligned}$$

(24)

Here, $\Lambda ^b_a$ represents the set of neighbouring particles b that lie inside a sphere of a given radius around $\varvec{X}_a$, $\Omega _0^b$ represents the volume associated to particle b and $\varvec{\mathcal {U}}_b (t)$ and $\delta \varvec{\mathcal {V}}_b$ represent the time varying variables and their virtual conjugate pairs stored at particle b, respectively. In addition, $\left[ \bullet \right] _a (t)$ represent the problem variables at position $\varvec{X}_a$ and time t, that is $\left[ \bullet \right] _a (t) = \left[ \bullet \right] (\varvec{X}_a, t)$. Additionally, the use of corrected SPH shape functions $\tilde{\mathrm {N}}$ ensures that both constant and linear functions are perfectly interpolated [54].

Fig. 2

The compact support of kernel approximation

In line with previous work [31], the material gradient $\varvec{\nabla }_0 (\bullet )$ of any arbitrary vector function $\varvec{f}$ can be approximated through the so-called corrected pseudo area vector $\tilde{\varvec{\mathrm {C}}}_{ab}$ (defined on the basis of corrected kernel gradient $\tilde{\varvec{\nabla }}_0$ [54]), described as

$$\begin{aligned}&\varvec{\nabla }_0 \varvec{f} (\varvec{X}_a) \approx \frac{1}{\Omega _0^a} \sum _{b \in \Lambda ^b_a} \frac{1}{2} \left( \varvec{f}_b - \varvec{f}_a \right) \otimes \tilde{\varvec{\mathrm {C}}}_{ab}; \nonumber \\&\tilde{\varvec{\mathrm {C}}}_{ab}:= 2 \Omega _0^a \varvec{\nabla }_0 \tilde{\mathrm {N}}_b (\varvec{X}_a); \nonumber \\&\varvec{\nabla }_0 \tilde{\mathrm {N}}_b (\varvec{X}_a) = \Omega _0^b \tilde{\varvec{\nabla }}_0 W_b (\varvec{X}_a). \end{aligned}$$

(25)

Specifically, the correction technique involves modifying (uncorrected) pseudo area vector $\varvec{\mathrm {C}}_{ab}$ (or kernel gradient $\varvec{\nabla }_0$) by introducing a particle based correction matrix $\varvec{L}_a$ [54]

$$\begin{aligned}&\tilde{\varvec{\mathrm {C}}}_{ab}= \varvec{L}_a \varvec{\mathrm {C}}_{ab}; \qquad \varvec{\mathrm {C}}_{ab}:= 2 \Omega _0^a \varvec{\nabla }_0 \mathrm {N}_b (\varvec{X}_a); \quad \varvec{\nabla }_0 \mathrm {N}_b (\varvec{X}_a) \nonumber \\&\quad = \Omega _0^b \varvec{\nabla }_0 W_b (\varvec{X}_a), \end{aligned}$$

(26)

from which $\varvec{L}_a$ is explicitly evaluated as [54, 55]

$$\begin{aligned} \varvec{L}_a = \left( \sum _{b \in \Lambda ^b_a} \Omega _0^b \varvec{\nabla }_0 W_b (\varvec{X}_a) \otimes (\varvec{X}_b - \varvec{X}_a) \right) ^{-1}. \end{aligned}$$

(27)

This type of kernel gradient correction [28, 54] ensures the gradient of any linear field distribution is exactly evaluated. Moreover, the term $- \varvec{f}_a$ is added in (25) in order to ensure the gradient vanishes for a constant (uniform) field [8]. Noticing that if no correction is applied (i.e. $\varvec{L}_a$ = $\varvec{I}$), the corrected pseudo area vector $\tilde{\varvec{\mathrm {C}}}_{ab}$ (26) degenerates to the usual SPH uncorrected area vector $\varvec{\mathrm {C}}_{ab}$ [8].

It is worth pointing out that the correction matrices (27) applied to a pair of interacting particles a and b are in general not the same, that is $\varvec{L}_a \ne \varvec{L}_b$. This would consequently destroy the anti-symmetric condition of the gradient evaluations (i.e. $\tilde{\varvec{\mathrm {C}}}_{ab}\ne - \tilde{\varvec{\mathrm {C}}}_{ba}$) [92]. One possible option to enforce such condition consists of skew-symmetrising the area vector $\tilde{\varvec{\mathrm {C}}}_{ab}^{\text {Skew}}$ by taking the difference between the area vectors of particle pairs to give

$$\begin{aligned} \begin{aligned} \tilde{\varvec{\mathrm {C}}}_{ab}^{\text {Skew}}&= \frac{1}{2} \left( \tilde{\varvec{\mathrm {C}}}_{ab}- \tilde{\varvec{\mathrm {C}}}_{ba}\right) \\&=\varvec{L}_{ab}^{\text {Ave}} \varvec{\mathrm {C}}_{ab}; \qquad \qquad \qquad \qquad \varvec{L}_{ab}^{\text {Ave}}=\frac{1}{2}(\varvec{L}_a + \varvec{L}_b). \end{aligned} \end{aligned}$$

(28)

4.1 Total Lagrangian SPH discrete formulation

Typically, in SPH type approaches [8, 14, 29,30,31, 54], the above weak statements (16) are under-integrated by using the cloud of particles as quadrature point [8, 28, 41, 54, 55, 93]. Doing this may potentially excite the instabilities (e.g. zero-energy modes and pressure instability) arising from rank-deficiency inherent to the use of particle (or nodal) integration approach. For instance, the error accumulation over time can eventually lead to the breakdown of a scheme. Additionally, in order to ensure the completeness conditions (low order polynomials are interpolated exactly) of a SPH algorithm, the type of correction technique described in expressions (24) and (25) are typically used when evaluating the integrals of (16). This correction method however is well known to destroy the desired anti-symmetric property of the pseudo area vector. In other words, pseudo area vectors between all pairwise interacting particles are no longer equal and opposite.

Now the question remains on how to construct a conservative SPH method capable of addressing persistent SPH numerical deficiencies (e.g. hour-glassing and pressure modes) described above, whilst still minimising the errors caused by the dissatisfaction of the anti-symmetric property as the result of gradient correction. To achieve this, and taking inspiration from previous work [31], a SPH method in conjunction with well established stabilisation procedure is presented. The SPH discretisation for the system $\{\varvec{p}, J\}$ described in (7) becomes

$$\begin{aligned} \Omega _0^a \frac{{\mathrm{d}} \varvec{p}_a}{{\mathrm{d}}t}&= \sum _{b \in \Lambda ^b_a} \varvec{T}_{ab} + A_0^a \varvec{t}_B^a + \Omega _0^a \varvec{f}^a_0 + \sum _{b \in \Lambda ^b_a} \varvec{\mathcal {D}}_{\varvec{v}}^{ab}+ \sum _{b \in \Lambda ^b_a} \varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab} ; \end{aligned}$$

(29a)

$$\begin{aligned} \Omega _0^a \frac{{\mathrm{d}} J_a}{{\mathrm{d}}t}&= \varvec{H}_{a} : \left[ \sum _{b \in \Lambda ^b_a} \frac{1}{2} \left( \varvec{v}_b - \varvec{v}_a \right) \otimes \tilde{\varvec{\mathrm {C}}}_{ab} \right] + \sum _{b \in \Lambda ^b_a} \mathcal {D}_p^{ab}. \end{aligned}$$

(29b)

In above expressions, the pairwise interacting force

$$\begin{aligned} \varvec{T}_{ab} = \frac{1}{2} \left[ \varvec{P}_a \tilde{\varvec{\mathrm {C}}}_{ab} - \varvec{P}_b \tilde{\varvec{\mathrm {C}}}_{ba} \right] , \end{aligned}$$

(30)

clearly satisfies local conservation by construction, namely $\varvec{T}_{ab} = - \varvec{T}_{ba}$. Moreover, the discrete cofactor $\varvec{H}_{a}=J_a \varvec{F}_a^{-T}$ (2) is evaluated using the discrete deformation gradient $\varvec{F}_{a}$ accurately approximated as

$$\begin{aligned} \varvec{F}_a = \frac{1}{\Omega _0^a} \sum _{b \in \Lambda ^b_a} \frac{1}{2} \left( \varvec{x}_b - \varvec{x}_a \right) \otimes \tilde{\varvec{\mathrm {C}}}_{ab}, \end{aligned}$$

(31)

where $\varvec{x}_a$ represent the current nodal geometry at particle a obtained from time integration of the velocity $\varvec{v}_a$.

Finally, the remaining terms to be defined in Eq. (29a, b) are the so-called pairwise stabilisation terms $\{ \varvec{\mathcal {D}}_{\varvec{v}}^{ab}, \mathcal {D}_p^{ab}, \varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab} \}$. In particular, the term $\varvec{\mathcal {D}}_{ab}^{\varvec{v}}$ in (29a) addresses spurious zero-energy (hourglass-like [33]) modes due to rank-deficiency, whereas the term $\mathcal {D}_{ab}^p$ (29b) removes pressure instabilities in near incompressibility regimes [80]. Another challenging aspect when designing a SPH numerical scheme is the ability to control errors arising from the violation of the skew-symmetric nature due to the gradient correction. One viable option is the introduction of an extra stabilisation term $\varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab}$ into expression (29a), with the objective to penalise the difference between the usual uncorrected pseudo area vector $\varvec{\mathrm {C}}_{ab}$ and $\tilde{\varvec{\mathrm {C}}}_{ab}^{\text {Skew}}$ (28). Mathematically, this is described as $\varvec{\mathrm {C}}_{ab}- \tilde{\varvec{\mathrm {C}}}_{ab}^{\text {Skew}}$.

To ensure the global conservation of the SPH discretised Eq. (29), and following the style of (30), we choose to strongly enforce local conservation for the stabilisations involved, in such a way that $\varvec{\mathcal {D}}_{\varvec{v}}^{ab} = - \varvec{\mathcal {D}}_{\varvec{v}}^{ba}$, $\mathcal {D}_p^{ab} = - \mathcal {D}_p^{ba}$ and $\varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab} = - \varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ba}$. Definition of these terms follows via the semi-discrete version of the Colemann-Noll procedure [94] in order to ensure the production of numerical entropy. This will be discussed in the following section.

4.2 Numerical entropy production

In this section, we present a procedure to obtain consistent and locally conservative stabilisation terms by utilising the concept of the rate of entropy production [95,96,97], which can be understood as a semi-discrete version of the classical Coleman–Noll procedure. Specifically, the semi-discrete counterpart of (17) is

$$\begin{aligned} \begin{aligned} \sum _a \Omega _0^a \frac{{\mathrm{d}} \mathcal {H}_a}{{\mathrm{d}} t}&= \sum _a \Omega _0^a \varvec{\mathcal {V}}_a \cdot \frac{{\mathrm{d}} \varvec{\mathcal {U}}_a}{{\mathrm{d}}t} \\&= \sum _a \Omega _0^a \left[ \varvec{v}_a \cdot \frac{{\mathrm{d}} \varvec{p}_a}{{\mathrm{d}} t} - p_a \frac{{\mathrm{d}} J_a}{{\mathrm{d}} t} \right] . \end{aligned} \end{aligned}$$

(32)

Subsequently, we can substitute the linear momentum equation (29a) and the Jacobian equation (29b) into (32), and after some algebra, gives

$$\begin{aligned} \begin{aligned} \sum _a \Omega _0^a \frac{{\mathrm{d}} \mathcal {H}_a}{{\mathrm{d}}t}&= \sum _a \sum _{b \in \Lambda ^b_a} \frac{1}{2} \left[ \left( \varvec{P}_a \tilde{\varvec{\mathrm {C}}}_{ab}\right) \cdot \varvec{v}_b\right. \\&\qquad -\,\left. \left( \varvec{P}_b \tilde{\varvec{\mathrm {C}}}_{ba}\right) \cdot \varvec{v}_a \right] + \dot{\Pi }_{\text {ext}} \\&\qquad +\,\underbrace{\sum _a \sum _{b \in \Lambda ^b_a} \left( \varvec{v}_a \cdot \varvec{\mathcal {D}}_{\varvec{v}}^{ab} - p_a \mathcal {D}_p^{ab} + \varvec{v}_a \cdot \varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab} \right) }_{- \mathcal {D}_{\text {total}}}. \end{aligned} \end{aligned}$$

(33)

Here, $\dot{\Pi }_{\text {ext}}$ denotes the semi-discrete power introduced by external forces expressed as

$$\begin{aligned} \dot{\Pi }_{\text {ext}} = \sum _a \left( \Omega _0^a \varvec{v}_a \cdot \varvec{f}_0^a + A_0^a \varvec{v}_B^a \cdot \varvec{t}_B^a \right) . \end{aligned}$$

(34)

Additionally, due to the anti-symmetric nature of the first term in the first line of (33), we can conclude that this term cancel and thus (33) reduces to

$$\begin{aligned} \begin{aligned} \sum _a \Omega _0^a \frac{{\mathrm{d}} \mathcal {H}_a}{{\mathrm{d}}t} - \dot{\Pi }_{\text {ext}}&= - \mathcal {D}_{\text {total}}. \end{aligned} \end{aligned}$$

(35)

It is now the objective to demonstrate that the term on the right hand side of (35) is always non-positive, that is $\mathcal {D}_{\text {total}} \ge 0$ (to be in agreement with (22)). Specifically,

$$\begin{aligned} \begin{aligned} \mathcal {D}_{\text {total}}&= - \sum _a \sum _{b \in \Lambda ^b_a} \left( \varvec{v}_a \cdot \varvec{\mathcal {D}}_{\varvec{v}}^{ab} - p_a \mathcal {D}_p^{ab} + \varvec{v}_a \cdot \varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab} \right) \\&= - \sum _a \sum _{b \in \Lambda ^b_a} \left( \varvec{v}_b \cdot \varvec{\mathcal {D}}_{\varvec{v}}^{ba} - p_b \mathcal {D}_p^{ba} + \varvec{v}_b \cdot \varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ba} \right) \end{aligned} \end{aligned}$$

(36)

Adding the first line and the second line of the equation above, and again noting the anti-symmetric nature of the stabilisation terms (e.g. $\varvec{\mathcal {D}}_{\varvec{v}}^{ab} = - \varvec{\mathcal {D}}_{\varvec{v}}^{ba}$, $\mathcal {D}_p^{ab} = - \mathcal {D}_p^{ba}$, $\varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab} = - \varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ba}$), an alternative expression for $\mathcal {D}_{\text {total}}$ is

$$\begin{aligned} \mathcal {D}_{\text {total}} = \frac{1}{2} \sum _a \sum _{b \in \Lambda ^b_a} \mathcal {D}_{\text {total}}^{ab} \end{aligned}$$

(37)

with

$$\begin{aligned} \mathcal {D}_{\text {total}}^{ab} = \left( \varvec{v}_b - \varvec{v}_a \right) \cdot \varvec{\mathcal {D}}_{\varvec{v}}^{ab} - ( p_b - p_a ) \mathcal {D}_p^{ab} + \left( \varvec{v}_b - \varvec{v}_a \right) \cdot \varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab} . \end{aligned}$$

(38)

Dissipation terms $\{\varvec{\mathcal {D}}_{\varvec{v}}^{ab},\mathcal {D}_{p}^{ab}, \varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab} \}$ remain to be defined in order to ensure non-negative entropy production, namely, $\mathcal {D}_{\text {total}}>0$. In this work, the term $\varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab}$ is defined as

$$\begin{aligned} \varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab} = \frac{1}{2} \left( \varvec{P}_a + \varvec{P}_b \right) ( \varvec{\mathrm {C}}_{ab}- \tilde{\varvec{\mathrm {C}}}_{ab}^{\text {Skew}} ), \end{aligned}$$

(39)

which is designed in order to penalise the violation of the skew-symmetric nature of the gradient kernel correction. Note that the above definition of $\varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab}$ ensures the skew symmetric nature of this term. As for the remaining stabilisation terms, we propose the following definitions

$$\begin{aligned} \begin{aligned} \varvec{\mathcal {D}}_{\varvec{v}}^{ab}&= \varvec{S}_{ab}^{\varvec{v}} \left( \varvec{v}_b - \varvec{v}_a \right) ; \\ \mathcal {D}_p^{ab}&= - S_{ab}^p (p_b - p_a) , \end{aligned} \end{aligned}$$

(40)

where $\varvec{S}_{ab}^{\varvec{v}}=\varvec{S}_{ba}^{\varvec{v}}$ is a positive definite stabilisation matrix and $S_{ab}^p=S_{ba}^p > 0$. Notice how the dissipation terms $\{\varvec{\mathcal {D}}^{ab}_{\varvec{v}}, \mathcal {D}^{ab}_{p} \}$ defined in (40) are skew symmetric and proportional to the difference (or jump) in velocity $(\varvec{v}_b-\varvec{v}_a)$ and pressure $(p_b - p_a)$ between pairwise particles, typical of Riemann solver based upwinding terms [78]. Moreover, with these definitions for $\{\varvec{\mathcal {D}}^{ab}_{\varvec{v}}, \mathcal {D}^{ab}_{p} \}$ (40), the first two terms on the right hand side of (38) are always strictly positive, thus counterbalancing possible negative values taken by the third term. In this paper, the dissipation terms used are specifically chosen as

$$\begin{aligned} \varvec{S}_{ab}^{\varvec{v}} = \frac{1}{2}\rho _{R, ab}^{\text {Ave}} \Vert \tilde{\varvec{\mathrm {C}}}_{ab}^{\text {Skew}} \Vert c_{p, ab}^{\text {Ave}} \varvec{I};\qquad S_{ab}^p = \frac{\tilde{\varvec{\mathrm {c}}}_{ab}^{\text {Skew}} \cdot \tilde{\varvec{\mathrm {c}}}_{ab}^{\text {Skew}}}{2 \rho _{R, ab}^{\text {Ave}} c_{p,ab}^{\text {Ave}} \Vert \tilde{\varvec{\mathrm {C}}}_{ab}^{\text {Skew}} \Vert }, \end{aligned}$$

(41)

with

$$\begin{aligned} \rho ^{\text {Ave}}_{R,ab} = \frac{1}{2} \left( \rho _{R,a} + \rho _{R,b} \right) ;\qquad c^{\text {Ave}}_{p,ab} = \frac{1}{2} \left( c_{p,a} + c_{p,b} \right) , \end{aligned}$$

(42)

respectively. Here, $\rho _{R,a}$ denotes the material density evaluated at particle a and $c_{p,a}$ represents the pressure wave speed (58) computed at particle a and $\tilde{\varvec{\mathrm {c}}}_{ab}^{\text {Skew}} = \frac{1}{2} \left( \varvec{H}_a \tilde{\varvec{\mathrm {C}}}_{ab} - \varvec{H}_b \tilde{\varvec{\mathrm {C}}}_{ba} \right)$.

Remark 2

Notice that it is always possible to ensure ab initio non-negative entropy production by selecting an appropriate numerical value for the wave speed, possibly different to $c_{p,ab}^{\text {Ave}}$ in (42). For instance, selecting $S^p_{ab}=0$, the pairwise numerical dissipation becomes

$$\begin{aligned} \begin{aligned} \mathcal {D}_{\text {total}}^{ab}&= \left( \varvec{v}_b - \varvec{v}_a \right) \cdot \varvec{\mathcal {D}}_{\varvec{v}}^{ab} + \left( \varvec{v}_b - \varvec{v}_a \right) \cdot \varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab} \\&=\frac{1}{2}\rho ^{\text {Ave}}_{R,ab} \;\tilde{c}_{p, ab} \; \Vert \tilde{\varvec{\mathrm {C}}}_{ab}^{\text {Skew}} \Vert \; \Vert \varvec{v}_b-\varvec{v}_a \Vert ^2 \\&+ (\varvec{v}_b - \varvec{v}_a) \cdot \varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab} \ge 0, \end{aligned} \end{aligned}$$

(43)

where $\tilde{c}_{p, ab}$ is an appropriate wave speed which can be selected as

$$\begin{aligned} \tilde{c}_{p, ab} =\text {max} \left\{ c_{p, ab}^{\text {Ave}}, \frac{-2 (\varvec{v}_b - \varvec{v}_a) \cdot \varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab}}{\rho ^{\text {Ave}}_{R,ab} \Vert \tilde{\varvec{\mathrm {C}}}_{ab}^{\text {Skew}} \Vert \Vert \varvec{v}_b-\varvec{v}_a \Vert ^{2} }\right\} , \end{aligned}$$

(44)

thus ensuring positive entropy production.

5 Time integration

With respect to the time integration of the governing equations, and having in mind a fast and efficient algorithm, we advocate for an explicit type of time integrator. For simplicity, an explicit one-step two-stage Total Variation Diminishing Runge-Kutta (TVD-RK) scheme has been preferred [29,30,31, 78]. This is described by the following time update equations from time step $t^n$ to $t^{n+1}$ as

$$\begin{aligned} \begin{aligned}&\varvec{\mathcal {U}}_a^{\star } = \varvec{\mathcal {U}}_a^n + \Delta t\ \dot{\varvec{\mathcal {U}}}_a^n (\varvec{\mathcal {U}}_a^n, t^n); \\&\varvec{\mathcal {U}}_a^{\star \star } = \varvec{\mathcal {U}}_a^{\star } + \Delta t\ \dot{\varvec{\mathcal {U}}}_a^{\star } (\varvec{\mathcal {U}}_a^{\star }, t^{n+1}); \\&\varvec{\mathcal {U}}_a^{n+1} = \frac{1}{2}(\varvec{\mathcal {U}}_a^n + \varvec{\mathcal {U}}_a^{\star \star }). \end{aligned} \end{aligned}$$

(45)

In this manuscript, the geometry (spatial locations of the particles) is also updated through the above TVD-RK algorithm [29,30,31, 83, 84]. This results in a monolithic time integration procedure where the vector of particle unknowns $\{\varvec{p}_a, J_a, \varvec{x}_a\}$ ($a = 1, \ldots , N$, where N represents the total number of particles across the computational domain) are all updated via (45). The maximum time step $\Delta t$ is governed by the standard Courant–Friedrichs–Lewy (CFL) condition given as

$$\begin{aligned} \Delta t = \alpha _{{\mathrm{CFL}}}\frac{h_{{\mathrm{min}}}}{c_{p,{\mathrm{max}}}}, \end{aligned}$$

(46)

where $c_{p,{\mathrm{max}}}$ is the maximum p-wave speed, $h_{{\mathrm{min}}}$ is the minimum particle spacing within the computational domain and $\alpha _{{\mathrm{CFL}}}$ is the CFL stability number. For the numerical computations presented in this paper, a value of $\alpha _{{\mathrm{CFL}}} = 0.3$, unless otherwise stated, has been chosen to ensure both accuracy and stability [29, 31] of the algorithm.

Given the fact that the interacting force $\varvec{T}_{ab}$ (30) is not co-linear with the position vector $\left( \varvec{x}_b - \varvec{x}_a \right)$ [54], the SPH algorithm described above (29a,b) does not intrinsically fulfil conservation of angular momentum. To address this issue, we adapt the global projection algorithm introduced by Lee et al. [31]. Specifically, the internal nodal force $\sum _{b \in \Lambda ^b_a} \varvec{T}_{ab}$ and the particle based stabilisation terms, namely $\sum _{b \in \Lambda ^b_a} \varvec{\mathcal {D}}^{ab}_{\varvec{v}}$ and $\sum _{b \in \Lambda ^b_a} \varvec{\mathcal {D}}^{ab}_{\Delta \varvec{\mathrm {C}}}$, described in (29a) are suitably modified in a least-square sense in order to preserve the total angular momentum, whilst still ensuring the global conservation of linear momentum.

6 Algorithmic description

For ease of understanding, Algorithm 1 summarises the complete algorithmic description of the mixed-based $\{\varvec{p}, J \}$ Upwind Smooth Particle Hydrodynamics (Upwind-SPH) methodology, with all the necessary numerical ingredients.

7 Numerical examples

In this section, a series of two dimensional (2D) and three dimensional (3D) numerical examples are presented in order to assess the robustness, effectiveness and applicability of the new SPH computational framework. Notice that as the numerical examples included in this paper are restricted to inviscid flow dynamics (and importantly, without severe changes in topology), strictly speaking, no physical dissipation, but only numerical, will be present in our simulations.

For post-processing purposes, one popular option is to plot the solutions directly based on the particle values. This however may not be viable when a fluid patch experiences extremely large deformation and the distance between particles can become excessively large, unless a prohibitively large number of particles is used in the simulation, which is not the most preferred option in industry. In order to circumvent this shortcoming, we present another simple option to visualise the results. This is achieved by introducing a secondary set of particles, with their values approximated using an appropriate kernel interpolation procedure [54] (ensuring zeroth- and first-order completeness). Advantages of this alternative visualisation procedure are demonstrated in the subsequent examples.

7.1 Wave-like cubical fluid patch

As it is well known, one of the persistent drawbacks of the SPH method is the difficulty to ensure a proper order of convergence of the overall algorithm. Spatial order of convergence is still regarded as one of the key challenges by SPH developers [4].

The aim of this example is to examine the spatial convergence of the proposed $\{\varvec{p}, J\}$ Upwind-SPH methodology for the unknown fields velocity and pressure. A three dimensional cube of length $L=2$ m, as shown in Fig. 3a, is analysed with symmetric boundary conditions for all the surfaces. When small deformations are considered, this example has a closed-form solution for the velocity field described as

$$\begin{aligned}&\varvec{v}(\varvec{X},t) = \omega A_0 \cos \left( \omega t\right) \varvec{V}; \nonumber \\&\varvec{V} = \begin{bmatrix} A\sin \left( \frac{\pi X}{L}\right) \cos \left( \frac{\pi Y}{L}\right) \cos \left( \frac{\pi Z}{L}\right) \\ B\cos \left( \frac{\pi X}{L}\right) \sin \left( \frac{\pi Y}{L}\right) \cos \left( \frac{\pi Z}{L}\right) \\ C\cos \left( \frac{\pi X}{L}\right) \cos \left( \frac{\pi Y}{L}\right) \sin \left( \frac{\pi Z}{L}\right) \end{bmatrix}, \end{aligned}$$

(47)

where $\omega =2\pi$ s$^{-1}$ and $A_0=5\times 10^{-4}$ m. The parameters $\{A,B,C\}$ are set as $A=-2$, $B=C=1$ to ensure the existence of an incompressible flow. The problem is initialised with a zero pressure field (see Fig. 3b) together with an initial velocity field defined by using $t = 0$ in (47) (see Fig. 3c). For the fluid patch to be in equilibrium, a specific representation of the body force term is required

$$\begin{aligned} \varvec{f}(\varvec{X},t) =- \omega ^2 A_0 \sin \left( \omega t\right) \varvec{V}. \end{aligned}$$

(48)

An elastic fluid type of constitutive law is used with density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 10$ MPa and $\gamma =1$. As presented in 58, the initial pressure wave speed used in this case is $c_p^0 = \sqrt{\frac{\kappa \gamma }{\rho _0}} = 100$ m/s.

Fig. 3

Wave-like cubical fluid patch: problem setup. (Color figure online)

Convergence analysis is carried out by measuring both the $L_1$ and $L_2$ error norms between the numerical and analytical solutions. Fig. 4a, b show equal second-order convergence for both the velocity field and the pressure field, which is one of the main advantages of the proposed SPH methodology. In this case, mesh refinement improves at the same rate results for both velocity and pressure.

Fig. 4

Wave-like cubical fluid patch: a $L_1$ and (2) $L_2$ norm convergence of components of pressure, X-velocity, Y-velocity and Z-velocity using the proposed {${\varvec{p}, J}$} Upwind-SPH scheme. Results are obtained with $A=-2$, $B=C=1$, $\omega =2\pi$ s$^{-1}$ and $A_0=5\times 10^{-4}$ m at time $t=0.9375$ ms. An elastic fluid is used with density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 10$ MPa, material coefficient $\gamma =1$ and $\Delta t=9.375 \times 10^{-5}$ s

7.2 Stretching of a circular fluid patch

This two dimensional example was first proposed by Monaghan [9], and later explored by several authors in References [54, 92, 98, 99]. The main objective of this example is to verify the accuracy of the developed SPH algorithm and to show the preservation of total linear momentum. In this example, analytical expressions for the velocity (or displacement) and pressure fields are provided. Their detailed derivations can be found in [9, 54, 99].

The circular fluid patch has a radius R of 1 m and is free from any external forces, as illustrated in Fig. 5a. Both the initial velocity field $\varvec{v}_0$ and the initial pressure field $p_0$ are, respectively, described as follows

$$\begin{aligned} \varvec{v}_0=\left[ \begin{array}{c} -A_0X \\ A_0Y \\ 0 \end{array} \right] ; \qquad p_0 = \frac{1}{2}\rho _0 A_0^2 \left[ R^2 - \left( X^2 + Y^2\right) \right] , \end{aligned}$$

(49)

where $A_0 = 100$ s$^{-1}$. For visualisation purposes, contour plots of the above initial conditions are also displayed in Fig. 5b, c. The fluid properties used in this example are density $\rho _0=1000$ kg/m$^3$ and bulk modulus $\kappa = 1.96$ GPa. Utilising the value of $\gamma =1$ in the equation of state, the initial pressure wave speed of the fluid patch can be computed as $c_p^0 := \sqrt{\frac{\kappa \gamma }{\rho _0}} = 1400$ m/s. The domain is spatially discretised with a uniformly distributed particle arrangement. Moreover, the outermost radial particles are set free in order to allow the boundary surface to deform freely.

Fig. 5

Stretching of a circular fluid patch: problem setup

We start off by benchmarking the positions of the particles, situated at the highest and lowest points along the semi-major axis, against the published results reported in References [9, 54]. Due to the symmetric nature of the problem, the solutions of those two particles are found to be identical, and hence only the results based on the highest particle are presented (see Table 1). It can be seen that the results obtained with the proposed method agree extremely well with those published. The maximum difference with respect to the analytical solution is calculated to be 0.75%.

Table 1 Stretching of a circular fluid patch: comparison of position for the highest particle along the semi-major axis from analytical, standard SPH, corrected SPH (CSPH) and $\{ \varvec{p}, J \}$ Upwind-SPH

Table 2 Stretching of a circular fluid patch: comparison of position for farthest particle along the ($\varvec{X}=\left[ 1, 0, 0\right] ^{\mathrm{T}}$ m) and semi-major ($\varvec{X}=\left[ 0, 1, 0\right] ^{\mathrm{T}}$ m) axes between sequential particle density using $\{ \varvec{p}, J \}$ Upwind-SPH against analytical solution at $tA_0=1.294$

Fig. 6

Stretching of a circular fluid patch: Particle refinement of deformed fluid patch alongside with pressure distribution plotted as particles (top) and contour (bottom) at $tA_0=1.294$ using of a 317, b 1257, c 5025, and d 20,081 particles. Red dashed lines represent the analytical deformation. Results obtained using $\{ \varvec{p}, J \}$ Upwind-SPH with an initial velocity and pressure (refer to Eq. 49). An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 1.96$ GPa, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$.

To verify the convergence of the problem, a particle refinement study is examined (see Fig. 6). Particles are sequentially refined in order to compare the resolution of the deformed fluid patch at $tA_0=1.294$. Notice that the particles remain stable (without any oscillations) whilst being stretched. The free surface is accurately captured and no particle clumping is observed. The clumping of SPH particles has been reported in Reference [99] when the particles are distributed along the axes of a Cartesian coordinate system. This shortcoming can also be circumvented by introducing either a boundary conforming particle distribution [54] or a priori particle redistribution algorithm [99].

Fig. 7

Time evolution of a numerical dissipation using a discretisation of 317, 709, 1257 and 2821 particles; and b components of linear momentum using a discretisation of 20,081 particles. Results obtained using $\{ \varvec{p}, J \}$ Upwind-SPH with an initial velocity and pressure (refer to Eq. 49). An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 1.96$ GPa, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$

It can be seen that the pressure plot based on particle values is no longer visible (refer to the top row of Fig. 6) when the fluid patch experiences a very large stretch. For this reason, the visualisation is carried out with the help of a secondary set of particles with values interpolated via SPH shape functions (24c). As shown in the bottom row of Fig. 6, the deformed shapes of the fluid patch are practically identical. Pressure resolution is further enhanced when using a larger number of particles. For qualitative comparison purposes, Table 2 shows the positions of the semi-minor ($\varvec{X}=\left[ 1, 0, 0\right] ^{\mathrm{T}}$ m) and semi-major ($\varvec{X}=\left[ 0, 1, 0\right] ^{\mathrm{T}}$ m) axes at various time instants. As expected, the obtained results converge to the available analytical solution.

Figure 7a shows the time evolution of the accumulated numerical entropy (dissipation) present in the algorithm. This is achieved by integrating the Hamiltonian energy of the system described in (35) over time. The total numerical dissipation is reduced when successively increasing the particle density. Furthermore, Fig. 7b depicts the time history of the components of global linear momentum of the system. As expected, the global linear momentum fluctuates around zero machine accuracy in the absence of external forces.

The predicted solution is expected to oscillate around and then gradually converge towards the incompressible flow solution. The incompressibility limit is governed by the bulk modulus of the fluid, which in turn affects the corresponding pressure wave speed. Notice here that an increase in the bulk modulus of the fluid patch would lead to a smaller time increment used in the simulation. To assess the effectiveness of the algorithm in the near incompressibility regime, the pressure evolution is monitored at position $\varvec{X}=[0,0,0]^{\mathrm{T}}$ using three different values of the bulk modulus, namely $\kappa$, 4$\kappa$ and 16$\kappa$. As shown in Fig. 8, all the results converge towards the incompressible solution without the appearance of non-physical pressure instabilities.

Fig. 8

Stretching of a circular fluid patch: time evolution of non-dimensional pressure at $X=[0,0,0]^{\mathrm{T}}$ m with varying bulk modulus $\kappa$. Results obtained using 1257 particles with an initial velocity and pressure (refer to Eq. 49). An elastic fluid of density $\rho _0=1000$ kg/m$^3$, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$

For completeness, Fig. 9 shows a sequence of snapshots and their pressure contours. Majority of the authors [9, 54, 92, 98, 99] presented the deformation up to $tA_0 =2$. To demonstrate the stability of the algorithm in a long-term response, we continue running the simulation up to $tA_0 =4$. Two particles, one located at the semi-major axis and the other one located at the semi-minor axis, are monitored and compared against the analytical solutions (see Table 3). The particles at both axes are stretched to approximately 6.5 times and 2.9 times of their original positions.

Fig. 9

Stretching of a circular fluid patch: Evolution of the deformation alongside with pressure distribution using $\{ \varvec{p}, J \}$ Upwind-SPH at (left to right) $tA_0= \{ 0.0$, 0.228, 0.421, 0.592, 0.747, 1.019, 1.141, 1.367, 1.472, 1.759, 2.015, 2.467, 3.042, 3.536, $4.0 \}$. Results obtained using 1257 particles with an initial velocity and pressure (refer to Eq. 49). An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 1.96$ GPa, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$.

Table 3 Stretching of a circular fluid patch: Time evolution of farthest particle along the semi-major and semi-minor axes using $\{ \varvec{p}, J \}$ Upwind-SPH

Fig. 10

Rotation of a square fluid patch: problem setup

Fig. 11

Rotation of a square fluid patch: Particle refinement of deformed fluid patch alongside with pressure distribution plotted at particles (top) and as a contour plot (bottom) at $tA_0=1.1$ using a 441, b 1681, c 6561 and d 25,921 particles. Grey dashed line represents the initial position, black dashed lines are the trajectories and red dashed line represents the free surface deformation from ANSYS Fluent using an ultrafine reference mesh. Results obtained using $\{ \varvec{p}, J \}$ Upwind-SPH with an initial velocity and pressure distributions (Refer to Eq. 50). An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 1.96$ GPa, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$

Fig. 12

Rotation of a square fluid patch: time evolution of a trajectory of the corner particle at the material point $\varvec{X}=\left[ -0.5, 0.5, 0\right] ^{\mathrm{T}}$ m; and b numerical dissipation using 441, 1681, 6561 and 25,921 particles. Results obtained using $\{ \varvec{p}, J \}$ Upwind-SPH with an initial velocity and pressure (refer to Eq. 50). An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 1.96$ GPa, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$

Fig. 13

Rotation of a square fluid patch: time evolution of components a global linear momentum; and b global angular momentum using a discretisation of 6561 particles. Results obtained using $\{ \varvec{p}, J \}$ Upwind-SPH with an initial velocity and pressure (refer to Eq. 50). An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 1.96$ GPa, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$

Fig. 14

Rotation of a square fluid patch: deformed fluid patch at $tA_0=1.1$ using a $\{\varvec{p}, J\}$ Upwind-SPH with $\varvec{F}$ ensures zero-th order completeness and b $\{\varvec{p}, J\}$ Upwind-SPH with $\varvec{F}$ ensures zero-th and first-order completeness. Results obtained using 6561 particles with an initial velocity and pressure (refer to Eq. 50). An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 1.96$ GPa, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$.

Fig. 15

Rotation of a square fluid patch: Deformed fluid patch using a $\{\varvec{p}, J\}$ Upwind-SPH without dissipation terms (e.g. $\varvec{\mathcal {D}}_{\varvec{v}}^{ab} = \varvec{\mathcal {D}}_{\Delta \varvec{C}}^{ab} = \varvec{0}$ and $\mathcal {D}_p^{ab} = 0$) and b $\{\varvec{p}, J\}$ Upwind-SPH at $tA_0=0.324$. Results obtained using 6561 particles with an initial velocity and pressure (refer to Eq. 50). An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 1.96$ GPa, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$

7.3 Rotation of a square fluid patch

Similar to the example presented in Sect. 7.2, a square patch of inviscid fluid without the imposition of any external forces is considered. This problem was first proposed by Colagrossi [14], and this has since been used as an interesting benchmark validation example in the SPH community [68, 100, 101]. This problem is particularly challenging as spurious zero energy modes[8, 92] would accumulate and, eventually, lead to the breakdown of the scheme. One viable option to remove this instability is the use of kernel derivatives that is fixed in the reference configuration (inherent to the use of Lagrangian SPH formalism [41]), which is exactly one of the key motivations of this work.

The fluid patch, of length $L=1$ m, is subject to an initial velocity and pressure distribution described as

$$\begin{aligned}&\varvec{v}_0=\left[ \begin{array}{c} A_0Y \\ - A_0X \\ 0 \end{array} \right] ; \nonumber \\&p_0 = \rho _0 \sum _m^{\infty } \sum _n^{\infty } -\frac{32\omega ^2}{\left[ mn\pi ^2\left( \frac{m\pi }{L}\right) ^2 + \left( \frac{n\pi }{L}\right) ^2\right] } \nonumber \\&\quad \sin \left( \frac{m\pi X^*}{L}\right) \sin \left( \frac{n\pi Y^*}{L}\right) , \end{aligned}$$

(50)

where $A_0 = 200$ s$^{-1}$, $X^* = X + \frac{L}{2}$ and $Y^* = Y + \frac{L}{2}$. The parameters m and n are odd numbers and the pressure field (50b) converges rapidly when $\{ m, n \}$ $\ge$ 3. The centre of mass of this fluid patch is located at $\varvec{X}=\left[ 0, 0, 0\right] ^{\mathrm{T}}$ m and the initial density is taken as $\rho _0=1000$ kg/m$^3$. Using both the bulk modulus of $\kappa = 1.96$ GPa and the coefficient of $\gamma =1$, the pressure wave speed at time $t = 0$ becomes $c_p^0 = 1400$ m/s.

Fig. 16

Rotation of a square fluid patch: evolution of the deformation alongside with pressure distribution using $\{ \varvec{p}, J \}$ Upwind-SPH at (left to right) $tA_0= \{ 0.085$, 0.268, 0.594, 0.903, 1.466, 2.416, $3.570 \}$. Results obtained using 6561 particles with an initial velocity and pressure (refer to Eq. 50). An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 1.96$ GPa, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$

Fig. 17

Swirling of a square fluid patch: geometry configuration and initial conditions

Fig. 18

Swirling of a square fluid patch: particle refinement of deformed fluid patch alongside with pressure distribution plotted as particles (top) and contour (bottom) at $tA_0=0.2$ using a 441, b 1681, c 6561 and d 25,921 particles. Grey dashed line represents the initial position and red dashed line refers to the free surface deformation from ANSYS Fluent. Results obtained using $\{ \varvec{p}, J \}$ Upwind-SPH with an initial velocity and pressure distributions (refer to Eq. 51). An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 1.96$ GPa, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$. (Color figure online)

Fig. 19

Time evolution of a relative deformation of the non-dimensional corner particle at the material point $\varvec{X}=\left[ 0.5, 0.5, 0\right] ^{\mathrm{T}}$ m with its position shown in ‘$\square$’ up until $tA_0=0.5$; and b numerical dissipation 441, 1681, 6561 and 25,921 particles. Results obtained using $\{ \varvec{p}, J \}$ Upwind-SPH with an initial velocity and pressure (refer to Eq. 51). An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 1.96$ GPa, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$

Fig. 20

Swirling of a square fluid patch: Deformed fluid patch at $tA_0=0.5$ using a $\{\varvec{p}, J\}$ Upwind-SPH with $\varvec{F}$ ensures only zeroth order completeness and b $\{\varvec{p}, J\}$ Upwind-SPH with $\varvec{F}$ ensures both zeroth- and first-order completeness (31) Upwind-SPH. Results obtained using 6561 particles with an initial velocity and pressure distributions (refer to Eq. 51). An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 1.96$ GPa, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$

Fig. 21

Swirling of a square fluid patch: Deformed fluid patch using a $\{\varvec{p}, J\}$ Upwind-SPH by setting $\varvec{\mathcal {D}}_{\varvec{v}}^{ab} = \varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab} = \varvec{0}$ and $\mathcal {D}_p^{ab} = 0$, b $\{\varvec{p}, J\}$ Upwind-SPH by setting $\mathcal {D}_{p}^{ab} = 0$ and c $\{\varvec{p}, J\}$ Upwind-SPH at $tA_0=0.256$. Results obtained using 6561 particles with an initial velocity and pressure distributions (refer to Eq. 51). An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 1.96$ GPa, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$

Fig. 22

Swirling of a square fluid patch: Deformed fluid patch using a $\{\varvec{p}, J\}$ Upwind-SPH with $\mathcal {D}_p^{ab} = 0$ and b $\{\varvec{p}, J\}$ Upwind-SPH at $tA_0=1.1$. Results obtained using 6561 particles with an initial velocity and pressure distributions (refer to Eq. 51). An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 1.96$ GPa, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$

Fig. 23

Swirling of a square fluid patch: evolution of the deformation alongside with pressure distribution using $\{ \varvec{p}, J \}$ Upwind-SPH at (left to right) $tA_0= \{0.0$, 0.20, 0.51, 0.80, 1.01, $1.30 \}$. Results obtained using 25,921 particles with an initial velocity and pressure (refer to Eq. 51). An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 1.96$ GPa, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$

Fig. 24

Collapsing water column: problem setup

Fig. 25

Collapsing water column: particle refinement of deformed fluid patch alongside with pressure distribution plotted as contour at $t \sqrt{g/L}_0=0.9$ using a 1331, b 9261 and c 68,921 particles. Red dashed line represents the free surface deformation from ANSYS Fluent and the square markers represent the free surface solutions from Kelecy and Pletcher [105]. Results obtained using $\{ \varvec{p}, J \}$ Upwind-SPH with an initial hydrostatic pressure. An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 16$ MPa, material coefficient $\gamma =7$ and $\alpha _\text {CFL}=0.3$. (Color figure online)

By imposing zero pressure gradient at the four corners of the fluid patch, the particles initially located at these corners would evolve in time in the direction of the velocity field [102]. To examine this, a particle refinement analysis at time $tA_0=1.1$ is presented in Fig. 11. Both the free surface deformation and the pressure representation for every level of particle refinement are displayed. For benchmarking purposes, we simulate the exact same problem using the commercial software package ANSYS Fluent via an ultrafine discretisation of 2 million elements, which can then be treated as a reference solution. The proposed method accurately captures the location of the corners, as well as the bending of four arms, even with the utilisation of very few particles. No spurious mechanism is observed, which typically appears between $tA_0=0.6$ and $tA_0=1.2$ as reported in [68, 102].

Figure 12a monitors the trajectory path of the particle at $\varvec{X}=\left[ -0.5, 0.5, 0\right] ^{\mathrm{T}}$ m. After a relatively long time, the position of the particle becomes inaccurate due to the accumulation of the numerical dissipation. This is commonly found in the work of the other authors [14, 103] who attempted to simulate this example via a stabilisation procedure. However, dramatic improvement is seen with increased particle density. For completeness, the total numerical dissipation of the algorithm is measured. Figure 12b shows a clear reduction in the numerical dissipation by increasing the level of refinement.

Figure 13a, b show the time history of the components of the global linear and angular momenta of the system. As expected, the angular momentum components are perfectly conserved and the linear momentum components fluctuate around zero machine accuracy. Crucially, Fig. 14 shows the necessity for accurate evaluation of the deformation gradient $\varvec{F}$. Wrong deformation path is observed when employing the standard kernel gradient approximation $\varvec{\mathrm {C}}_{ab}$ [8] that only ensures zeroth-order consistency. Such inaccurate deformation can be rectified by utilising the correction technique to evaluate the deformation $\varvec{F}$ (31) (refer to Fig. 14b). Moreover, Fig. 15 highlights the importance of adding sufficient amount of numerical dissipation to the algorithm. As displayed in Fig. 15a, four corners of the fluid patch exhibit spurious mechanisms when using the $\{ \varvec{p}, J\}$ SPH algorithm without introducing any sort of stabilisation terms (i.e. set $\varvec{\mathcal {D}}_{\varvec{v}}^{ab} = \varvec{\mathcal {D}}_{\Delta \varvec{C}}^{ab} = \varvec{0}$ and $\mathcal {D}_p^{ab} = 0$). This shortcoming can be alleviated using suitable dissipation terms (see Fig. 15b), consistently derived from the satisfaction of global Hamiltonian energy.

For visualisation purposes, Fig. 16 shows a sequence of deformed states for the rotating fluid patch. Given the initial conditions described in (50), the fluid patch initially is dominated by a centrifugal force [99]. This would allow the fluid patch first to rotate and then to start developing the fluid arms at those four corners. The arms would continue to elongate and deform, and eventually creating four extremely thin bent arms at a later time. Very smooth pressure field is observed.

7.4 Swirling of a square fluid patch

Using exactly the same geometry as shown in Fig. 10a, another type of square fluid patch is assessed [98, 102, 104]. In this case, the fluid patch is subject to a specific representation for the velocity field described as follows

$$\begin{aligned} \varvec{v}_0= - A_0 L \left[ \begin{array}{c} e^{-4} - e^{-16(Y/L)^2} \\ e^{-4} + e^{-16(X/L)^2} \\ 0\\ \end{array} \right] , \end{aligned}$$

(51)

where $A_0 = 100$ s$^{-1}$. This particular type of velocity distribution enables the fluid patch to undergo extremely massive distortion. The adopted fluid properties are material density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 1.96$ GPa and the fluid coefficient of $\gamma =1$.

At time $t = 0$, the fluid domain is subjected to a pair of rotating vorticity fields pushing the fluid patch to deform along the flow direction. The two corner particles located at the centre of the vortices remain always at rest. The remaining corners of the domain are located at the flow symmetry direction, and hence known as free boundary particles. In order to show the convergence of the overall algorithm, a particle refinement study is carried out at time $tA_0=0.2$ and is shown in Fig. 18. Both the top and left surfaces of the fluid patch are being pushed towards the domain, whereas all the corner particles remain at rest. As time evolves, the free boundary particle located at position $\varvec{X}=[-0.5,-0.5,0]^{\mathrm{T}}$ m (i.e. the right bottom corner of the domain) gradually moves away from its original location. Insofar as the symmetric vorticity field is considered, the boundary particle propagates in a linear fashion. The movement evolution relative to its original position is depicted in Fig. 19a. Interestingly, the particle travels linearly with respect to its initial position for all particle refinements. The position of the particle converges with increasing particle density, which can be observed through the location of markers. Moreover, significant reduction in the amount of numerical dissipation is observed in Fig. 19b. This indeed shows that the numerical dissipation term is consistent and would vanish by increasing the number of particles.

Figure 20 demonstrates the importance of the proper evaluation of the deformation $\varvec{F}$. Clearly, the deformation path is wrongly captured without imposing any kernel gradient correction (refer to Fig. 20a) unless excessively fine meshes are used. As shown in Fig. 20b, this inaccurate deformed shape is corrected by using the enhanced kernel gradient approximation. It is well known that, sufficient amount of numerical stabilisation that needs to be added to the algorithm remains a persistent numerical issue in the SPH community. The SPH algorithm clearly exhibits pressure fluctuations without the introduction of any sort of stabilisation by setting the values of $\varvec{\mathcal {D}}_{\varvec{v}}^{ab} = \varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab} = \varvec{0}$ and $\mathcal {D}_p^{ab} = 0$. Addition of $\varvec{\mathcal {D}}^{ab}_{\varvec{v}}$ and $\varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab}$ to the linear momentum evolution alone would partially alleviate the mechansims at the beginning of the simulation Fig. 21b, c. The mechanism however would accumulate and re-emerge at a later time. As displayed in Fig. 22, this mode can be entirely removed by using the proposed SPH algorithm where $\{ \varvec{\mathcal {D}}^{ab}_{\varvec{v}}, \mathcal {D}^{ab}_p, \varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab} \}$ are appropriately introduced without the need to tune any user defined stabilisation parameters.

Finally, Fig. 23 shows a series of snapshots for the fluid patch, displaying the deformation patterns and the pressure field. A very smooth and stable pressure pattern is observed.

7.5 Collapsing water column

In this section, a classical three dimensional benchmark example is presented (see Fig. 24). A cubic water column, of length L = 0.05715 m, in hydrostatic equilibrium is confined between a wall and a gate. The oscillating water column is subject to a gravitational acceleration of $g=9.81$ m/s$^2$ acting vertically downwards. The material properties under consideration are density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 16$ kPa and fluid coefficient $\gamma =7$. Symmetric boundary conditions are imposed on the solid walls.

In this example, the gate is removed instantaneously and the fluid column collapses under the gravitational force. Specifically, a number of experimental and numerical results are available in the literature for benchmarking purposes. Figure 25 shows a refinement analysis at a particular time $t \sqrt{g/L}_0=0.9$, with special emphasis on the deformation patterns of the free surface alongside pressure distribution. With respect to the free surface deformation, the proposed upwind SPH algorithm agrees well with other published results (e.g. Kelecy and Pletcher [105] and ANSYS Fluent). Furthermore, the time varying positions at both the surge front and the column height of the wall are monitored and compared [2, 54, 105,106,107,108,109,110] in Fig. 26a, b. Again, very good agreement is observed even with very few particles.

Fig. 26

Collapsing water column: time evolution of a surge front and b column height using 1331, 9261 and 68,921 particles. Markers represent experimental results and dashed lines represent numerical solutions. Results obtained using $\{ \varvec{p}, J \}$ Upwind-SPH with an initial hydrostatic pressure. An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 16$ MPa, material coefficient $\gamma =7$ and $\alpha _\text {CFL}=0.3$

Essentially, an appropriate evaluation for the deformation map $\varvec{F}$ prevents an unrealistic collapsing trend of the water column (refer to Fig. 27). Most of the particles belonging to the surge regions are seen to be clumped together and some are seen to be protruding across the bottom wall. In addition, a dip is observed at column height.

Fig. 27

Collapsing water column: deformed fluid patch at $t\sqrt{g/L}_0=1.965$ using a $\{\varvec{p}, J\}$ Upwind-SPH with $\varvec{F}$ only ensure zeroth-order completeness and b $\{\varvec{p}, J\}$ Upwind-SPH with $\varvec{F}$ ensures both zeroth- and first-order completeness (31). Results obtained using 9261 particles with an initial hydrostatic pressure field. An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 16$ MPa, material coefficient $\gamma =7$ and $\alpha _\text {CFL}=0.3$

Figure 28 illustrates the deformation history of the water column, displaying pressure contour. As already reported in [47, 54, 80], pressure instabilities would be excited without incorporating any numerical damping to the Jacobian evolution (see Fig. 28a,b). This indeed can be eliminated by introducing a consistent Riemann-based dissipation term $\mathcal {D}^{ab}_p$ (refer to Fig. 28c), without affecting the order of convergence of the algorithm. Another popular option to address this drawback is to employ the $\delta$-SPH method [46, 67], and which requires a user defined stabilisation parameter for the evaluation of the Laplacian operator.

Fig. 28

Collapsing water column: evolution of the deformation (top to bottom) alongside pressure distribution using a $\{\varvec{p}, J\}$ Upwind-SPH by setting $\varvec{\mathcal {D}}_{\varvec{v}}^{ab} = \varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab} = \varvec{0}$ and $\mathcal {D}_p^{ab} = 0$, b $\{\varvec{p}, J\}$ Upwind-SPH with $\mathcal {D}_{p}^{ab} = 0$ and c $\{\varvec{p}, J\}$ Upwind-SPH at $t \sqrt{g/L}_0= \{ 0.39$, 0.59, 0.98, 1.46, 1.96, $2.49 \}$. Results obtained using 9261 particles with an hydrostatic pressure. An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 16$ MPa, material coefficient $\gamma =7$ and $\alpha _\text {CFL}=0.3$

7.6 Impact of a fluid patch on a wall

The final example of this work is the numerical analysis of an inviscid fluid patch impacting on a wall. The setup of this problem is very similar to the work reported in References [67, 68, 77] when considering impact of two identical fluid patches. A cuboid, with a height of $H=1.5$ m (or $L = 1$ m in Fig. 29a) and a unit cross section, is subject to an initial dropping velocity field

$$\begin{aligned} \varvec{v}_0= A_0 L \left[ \begin{array}{c} 0 \\ -1 \\ 0 \\ \end{array} \right] , \end{aligned}$$

(52)

where $A_0 = 1$ s$^{-1}$. The fluid properties used in this example are density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 10$ MPa and the fluid coefficient $\gamma =1$, which in turn yields the corresponding pressure wave speed $c_p^0 = 100$ m/s. Moreover, symmetric boundary conditions are enforced at the bottom surface, whereas the remaining boundary surfaces are set free.

Fig. 29

Impact of a cuboid fluid patch on a wall: problem setup

We begin this example by performing a particle refinement analysis at time $tA_0 = 0.098$ (see Fig. 30). The deformed shape of the free surface predicted by the proposed SPH method agrees extremely well with the ANSYS Fluent solutions via an ultrafine discretisation of 1.6 million elements. Pressure resolution is improved by increasing the number of particles. The results are further analysed by monitoring the displacement at two different positions, one at a surge front $\varvec{X} = [0.5,0,0]^{\mathrm{T}}$ m (see Fig. 31a) and the other one at a column height $\varvec{X} = [0,1.5,0]^{\mathrm{T}}$ m (see Fig. 31b). It is seen that the time history of the horizontal displacement at the surge front location converges to the reference solution (i.e. ANSYS Fluent) when increasing the number of particles. More interestingly, the vertical displacement at the column height is converged even with the use of a small number of particles, showing optimal convergence of the algorithm. Very small oscillations in column height being observed due to the nature of the compressible algorithm.

Fig. 30

Impact of a cuboid fluid patch on a wall: particle refinement of deformed fluid patch alongside with pressure distribution plotted as contour at $t A_0=0.098$ using a 2501, b 9801, c 38,801 and d 154,401 particles. Red dashed line represents the free surface deformation from ANSYS Fluent. Results obtained using $\{ \varvec{p}, J \}$ Upwind-SPH with an initial velocity $\varvec{v_0} = A_0[0,-1,0]^{\mathrm{T}}$ m/s using $A_0=1$ $s^{-1}$ and pressure $p_0 = 0$ Pa. An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 10$ MPa, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$. (Color figure online)

Fig. 31

Impact of a cuboid fluid patch on a wall: time evolution of a surge front and b column height using 2501, 9801, 38,801 and 154,401 particles. Black dashed line represents the numerical solutions for validation purposes. Results obtained using $\{ \varvec{p}, J \}$ Upwind-SPH with an initial velocity $\varvec{v_0} = A_0[0,-1,0]^{\mathrm{T}}$ m/s using $A_0=1$ s$^{-1}$ and pressure $p_0 = 0$ Pa. An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 10$ MPa, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$

Figure 32 depicts a series of snapshots of the pressure evolution simulated using the proposed Upwind-SPH method. These snapshots are captured at every time interval of $\Delta tA_0 = 0.001$ up until $tA_0 = 0.03$. After the fluid patch is in contact with the wall, a compressive shock wave is instantaneously generated and then followed by a tensile wave. These waves would then propagate towards the top free surface and reflect back to the fluid domain. A smooth pressure distribution is observed throughout the evolution without the need to introduce extra artificial (or non-consistent) viscosity into the algorithm. To check if the method suffers from long-term instability, we continue running the simulation for a longer time. As shown in Fig. 33, extremely good results are obtained via a Riemann based SPH algorithm.

Fig. 32

Impact of a cuboid fluid patch on a wall: time evolution of pressure monitored at material point $\varvec{X} = [0,0,0]^{\mathrm{T}}$ m. The blue dots represent the time instant (every 0.001 s) of the deformation alongside with pressure distribution (left to right; top to bottom) using 102,541 particles. Results obtained using $\{ \varvec{p}, J \}$ Upwind-SPH with an initial velocity $\varvec{v_0} = A_0[0,-1,0]^{\mathrm{T}}$ m/s using $A_0=1$ $s^{-1}$ and pressure $p_0 = 0$ Pa. An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 10$ MPa, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$. (Color figure online)

Fig. 33

Impact of a cuboid fluid patch on a wall: evolution of the deformation (top to bottom) alongside with pressure distribution using a $\{\varvec{p}, J\}$ Upwind-SPH with $\varvec{\mathcal {D}}_{\varvec{v}}^{ab} = \varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab} = \varvec{0}$ and $\mathcal {D}_p^{ab} = 0$, b $\{\varvec{p}, J\}$ Upwind-SPH with $\mathcal {D}_p^{ab} = 0$ and c $\{ \varvec{p}, J \}$ Upwind-SPH at $t A_0= \{0.00487$, 0.0355, 0.1355, 0.255, 0.3672, 0.549, $0.7935 \}$. Results obtained using 102,541 particles with an initial velocity $\varvec{v_0} = A_0[0,-1,0]^{\mathrm{T}}$ m/s using $A_0=1$ s$^{-1}$ and pressure $p_0 = 0$ Pa. An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 10$ MPa, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$

To further examine the robustness of the algorithm, we simulate the exact same problem but this time applied on a hemisphere. Figure 34 shows the importance of adding sufficient amount of numerical dissipation to the volume map conservation law. The proposed upwind SPH algorithm performs extremely well in this case without any numerical difficulties.

Fig. 34

Impact of a hemisphere fluid patch on a wall: evolution of the deformation (top to bottom) alongside with pressure distribution using a $\{\varvec{p}, J\}$ Upwind-SPH with $\varvec{\mathcal {D}}_{\varvec{v}}^{ab} = \varvec{\mathcal {D}}_{\Delta \varvec{\mathrm {C}}}^{ab} = \varvec{0}$ and $\mathcal {D}_p^{ab} = 0$, b $\{\varvec{p}, J\}$ Upwind-SPH with $\mathcal {D}_p^{ab} = 0$ and c $\{ \varvec{p}, J \}$ Upwind-SPH at $t A_0= \{0.003$, 0.0607, 0.1310, 0.22, 0.303, 0.395, $0.48 \}$. Results obtained using 136,393 particles with an initial velocity $\varvec{v_0} = A_0[0,-1,0]^{\mathrm{T}}$ m/s using $A_0=1$ s$^{-1}$ and pressure $p_0 = 0$ Pa. An elastic fluid of density $\rho _0=1000$ kg/m$^3$, bulk modulus $\kappa = 10$ MPa, material coefficient $\gamma =1$ and $\alpha _\text {CFL}=0.3$

8 Conclusions

In this paper, a new Total Lagrangian Smooth Particle Hydrodynamics computational framework has been introduced for the numerical analysis of inviscid fluid flows. The methodology is established starting from a mixed-based system of first-order hyperbolic conservation laws, where the linear momentum $\varvec{p}$ conservation law is solved along with the conservation equation for the Jacobian of the deformation J. For problems without experiencing severe changes in topological surface, the use of a Total Lagrangian description removes the need for continuous (updated) particle neighbouring search (either removing it altogether or delaying it until an extremely large number of computational time steps has taken place), yielding an extremely efficient scheme still capable of handling massive deformations.

From the spatial discretisation point of view, an entropy-stable Smooth Particle Hydrodynamics is presented where consistent Riemann-based (upwinding) numerical dissipation is suitably introduced guaranteeing global numerical entropy production, the latter demonstrated by monitoring the so-called Hamiltonian energy of the system. Such numerical dissipation is very well designed from a mathematical standpoint (taking into consideration of the characteristic structure of the hyperbolic system) and, crucially, does not rely on the use of any user defined artificial stabilisation parameters. From the temporal discretisation point of view, an explicit one-step two-stage Total Variation Diminishing Runge Kutta time integrator combined with an angular momentum preserving algorithm is presented.

Finally, a comprehensive list of challenging prototype problems has been presented in order to benchmark the results obtained against alternative numerical strategies available in the literature, including the standard SPH [8, 32, 40], the Corrected SPH [54, 55] and ANSYS Fluent [110]. It has been shown that the resulting SPH algorithm overcomes a number of numerical difficulties posed by SPH methods when simulating slightly compressible fluids, namely non-physical hydrostatic pressure fluctuations, hour-glassing and numerical errors associated with global conservation, completeness, long-term instability and convergence. Very importantly, both velocities and pressure (or the volumetric strain) display equal second order of convergence. As a result, the new SPH algorithm provides a good balance between accuracy and speed of computation.

The next step of our work is the adaptation of the current Total Lagrangian SPH algorithm to handle violent free-surface flows typically accompanied with large topological changes, by incorporating alternative Updated Lagrangian formalism in the style of [28].

Hyperbolicity

This Appendix is included to demonstrate the hyperbolicity [85, 86] of system (8). In addition, the computation of the maximum (pressure) wave speed is necessary for the evaluation of the maximum time step of the explicit time integrator. This can be achieved by, first recalling that $\varvec{P} = - p \varvec{H}$ (10) and, combining Eq. (6) multiplied by $\varvec{H}^{\text {Ave}} \varvec{N}$ with Eq. (4) to give

$$\begin{aligned} U^2 = - \frac{\Lambda _A^2 }{\rho _0} \frac{\llbracket p \rrbracket }{\llbracket J \rrbracket }; \qquad \Lambda _A^2 = \left( \varvec{H}^{\text {Ave}} \varvec{N} \right) \cdot \left( \varvec{H}^{\text {Ave}} \varvec{N} \right) . \end{aligned}$$

(53)

In this paper, we consider the simple case of an “elastic fluid” type of model described in (10), that is $p = \kappa (J^{-\gamma } - 1)$. In this case

$$\begin{aligned} \llbracket p \rrbracket = \kappa \llbracket J^{-\gamma } \rrbracket . \end{aligned}$$

(54)

In the case of slightly compressible fluid, it is instructive to consider small jumps in $\llbracket J^{-\gamma } \rrbracket$ by approximating it using a Taylor series expansion as

$$\begin{aligned} \llbracket J^{-\gamma } \rrbracket \approx -\frac{\gamma \llbracket J \rrbracket }{(J^{\text {Ave}})^{\gamma +1}};\qquad J^{\text {Ave}}=\frac{1}{2}(J^{+}+J^{-}). \end{aligned}$$

(55)

Combining (55) with the second line of (54) would give an alternative expression for the jump in pressure described as

$$\begin{aligned} \llbracket p \rrbracket = -\frac{\kappa \gamma }{(J^{\text {Ave}})^{\gamma +1}}\llbracket J \rrbracket . \end{aligned}$$

(56)

Finally, the determination of pressure wave speed can now follow by substituting Eqs. (56) into (53), which results in

$$\begin{aligned} U = \pm c_p; \qquad c_p = \Lambda _A \sqrt{ \frac{\gamma \kappa }{ \rho _0 (J^{\text {Ave}})^{\gamma +1}}} . \end{aligned}$$

(57)

It is worth noticing that expression above (57) degenerates to the so-called linear acoustic pressure wave speed $c_{p}^{\text {Lin}}$ by evaluating the nonlinear pressure wave speed $c_p$ (57) at the initial undeformed configuration (that is, $\varvec{H} = \varvec{I}$ and $J = 1$)

$$\begin{aligned} c_{p}^{\text {Lin}} = \sqrt{\frac{\gamma \kappa }{\rho _0}}. \end{aligned}$$

(58)