High-Order Eulerian Incompressible Smoothed Particle Hydrodynamics with Transition to Lagrangian Free-Surface Motion

5 The incompressible Smoothed Particle Hydrodynamics (ISPH) method is derived in Eulerian form with high-order smoothing kernels to provide increased accuracy for a range of steady and transient internal flows. Periodic transient flows, in particular, demonstrate high-order convergence and accuracies approaching, for example, spectral mesh-based methods. The improved accuracies are achieved through new high-order Gaussian kernels applied over regular particle distributions with time stepping formally up to 2 order for transient flows. The approach can be easily extended to model free surface flows by merging from Eulerian to Lagrangian regions in an Arbitrary-Lagrangian-Eulerian (ALE) fashion, and a demonstration with periodic wave propagation is presented. In the long term, it is envisaged that the method will greatly increase the accuracy and efficiency of SPH methods, while retaining the flexibility of SPH in modelling free surface and multiphase flows.

drodynamics, originally due to Dilts [20], minimise a least squares functional to reconstruct (potentially 48 high-order) polynomials based on surrounding particle values. The approach is once again implicit (and 49 likely to be costly at higher orders and dimensions), but notable improvements in accuracy have been 50 demonstrated with this approach in weakly compressible SPH when combined with a WENO scheme 51 and an appropriate Riemann solver [3]. Particle regularisation procedures, such as shifting at each time 52 step [66, 38] and particle re-meshing [14, 17] help by improving particle distribution and thereby reduce 53 error due to non-uniformity, but perfect convergence is not recovered (when it has been considered). 54 For example, studies often show observed convergence rates in the primitive variables lying between 1.5 55 and (less than) 2 [38]. A review of particle methods with some focus on re-meshing (and over multiple 56 length scales) is given in [31]. Very recent work by Litvinov et al. [39] undertakes an iterative particle 57 regularisation procedure (akin to shifting) that is able to recover particle distributions that show the ideal 58 convergence in the discretisation error. However, such an approach is likely to be costly for transient 59 simulations as iterative regularisation will need to be undertaken at each time step. As an alternative to 60 the aforementioned techniques, this paper will demonstrate that Eulerian SPH with high-order kernels 61 can offer straightforward improvements in convergence and accuracy quite efficiently for internal flows. 62 Importantly, use of ESPH does not mean the attractive features of Lagrangian SPH are lost. It shall also 63 be shown that the fully Lagrangian form can be coupled in a straightforward way where necessary (e.g. 64 near free surfaces).

SPH Interpolation and High-Order Kernels
In SPH, a variable A at a point r i is approximated by a convolution product of the variable A with a 81 smoothing kernel function ω h (| r i − r |), with a smoothing length h, and is written as where Ω is the supporting domain. In a discretised format, the interpolation can be written as where V j is the particle volume, and r ij the distance between particle i and j.
Hereafter ω h (r ij ) will 84 be simply written as ω ij or ω in the discrete and continuous cases, respectively. In addition to the 85 requirement that ω tends to the Dirac delta function as h → 0, standard SPH usually advocates the where A i = e r ·∇A i denotes the directional derivative along r at point i, with e r the unit direction vector.

93
Substituting (5) into (1) and non-dimensionalising kernel length scales with respect to the characteristic 94 smoothing length, h, yields where s = r/h. If condition 1 is obeyed and the kernel is normalised then the first term in expansion (6) 96 is identically A i . If the kernel is symmetric then all terms where ω is multiplied by an odd power of s 97 are zero as the integrand is an odd function with an integral of zero over a radially symmetric domain.

98
Consequently the error associated with smoothing/interpolation is O(h 2 ). Removal of condition 2 (pos-99 itivity) allows one to attain arbitrary (but even) orders of accuracy in smoothing error by constructing 100 extended kernels that satisfy the higher order terms (or moments) in the Taylor expansion. As indicated 101 by [46], removal of the 3 rd term on the RHS of (6) can be achieved with a kernel of the form: 102 where the subscript 4 denotes a kernel with a smoothing error that is now O(h 4 ). Indeed, the argument 103 can be generalised such that all even moments up to 2n vanish with a kernel of the form: Gaussian kernel (on a regular particle distribution) results in a discretisation error that is small compared 115 to the smoothing error. For kernels with compact support and different decay properties, Quinlan et al.

