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

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.


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].
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 δ-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 1 is presented in terms of the linear momentum p and the Jacobian J of the deformation. The mixed { 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 { 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.

Kinematics and conservation of volume map
Consider the three dimensional deformation of a continuum moving from its initial (material) configuration occupying a volume 0 , of boundary ∂ 0 with outward unit normal N, into a current configuration at time t occupying a volume (t) of boundary ∂ (t) with outward unit normal n, see Fig. 1. The standard notation and definitions for the deformation gradient tensor F (used to map differential fibre elements d x = Fd X) and its determinant J (used to map differential volume elements d = J d 0 ) are introduced where x represents the current position of a particle originally at X and ∇ 0 denotes the gradient with respect to the material configuration. Similarly, material differential area vectors d A (co-linear with N) are mapped to current differential area vectors d a (co-linear with n) through the co-factor H of the deformation, which is related to the deformation gradient tensor F via the Nanson's rule [86] 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, where v = dx dt is the velocity field. The equivalent local conservation law and associated jump condition across a discontinuity surface with normal N, propagating with speed U in the reference space [78,83,84], are In this expression, DIV represents the material divergence operator, H Ave = 1 2 H − + H + represents an average state of the co-factor of the deformation between the left and right states across the discontinuity surface and · = [·] + − [·] − represents the jump operator.

Conservation of linear momentum
The conservation of linear momentum per unit undeformed volume p = ρ 0 v (with ρ 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 0 by where f 0 is a body force per unit of undeformed volume and t = P N is the traction vector associated with the material outward unit normal surface vector N with 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

Combined equations
Combining expressions (3) and (5), above global conservation laws can now be summarised in a concise manner as representing a system of Total Lagrangian conservation laws where U denotes the set of conservation variables, S the source term and F N the surface flux vector, described as follows where H is the cofactor of the deformation previously defined in (2). The computation of 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 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 P with the Jacobian of the deformation J . In this specific case, a strain energy per unit undeformed volume is introduced as with {γ, κ} > 0 the positive material constants. 2 Thus, the first Piola-Kirchhoff stress tensor P can be expressed as 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 P v , where the combined stress becomes 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 P v = 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 DIVH = 0, into Eq. (4a) to give

Weak form statements
In general, a standard weak variational statement for the system (8) is established by multiplying the local differential equations U (8) with their appropriate work conjugate virtual fields δV, and integrating over the reference domain 0 of the body, to give where the normal fluxes F N = F I N I with N I being the material outward normal in I -th material direction and the symbol • 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 δ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 H = H(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 H is defined by which represents the summation of the kinetic energy (i.e. the first term on the right hand side of (14)) and Helmholtz free energy per unit of undeformed volume. Note here that H(X, t) andĤ(U ) represent alternative functional representations of the same quantity. With expression (14) at hand, the work conjugates V can then be obtained by taking simple differentiation to give Finally, upon substitution of (15) into (13), individual variational statement for the linear momentum p and the Jacobian of the deformation J are Here, both δv and δ p represent the virtual velocity and pressure fields, respectively.

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 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 Recalling that v · DIV P + P : ∇ 0 v = DIV( P T v), above equation reduces to By performing integration by parts of the DIV term in Eq. (19) and after some re-arrangement, it yields where˙ ext denotes the power introduced by external forces, defined aṡ In expression (20) , the term on the right hand side represents possible physical dissipation and is always non-positive (namely, − P v : ∇ 0 v ≤ 0) and, consequently, Eq. (20) can be transformed into the following inequality 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.˙ ext = 0), inequality (22) reduces to 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.

Spatial discretisation
For the case of corrected SPH methods [54], both the problem variables U and the conjugate pairs δV are in general interpolated at any given position via corrected SPH shape functions N (or kernel functionsW ) with a given compact support of radius 2h around every particle (see Fig. 2). Specifically, for a given position X a , both U and δV can be approximated as Here, b a represents the set of neighbouring particles b that lie inside a sphere of a given radius around X a , b 0 represents the volume associated to particle b and U b (t) and δV b represent the time varying variables and their virtual conjugate  t). Additionally, the use of corrected SPH shape functionsÑ ensures that both constant and linear functions are perfectly interpolated [54].
In line with previous work [31], the material gradient ∇ 0 (•) of any arbitrary vector function f can be approximated through the so-called corrected pseudo area vector C ab (defined on the basis of corrected kernel gradient∇ 0 [54]), described as Specifically, the correction technique involves modifying (uncorrected) pseudo area vector C ab (or kernel gradient ∇ 0 ) by introducing a particle based correction matrix L a [54] C ab = L a C ab ; from which L a is explicitly evaluated as [54,55] This type of kernel gradient correction [28,54] ensures the gradient of any linear field distribution is exactly evaluated. Moreover, the term − 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. L a = I), the corrected pseudo area vectorC ab (26) degenerates to the usual SPH uncorrected area vector 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 L a = L b . This would consequently destroy the anti-symmetric condition of the gradient evaluations (i.e.C ab = −C ba ) [92]. One possible option to enforce such condition consists of skew-symmetrising the area vec-torC Skew ab by taking the difference between the area vectors of particle pairs to givẽ

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 { p, J } described in (7) becomes In above expressions, the pairwise interacting force clearly satisfies local conservation by construction, namely (2) is evaluated using the discrete deformation gradient F a accurately approximated as where x a represent the current nodal geometry at particle a obtained from time integration of the velocity v a . Finally, the remaining terms to be defined in Eq. (29a, b) are the so-called pairwise stabilisation terms In particular, the term D v ab in (29a) addresses spurious zero-energy (hourglass-like [33]) modes due to rank-deficiency, whereas the term D p ab (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 D ab C into expression (29a), with the objective to penalise the difference between the usual uncorrected pseudo area vector C ab andC Skew ab (28). Mathematically, this is described as C ab −C Skew ab . 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 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.

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 Subsequently, we can substitute the linear momentum equation (29a) and the Jacobian equation (29b) into (32), and after some algebra, gives Here,˙ ext denotes the semi-discrete power introduced by external forces expressed aṡ 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 It is now the objective to demonstrate that the term on the right hand side of (35) is always non-positive, that is D total ≥ 0 (to be in agreement with (22)). Specifically, 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.
, an alternative expression for D total is Dissipation terms {D ab v , D ab p , D ab C } remain to be defined in order to ensure non-negative entropy production, namely, D total > 0. In this work, the term D ab C is defined as 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 D ab C ensures the skew symmetric nature of this term. As for the remaining stabilisation terms, we propose the following definitions (40) are skew symmetric and proportional to the difference (or jump) in velocity (v b − v a ) and pressure ( p b − p a ) between pairwise particles, typical of Riemann solver based upwinding terms [78]. Moreover, with these definitions for {D ab v , 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 with respectively. Here, ρ R,a denotes the material density evaluated at particle a and c p,a represents the pressure wave speed (58) computed at particle a andc Skew 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 Ave p,ab in (42). For instance, selecting S p ab = 0, the pairwise numerical dissipation becomes wherec p,ab is an appropriate wave speed which can be selected as thus ensuring positive entropy production.

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 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 { p a , J a , x a } (a = 1, . . . , N , where N represents the total number of particles across the computational domain) are all updated via (45). The maximum time step t is governed by the standard Courant-Friedrichs-Lewy (CFL) condition given as where c p,max is the maximum p-wave speed, h min is the minimum particle spacing within the computational domain and α CFL is the CFL stability number. For the numerical computations presented in this paper, a value of α 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 T ab (30) is not co-linear with the position vector (x b − x a ) [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 b∈ b a T ab and the particle based stabilisation terms, namely b∈ b a D ab v and b∈ b a D ab 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. (4) COMPUTE right-hand-side of the mixed-based system: p a (29a) (with T ab described in (30)) andJ a (29b) (5) APPLY discrete angular momentum preserving algorithm (6) EVOLVE {U a , x a } via TVD-RK (45) (7) COMPUTE first Piola P a ((10)) end

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.

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 { 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 v(X, t) = ω A 0 cos (ωt) V ; where ω = 2π s −1 and A 0 = 5 × 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 An elastic fluid type of constitutive law is used with density ρ 0 = 1000 kg/m 3 , bulk modulus κ = 10 MPa and γ = 1.
As presented in 58, the initial pressure wave speed used in this case is c 0 p = κγ ρ 0 = 100 m/s. 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.

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 v 0 and the initial pressure field p 0 are, respectively, described as follows 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 ρ 0 = 1000 kg/m 3 and bulk modulus κ = 1.96 GPa.
Utilising the value of γ = 1 in the equation of state, the initial pressure wave speed of the fluid patch can be computed as c 0 p := κγ ρ 0 = 1400 m/s. The domain is spatially dis-cretised with a uniformly distributed particle arrangement. Moreover, the outermost radial particles are set free in order to allow the boundary surface to deform freely. 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%.
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 t A 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].
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 (X = [1, 0, 0] T m)    Table 2 Stretching of a circular fluid patch: comparison of position for farthest particle along the (X = [1, 0, 0] T m) and semi-major (X = [0, 1, 0] T m) axes between sequential particle density using { p, J } Upwind-SPH against analytical solution at t A 0 = The percentage difference against analytical solution for each particle density is stated in bracket  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 X = [0, 0, 0] T using three different values of the bulk modulus, namely κ, 4κ and 16κ. As shown in Fig. 8, all the results converge towards the incompressible solution without the appearance of non-physical pressure instabilities.
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 t A 0 = 2. To demonstrate the stability of the algorithm in a long-term response, we continue running the simulation up to t A 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.

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 where A 0 = 200 s −1 , X * = X + L 2 and Y * = Y + L 2 . The parameters m and n are odd numbers and the pressure field (50b) converges rapidly when {m, n} ≥ 3. The centre of mass of this fluid patch is located at X = [0, 0, 0] T m and the initial density is taken as ρ 0 = 1000 kg/m 3 . Using both the bulk modulus of κ = 1.96 GPa and the coefficient of γ = 1, the pressure wave speed at time t = 0 becomes c 0 p = 1400 m/s. By imposing zero pressure gradient at the four corners of the fluid patch, the particles initially located at these cor-    [102]. To examine this, a particle refinement analysis at time t A 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 t A 0 = 0.6 and t A 0 = 1.2 as reported in [68,102]. Figure 12a monitors the trajectory path of the particle at X = [−0.5, 0.5, 0] 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 F. Wrong deformation path is observed when employing the standard kernel gradient approximation C ab [8] that only ensures zeroth-order consistency. Such inaccurate deformation can be rectified by utilising the correction technique to evaluate the deformation F (31) (refer to Fig. 14b). More-over, 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 { p, J } SPH algorithm without introducing any sort of stabilisation terms (i.e. set D ab v = D ab C = 0 and D ab p = 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.

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 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 ρ 0 = 1000 kg/m 3 , bulk modulus κ = 1.96 GPa and the fluid coefficient of γ = 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 t A 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 X = [−0.5, −0.5, 0] 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 F. Clearly, the deformation path is wrongly captured without imposing any kernel gradient cor-  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 D ab v = D ab C = 0 and D ab p = 0. Addition of D ab v and D ab C 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 {D ab v , D ab p , D ab C } 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.

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 ρ 0 = 1000 kg/m 3 , bulk modulus κ = 16 kPa and fluid coefficient γ = 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 √ 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.
Essentially, an appropriate evaluation for the deformation map 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. 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 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 δ-SPH method [46,67], and which requires a user defined stabilisation parameter for the evaluation of the Laplacian operator.

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 where A 0 = 1 s −1 . The fluid properties used in this example are density ρ 0 = 1000 kg/m 3 , bulk modulus κ = 10 MPa and the fluid coefficient γ = 1, which in turn yields the corresponding pressure wave speed c 0 p = 100 m/s. Moreover, symmetric boundary conditions are enforced at the bottom surface, whereas the remaining boundary surfaces are set free.
We begin this example by performing a particle refinement analysis at time t A 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  Fig. 31a) and the other one at a column height X = [0, 1.5, 0] 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. 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 t A 0 = 0.001 up until t A 0 = 0.03. After the fluid patch is in contact with the wall, a compressive shock wave is instan-taneously 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.
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.

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 p conservation law is solved along with the con-servation 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 entropystable 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 freesurface flows typically accompanied with large topological changes, by incorporating alternative Updated Lagrangian formalism in the style of [28]. In this paper, we consider the simple case of an "elastic fluid" type of model described in (10), that is p = κ(J −γ − 1). In this case In the case of slightly compressible fluid, it is instructive to consider small jumps in J −γ by approximating it using a Taylor series expansion as Combining (55) with the second line of (54) would give an alternative expression for the jump in pressure described as Finally, the determination of pressure wave speed can now follow by substituting Eqs. (56) into (53), which results in It is worth noticing that expression above (57) degenerates to the so-called linear acoustic pressure wave speed c Lin p by evaluating the nonlinear pressure wave speed c p (57) at the initial undeformed configuration (that is, H = I and J = 1)