116
[51] showed that the discretisation error depends critically on the ratio of particle spacing to smoothing 117 length and on the boundary smoothness,β, of the kernel (β being the largest integer such that theβ th 118 derivative and all lower derivatives are zero at the edge of the kernel support). The Gaussian is infinitely 119 smooth, while popular kernels such as the cubic and quintic spline haveβ = 2 andβ = 4, respectively.

120
Compact kernels can then recover the ideal smoothing error, provided that the boundary smoothness is 121 sufficient and h is adjusted appropriately with particle spacing. Practically then, provided the particles 122 remain uniform and the first order error due to non-uniformity is removed [51], ideal convergence can 123 be achieved at higher-order. With the aforementioned discussion in mind, the Gaussian is chosen as our 124 "base kernel" from which we construct higher-order counterparts to use in simulations. Removal of the 125 2 nd moment from (6) results in the 4 th order Gaussian (G4) in 2D: Similarly, removal of the 4 th moment gives a two-dimensional 6 th order Gaussian kernel (G6): One could continue to construct higher-orders should one wish, but it will be demonstrated that the G6 128 is quite sufficient for attaining high accuracies at moderate particle spacing. SPH approximates a set of governing partial differential equations through substitution of appropriate 132 gradient (and, if required, Laplacian) operators. Two popular gradient approximations in SPH are: In standard SPH, the negatively signed version maintains zeroth order consistency as the gradient of 134 any constant function A is zero. Meanwhile, the positive version is symmetric in the particle index and 135 thereby conserves particle-particle momentum. Consequently, the positive version of (11) is favoured by 136 those who wish to conserve properties associated with the discrete particle system over the immediate 137 gains in accuracy at the continuum level that arise from the negative version. As ESPH treats particles 138 as fixed interpolation points the choice of sign in (11) becomes largely irrelevant, and both versions have 139 similar errors for uniform particle spacing [23]. To align with our previous work, we take the negative 140 version of (11). For a kernel of order O(h 2n ), the dominant error term (and so convergence properties) 141 of the gradient approximation can be determined through a Taylor expansion and integration by parts 142 to give: where e ij is the unit direction vector between points i and j. In a similar argument to that used for the 163 gradient, it can be shown that the lead smoothing error for the Morris Laplacian, E ∇ 2 , is of equal order 164 to that of the SPH interpolation and gradient. The Morris Laplacian is, therefore, suitable for high-order 165 calculations. In particular, for a kernel of order 2n: Despite remaining O(h 2n ), compared with (12), the Laplacian error depends on higher derivatives of 167 function A. It is a natural requirement that, in the extension to higher orders, a greater degree of conti-168 nuity is assumed in the solution, which may be disrupted at boundaries of the computational domain, at 169 the boundaries of compactly supported kernels, or the truncated Gaussian domain. We expect, therefore, 170 the Laplacian to be more prone to error, in general, than the first derivative.
the L 1 error in the gradient and the Laplacian of f will be studied. Note, we will take the horizontal com-176 ponent of the gradient to calculate the error, but the vertical component error magnitude and behaviour 177 is identical. Figure 1 shows convergence in the error of e x · ∇f for the 2 nd order (G2), 4 th order (G4) 178 and 6 th order (G6) Gaussian-based kernels. All three kernels demonstrate near optimum convergence.

179
Note that, as we are using Gaussian-based kernels (with associated exponential convergence in the dis-180 cretisation error), the ratio of h to dx was fixed at h = 2dx. Furthermore, to ensure that any error due the simulation. Therefore, despite a large support radius, the ESPH computation can be as fast as its

198
Consider the non-dimensional unsteady Navier-Stokes equations in Eulerian form: Of course, as the computation points are now fixed in space, the advection term is made explicit and 200 calculated at each particle using standard SPH techniques (Eqn. (11)). The Reynolds number Re is 201 defined in the usual way, where ρ and μ are the fluid density and viscosity, respectively. advantage of using a projection-based incompressible SPH method over weakly-compressible SPH is that 208 the resulting pressure field is smooth and non-oscillatory, provided the particle distribution is reasonably Estimations of the advection term contained in f at time step n + 1 can be found through a second order 221 extrapolation, e.g.
The divergence correction of u * is undertaken in the second step: where ∇ · u n+1 = 0 and δp n+1 = p n+1 − p n . Note that combining (18) with (20) yields an appropriate 224 temporal discretisation of (16). Defining q = (δp) n+1 + 1 Re ∇ · u * , the following Poisson equation can be with a consistent Neumann condition ∂q ∂n = 0 on the boundary of the flow domain. Upon solution of 227 (21), the pressure at time n + 1 can be found from p n+1 = q n+1 + p n − 1 Re ∇ · u * . The above pressure 228 correction scheme has the advantage that for no-penetration boundaries, boundary conditions for the 229 pressure are imposed correctly (i.e. they are consistent with the momentum equation (16)). Indeed, In particular, the following dimensionless inertial and viscous constraints apply [50]: In the finite difference study of [24], the values of constants C I and C μ are taken to be 2 and 1/6, respec-239 tively. Here we adopt conservative estimates and choose C I = 1 and C μ = 0.1 for most cases. Note that 240 Eqn. (22) is not the standard CFL condition, but the slightly stricter criterion derived from discretisa-241 tions of the Navier-Stokes which include the hyperbolic advection term without special treatment [54].

242
At the start of each simulation, the smaller of the two time step constraints is taken.

244
In summary, utilising the discrete SPH operators of Section 3, the fully discrete system of governing 245 equations reads: (24) • The Poisson Equation • The 2nd (velocity correction) step used depends on the test case, but its effect on error and convergence will be discussed in due course.

Taylor-Green Vortices
where b = − 8π 2 Re . The flow consists of counter-rotating vortices decaying temporally due to viscosity; 259 Figure 3 shows a typical pressure contour plot with streamlines over a fixed Cartesian particle distribution.

260
All results in this subsection use a support radius of 6h with h = 2dx. integration. Figure 4 plots the L 1 error convergence in the velocity for each of the studied kernels. Both 267 G2 and G4 demonstrate optimum convergence while G6 begins to converge optimally before being lim-268 ited by machine precision. This demonstrates the potential gains in accuracy with high-order kernels 269 and controlled particle distributions. Evidently, the more practical simulations that follow will be limited Dirichlet condition is prescribed to a single corner particle using the analytical solution above.
273 Figure 4: Convergence in the numerical horizontal velocity after one time step (Δt = 1 × 10 −8 ) given analytical initial conditions. 2 nd order (circles), 4 th order (triangles), 6 th order (squares). The straight lines depict theoretical convergence. in velocity consistent with the second order time integration scheme (see Table 1). However, the error 284 in the pressure remains limited by the iterative solver. Only when the tolerance of the solver is reduced 285 from 1 × 10 −11 to 1 × 10 −12 does the pressure error reduce accordingly. While improvements to the 286 pressure error can be achieved through decreases in tolerance, for the purpose of practical computation 287 such small tolerances can be quite expensive. For the cases studied, while a tolerance of 1 × 10 −11 can 288 take less than 100 iterations, 1 × 10 −12 can require between 500-1000. Application and preference will 289 determine whether one wishes to recover 6 th order convergence and accuracies (beyond O(10 −8 ) here) 290 in the pressure, given the increased computational expense. For further reductions in time step, Table   291 2 demonstrates the continued second order convergence in velocity. While the pressure demonstrates 292 similar second order convergence in the early stages (as predicted), it is once again limited by the solver 293 tolerance at smaller time steps (Table 2). The velocity data from Table 2 is plotted in Figure 7, which 294 shows how reductions in time step begin to recover the ideal 6 th order spatial convergence. Recovery of 295 ideal 6 th order convergence is possible with appropriate reductions in tolerance and time step for this 296 case, but, in practical simulations, given the exceptional accuracies already provided (O(10 −12 )), there 297 is arguably little to gain given the additional computational expense.

298
Considering the longer-term behaviour of the method, at t = 10 the convergence properties and accuracy 299 remain very good. Tables 3, 4 and 5 show L 1 error in velocity and pressure for a selection of particle 300 spacings and the different kernels. Table 3 shows results from G2 with ideal 2 nd order convergence being 301 demonstrated in both u and p. Table 4 shows errors for the 4 th order Gaussian G4: ideal 4 th order 302 Figure 5: L 1 error convergence in the horizontal velocity for 2 nd order (circles), 4 th order (triangles), and 6 th order (squares) Gaussian-based kernel functions. The straight lines denote the theoretical ideal convergence. The measurement is taken at time t = 0.1. Figure 6: L 1 error convergence in the pressure for 2 nd order (circles), 4 th order (triangles), and 6 th order (squares) Gaussian-based kernel functions. The straight lines denote the theoretical ideal convergence. The measurement is taken at time t = 0.1.

Time step & tolerance
u L 1 error p L 1 error dt = 1 × 10 −3 and tol=1 × 10 −11 1.31 × 10 −11 6.08 × 10 −8 dt = 5 × 10 −4 and tol=1 × 10 −11 3.50 × 10 −12 9.10 × 10 −8 dt = 5 × 10 −4 and tol=1 × 10 −12 3.25 × 10 −12 2.41 × 10 −9  convergence is seen in the velocity, while error order for the pressure is near ideal initially before de-303 creasing slightly as the error limits due to time stepping and solver tolerance are approached. For G6 304 (Table 5), the accuracy of the spatial approximation is such that the lower bound on error due to time 305 integration error and tolerance is reached quickly. Convergence in velocity is near ideal initially before 306 decreasing to around 4 between dx = 0.0125 and dx = 6.25 × 10 −3 . The error in the pressure reaches 307 its limit almost instantly, with errors restricted to O(10 −6 ). Even at t = 10, in the presence of accrued 308 time integration/solver errors, the G6 kernel provides results that are between 100 and 1000 times more 309 accurate than the standard 2 nd order Gaussian for a given particle spacing.  Table 6 shows the L 1 error in velocity and pressure for the three Gaussian kernels 313 (G2, G4, G6) at Re = 10 and t = 0.  Table 3: L 1 Error in horizontal velocity and pressure at t = 10 for different particle spacings with associated orders of convergence. Results are obtained with 2 nd order Gaussian, G2.

336
While an important test of spatial and temporal accuracy, the Taylor-Green vortex case is somewhat 337 idealised in having only periodic boundaries. To demonstrate the applicability of ESPH, subsequent test 338 cases will include basic rigid boundaries and other standard SPH numerical settings to assess ESPH in a 339 more practical context. In particular, a moderate solver tolerance of 1 × 10 −5 is now used within a first  3.92 1.00 × 10 −6 - Table 5: L 1 Error in horizontal velocity and pressure at t = 10 for different particle spacings with associated orders of convergence. Results are obtained with 6 th order Gaussian, G6.
over the domain. Note r j dθ j ≈ dx is the arclength associated with a particle volume centred at radius 347 r j . A simple analytical solution exists for the steady state azimuthal velocity given by, where ω is the angular velocity of the inner cylinder (radius R 1 = 1); the outer cylinder has a radius 349 R 2 = 2. Notably, this case allows assessment of ESPH over non-Cartesian particle distributions, as well 350 as in the presence of rigid boundaries.

351
One of the simplest SPH particle boundary conditions is implemented, so called dummy particles. Here where Δw max is the maximum difference in the velocity component between iterations over a maximum    With regard to increases in computational efficiency, Table 9 shows the CPU times required for an L 1  To demonstrate that high order convergence remains achievable for this flow and non-Cartesian geometry, 388 the analytical solution is now imposed everywhere initially to remove errors due to boundary conditions 389 and iteration to steady state. Table 10 shows the L 1 error in the calculated velocity magnitude for the 390 different kernels and particle resolutions after one iteration from steady state. Ideal spatial convergence 391 is recovered throughout.    here (as h (particle) refinement has been undertaken), but the gains in accuracy are still considerable, 415 as a close-up of Figure 12 shows (see Fig. 13). Evidently, both fourth and sixth order kernel solutions  dx Second Fourth Sixth 0.05 4.07 × 10 −3 (-) 8.02 × 10 −5 (-) 2.37 × 10 −6 (-) 0.025 1.05 × 10 −3 (1.95) 5.17 × 10 −6 (3.96) 3.74 × 10 −8 (5.99) 0.0125 2.69 × 10 −4 (1.96) 3.32 × 10 −7 (3.96) 6.01 × 10 −10 (5.96) Table 10: L 1 error in velocity for Taylor-Couette flow with exact initial and boundary conditions. Orders of convergence are given column-wise in brackets. Figure 13: A close-up of Figure 12 highlighting the discrepancy present when using second order kernel, G2.
rather than being given the actual boundary velocity (zero in this case). Here u bi is the velocity of  Table 11 shows 441 drag coefficient values calculated at Re = 15 for different kernels and particle resolutions, for the two 442 different boundary conditions. The drag coefficient, C d , is defined as where U m is the maximum velocity at the inlet and with integration over the surface of the cylinder with normal n. The error measurements in in error is now observed with decreasing particle spacing, but the error soon becomes limited at around 460 0.5% for the high-order kernels. As discussed in Section 5.2, the tolerance is sufficiently small so as to  Table 11.

473
To demonstrate that the gains in accuracy with high-order kernels remain for a range of Reynolds numbers 474 (and even for dummy boundary particles), Figure 17 shows the drag coefficients for β = 2 calculated for 475 several Reynolds numbers in the steady regime. ESPH predictions are plotted for the three Gaussian where α = 1 denotes a fully Lagrangian treatment (with α = 0 for fully Eulerian). Accordingly, ∂u ∂t | Xα is 507 Figure 19: Close up of particle distributions over transition region from Eulerian to Lagrangian frame for periodic progressive wave after one period. the simulation for one wave period from an initially regular particle distribution, Figure 19 shows the Lagrangian distribution typical of standard SPH. Furthermore, the wave profile and velocity contours at 521 t = 1T (Fig. 20) closely resemble the initial conditions (Fig. 18).

522
To further validate the ALE-ISPH approach, Figure 21 shows the horizontal velocity ( Fig. 21(a)) and 523 non-hydrostatic pressure ( Fig. 21(b)  where η is the wave elevation at a chosen horizontal location. It is a particularly sensitive flow measure, but 528 both p * and u remain in good agreement with irrotational stream-function theory (Fig. 21), despite ALE-529 ISPH including small viscous effects. Table 12   The average relative deviation of ALE-ISPH predictions for p * and u from the irrotational solution. Under either the crest or trough, the relative deviation is calculated at equispaced positions separated by 0.05 distance units along the vertical before the mean is taken.

533
Since convergence and accuracy in SPH benefits from regular particle distributions the Eulerian formu-534 lation has been thoroughly investigated. High-order spatial convergence and very good levels of accuracy 535 (or efficiency) have been demonstrated for fixed particles, while the capability of merging the Eulerian 536 form with the attractive Lagrangian features of SPH has been shown using an inviscid free-surface flow 537 test case. A detailed study into this new ALE-ISPH formulation will be undertaken for a range of free 538 surface flow test cases, and will investigate the effect of the transition region and the use of high-order ker-539 nels combined with variable particle resolution. The ALE-ISPH approach demonstrated herein is likely to 540 be of higher accuracy and more convenient than undertaking coupling of two different numerical methods

541
(one mesh based, possibly Eulerian or semi-Lagrangian, which may have good efficiency for large domains, 542 and the other Lagrangian SPH with good generality but high computational cost). Such couplings are 543 increasingly popular but are sensitive to numerical interface errors, inconsistent mathematical formula-544 tions (e.g. coupling potential and viscous flows), and are troublesome to implement for parallel processing. Even in a purely Eulerian form, the approach herein has a number of key features which may enable it 547 to compete with more established methods for many applications in the near future. The interpolative 548 nature of (E)SPH simply requires unordered summation over surrounding computation points; there is 549 no consideration of ordered node labelling (local or global), node connectivity, element/cell mappings, 550 shapes, orientations, or interface conditions. Accordingly, "mesh generation" in complex geometries has 551 the potential to be greatly simplified, requiring only a suitable distribution of unconnected particles.

552
Secondly, the simplicity of SPH (mathematically and in programming) greatly eases parallel implemen-553 tations and extension to 3D, as the discrete governing equations are essentially a series of summations, 554 with the fixed particles easily partitioned across parallel processes. Solid boundary conditions remain an 555 issue, and can restrict ideal convergence (as they do for standard 2nd order SPH). Nevertheless, even for 556 the basic (and easily implementable) dummy particle boundary conditions, there are order of magnitude 557 gains in accuracy available when using high order kernels and fixed particles. There may be options 558 available to recover higher order convergence with boundaries in practical situations in the near future: