The Finite Volume Particle Method (FVPM) is a meshless method that in- corporates features of both Smoothed Particle Hydrodynamics and the Finite Volume Method. Here, two new formulations are presented which enhance its performance in simulation of free surface ows. One is a method for determin- ing the velocity of the free surface, making use of a partial Riemann problem to analyse the ow between the particle barycentre and the geometric free surface. The second is a well-balanced formulation for gravity forces that enables hydro- static equilibrium to be preserved exactly. Results are presented for hydrostatic, 1D impact, dambreak and deep-water standing wave test cases. Computations display convergence and good agreement with experimental data. The new methods recover exact hydrostatic equilibrium, improve robustness and reduce acoustic pressure uctuations in dambreak ows, and reduce dissipation. 1. Introduction In this article, two modi cations are proposed to the Finite Volume Par- ticle Method (FVPM) to enhance advection of a dynamic free surface and to improve accuracy in hydrostatic simulations. FVPM, like Smoothed Particle Hydrodynamics (SPH), is a meshless method derived from a particle-based ker- Corresponding author. Email address: nathan.quinlan@nuigalway.ie (Nathan J. Quinlan) Preprint submitted to Elsevier November 22, 2018 nel smoothing operation. In FVPM, the smoothing may be understood as a distribution of the volume assigned to each particle, with well-de ned sharing of volume among overlapping neighbours. The method has strong similarities with mesh-based nite volume methods, and allows arbitrary choice of the par- ticle transport velocity. It can be implemented in fully Lagrangian form, like SPH, making it a candidate for applications to ows with dynamic interfaces and free surfaces. Jahanbakhsh et al. [1] analysed FVPM for free-surface ows and applied it to jet impingement, and later to ow in Pelton turbines [2]. However, to date, FVPM has not been validated or characterised in detail for benchmarks such as hydrostatic cases and dam-break ows. In the present article, two enhance- ments are proposed for simulation of free surfaces. The rst is a novel method to determine the arbitrary Lagrangian Eulerian (ALE) transport velocity for particles on the free surface. Purely Lagrangian particle transport is a natural choice for a particle method in free-surface ow. However, the free surface is not located at the centre of free-surface particles, and thus an error may be introduced by assuming the free-surface velocity equal to the particle velocity. In FVPM this distinction is rigorous. In the proposed new method, a re ned calculation of the velocity at each particle's free surface is used as the ALE particle transport velocity. The method is explained in detail in section 3.1. The second innovation concerns the discretisation of weight (or other body forces). The exact zero-order consistency of FVPM enables a well-balanced scheme in which weight and pressure gradient can be in exact equilibrium in a hydrostatic scenario. This method is presented in section 3.2. Results of numerical tests are presented in section 4. Tests are selected to evaluate the e ects of the new development and the performance of the method as a whole in free-surface ow, with comparisons against experimental and an- alytical benchmarks. The test cases are a hydrostatic tank, one-dimensional impact, the dambreak experiment of Lobovsk y et al. [3], and a standing wave. 2 2. The Finite Volume Particle Method 2.1. Principles FVPM was introduced independently by Hietel et al. [4] and Ismagilov [5], and further analysed by Keck and Hietel [6] and Junk [7]. Here, the derivation is brie y outlined. FVPM is applied to a conservation equation of the form @U @t+rF(U) = 0; (1) where Uis the vector of conserved variables and Fis the ux. The continuum is discretised by a set of particles, each one associated with a test function i(x;t) =Wi(x(t))P jWj(x(t)): (2) The compactly supported kernel-like particle weight function Wi(x;t) is zero outside i, the support of particle i. Where a point is covered by multiple overlapping particles, i(x;t) can be interpreted as the fraction of a local volume element assigned to particle i, sinceP i i(x;t) = 1. Therefore the total volume of particleiis de ned as Vi=Z i(x;t)dx: (3) Multiplying (1) by i(x;t) and integrating, after some manipulation, the following semi-discrete equation results: d dtZ i(x;t)Udx+X j ijGij+ b iGb i= 0: (4) Here, _xiis the transport velocity of particle i, which may be chosen arbitrarily. The numerical ux function Gij=G(~Uij;~Uji;_xi;_xj) approximates F(U(xij)), as in the classical nite volume method (FVM), and also accounts for transport between particles iandjdue to relative motion. The interface states ~Uijand ~Ujiare reconstructed from the barycentres biandbjrespectively, with the linear reconstruction ~Uij=Ui+~rUi(xijbi), (where ~rdenotes a meshless approximation to the gradient) or zero-order reconstruction ~Uij=Ui. The quantity ijdenotes ij ji, where ij=Z Wi(x)rWj(x) (P kWk(x))2dx: (5) 3 Boundary terms are denoted by the superscript band b i=Z @ i(x;t)dS; (6) where Sis outward-normal boundary area. The derivation of Eq. (4) is pre- sented in detail by Hietel et al. [4]. Equation (4) is similar in structure to the mesh-based FVM, with ijand b itaking the place of cell-cell and cell-boundary interface areas. Accordingly, ijis known as the interparticle area. De nition (5) impliesP j ij= 0, which ensures zero-order consistency, and ij= ji, which ensures conservation. Cell face areas in FVM satisfy the same conditions. The rate of change of particle volume is given by di erentiation of Eq. (3) as d dtVi=X j ij_ xj ji_ xi : (7) The method also de nes the centroid of particle volume, or barycentre, bi=Z i(x;t)xdx: (8) The barycentre is not generally equal to the particle centre xi, which is used to construct the weight function Wi. Typically, particle supports are circular, spherical, square or cubic, and xiis simply the geometric centre of the support. The quantityR i(x;t)Udxin Eq. (4) is the volume-weighted average of U over the particle, and is labelled Ui. Hietel et al. [4] showed that Uiis a second-order approximation to the value of Uat the barycentre bi. Most FVPM implementations use a Godunov-type method to determine the interparticle ux G(Ui;Uj). In the basic version, the ux is derived from the solution of a Riemann problem with initial states UiandUj, i.e. zero- order reconstruction from barycentres to particle interfaces. Nestor et al. [8] introduced a higher-order scheme based on the MUSCL nite volume approach [9]. Gradients of eld variables at particle barycentres are estimated using a meshless numerical operator, and a non-linear limiter is applied. The limited gradients are used to estimate eld values at particle interfaces to second-order 4 accuracy. These reconstructed values are the initial states used to determine interparticle uxes using a Riemann solver. In principle, the reconstruction operation can be performed using any meshless approximation of the required accuracy (speci c details of the present implementation are given in section 2.5 below). The FVPM de ned and named by Hietel et al. [4] is unrelated to several similarly named methods [10, 11, 12]. FVPM in the present work is characterised uniquely by its de nitions of interparticle area (5) and particle volume (3), which are rigorously derived from the weak form of the conservation PDE. 2.2. Evaluation of interparticle area In the rst versions of FVPM, ijwas calculated by Gaussian quadrature, which is extremely expensive and an additional source of error. Quinlan et al. [13, 14] de ned Wi(x) as the top-hat function, equal to 1 inside the support of particleiand 0 elsewhere. The integral of Eq. (5) can be computed analytically, quickly and exactly. This approach was extended to 3D by Jahanbakhsh et al. for cubic [1] and spherical [2] particles. Another approach for exact calculation of ij, using splines for Wi(x) in 1D, was proposed by Kaland [15]. 2.3. Consistency and convergence A minimum of zero-order consistency is ensured, sinceP j ij= 0. With a uniform value for ux Gijin Eq. (4), the time derivative evaluates to zero, and thus a uniform eld is preserved. The method is modelled on the second-order nite volume rst-order Godunov and second-order MUSCL schemes, with con- stant or linear reconstruction respectively. However, as yet there is no complete theoretical analysis of the order of consistency or convergence of FVPM. The approximation steps that in uence spatial accuracy are the FVPM gradient operatorrFP j ijFij[6], the Godunov-type calculation of uxes, the selection of the interface points xij[6], and weighted volume-averaging (proven second-order accurate [6]). Further work is required to analyse these e ects in detail. 5 In the present work, the main variant used is the linear reconstruction from barycentre to interface, similar to the second-order MUSCL nite volume scheme. Numerical tests of FVPM with linear reconstruction have previously displayed second-order convergence in Poiseuille ow [8] and Taylor-Green ow [8, 16, 14]. The zero-order reconstruction is used in some cases for comparison. 2.4. Free surfaces in FVPM In FVPM, the free surface is de ned as all regions of all particle support boundaries @ ithat are not covered by any other particle, as shown in Fig. 1. This gives an unambiguous de nition of the mechanical free surface as a geometric surface in 3D (curve in 2D). To facilitate the following discussion, a virtual void particle vis de ned for each uid particle i, as shown in Fig. 1. (This particle does not appear in the method in practice.) It overlaps only with one uid particle i, with in nitesimal volume of overlap. The free-surface region of @ i(the region @ inot covered by any other uid particle) is fully covered by v. To maintain these conditions, vmust be advected at the same velocity as i, i.e. _xv=_xi. 2.4.1. Mass and momentum At the interface between iandv, the free-surface dynamic boundary con- dition must be enforced. There should be zero ux of momentum between the void and the uid. Similarly there should be no mass ux. Since FVPM is formulated in terms of uxes, it is straightforward to implement this. The void simply makes no contribution to the sum of uxes in Eq. (4). There is no need to enforcep= 0 at the surface particles. To do so would result in a rst-order error, since p= 0 should apply at the particle's free surface, but eld values are de ned only at the barycentre. 2.4.2. Volume Although there is no exchange of mass and momentum bewteen the uid and the void, Jahanbakhsh et al.[1] have noted that there is a conservative exchange of volume. If particle iadvances into the void, for example, it gains the volume 6 ݅ݒvoid fluidΩ݅∂Ω݅ΩݒFigure 1: Schematic diagram of a particle ion the free surface with an implicit void particle v. swept out by the exposed segment of @ i(while also exchanging volume with its neighbour particles, depending on their relative motion). The transfer of volume from vtoiis given by Eq. (7), making use of _xv=_xias follows. dVi;v dt= [ iv_ xv vi_ xi] = iv_xi (9) SinceP j ij= 0, it follows that iv=P j6=v ij. Therefore, Eq. 9 can simply be written as dVi dt=X j ij_ xj ji_ xi X j ij_xi; (10) where the particle vis discarded and summations are understood to be over all uid neighbour particles. The void particles were introduced only for the sake of this argument, and need never be implemented in practice.P j ijis interpreted as the free surface of particle i(normal from the void into the particle), just as ijis interparticle area. 7 2.4.3. Interparticle area The above approach depends on enforcement of the consistency condition P j ij= 0. In early FVPM work, ijwas approximated by quadrature of Eq. (5). Later its accuracy was improved by corrections that enforce that condition [17, 18]. These corrections were found necessary to achieve an acceptable balance between overall accuracy and computation time in ows without free surfaces. At the free surface, however,P j ij= 0 does not hold, and should not be enforced arti cially. Thus, approximate ij, based on quadrature, cannot be corrected at the free surface using available methods. Consequently, exact evaluation of ij[14] is an essential enabling technique for free-surface modelling. It ensuresP j ij= 0 is satis ed exactly in the interior of the uid, wherever there are no gaps between particle supports. At the free surface, the non-zeroP j ijrepresents the surface area of a particle which is exposed to the void. 2.5. Numerical details In the present work, the top-hat particle weight function is used to allow exact calculation of ij[14]. Gradients are calculated using SPH operators with the B-spline kernel, with the correction of Bonet and Lok [19] for rst-order consistency. The Barth-Jespersen slope limiter [20] for unstructured meshes is applied and Liou's AUSM+-up approximate Riemann solver [21] is used. Time- stepping is by the second-order Runge-Kutta scheme. Fluid is modelled as weakly compressible, using the modi ed Tait equation of state with = 7. Particle motion is purely Lagrangian throughout the numerical tests presented in this work, i.e. _x=u, unless otherwise stated. Although ALE particle velocity correction (shifting) has been employed successfully in FVPM [8, 1], in the present work it is not used, to allow other e ects to be observed more clearly. 8 3. New extensions of FVPM for free-surface ow 3.1. Enhanced free-surface particle transport velocity As discussed in section 2.4, if the free surface in FVPM is modelled without any special treatment, the mass and momentum uxes between the particle and the void are zero. This enforces the correct dynamic boundary condition at the free surface. The material velocity uof free-surface particles is updated through the conservation equations. FVPM is an ALE method, and particle transport velocity _xican be chosen arbitrarily. A natural choice for free-surface ows is Lagrangian particle motion, i.e. _xi=ui. However, uiof a free-surface particle is not strictly equal to the velocity of the free surface. The particle's nominal material velocity uiis a weighted average of material velocity over the whole of particle i, and a second-order approxima- tion to the value at the barycentre [4]. The free surface has a clear geometric de nition, consisting of the exposed parts of particles' support boundaries, as shown in Fig. 2. The velocity of the free surface is the transport velocity _x, which is arbitrarily chosen. The barycentre is separated from the free surface by a distance approximately equal to the particle radius, 2 h. Thus, in pure Lagrangian mode, when _x=u, there is a rst-order error in the velocity of the free surface. The purpose of the present work is to calculate a better approximation to the free-surface velocity, which can be applied as transport velocity _xfor free-surface particles. This is entirely distinct from enforcement of the dynamic boundary condition, which imposes no conditions directly on velocity. A secondary aim is to improve robustness, since in fully Lagrangian mode, the present FVPM method sometimes fails in simulations of violent ows in regions of rapid free- surface deformation, such as breaking waves. The approach is to solve a partial Riemann problem for the ow normal to the free surface at sub-particle scale, relating the known velocity and pressure at the barycentre to the known pressure ( p= 0) and unknown velocity at the free surface. This is based on the method of Dubois [22], as adapted by Marongiu 9      Figure 2: Schematic diagram of circular particles bounded by walls and a free surface. The free surface is marked with a heavy curve, xiis the geometric centre of particle i(shaded), biis its barycentre, gis acceleration due to gravity,P j ijis the projected area of the free- surface region of the surface of particle i, and xfsis the projection of the barycentre onto the free surface. 10   =+ d d +  =0= = =  =0 = =  Figure 3: Space-time diagram of a characteristic wave in one dimension, propagating from the barycentre of particle iatx=bito the free surface at x=xfs. et al. [23] for walls in SPH, and already employed for wall boundary conditions in the present FVPM method. The principle is illustrated in 1D in Fig. 3. The Riemann invariant u+p=(c) is constant along the wave running at velocity u+cfrom the barycentre of particle ito the free surface, leading to  u+p c i= u+p c fs: (11) (For clarity, gravity is omitted from this discussion.) Since p= 0 holds at the free surface, Eq. (11) can be solved for the free-surface velocity ufs=ui+pi ici: (12) Finally, _xiis set toufs, so that the particle is transported at the correct velocity for the free surface. In the multidimensional case, a 1D problem is solved for ow normal to the free surface. The geometry is shown in Fig. 2. The sum of the interparticle areas,P j ij, is the projected area of the particle's free surface, directed inwards normal to the surface. This is used to de ne the outward unit normal, n. A representative point on the free surface, xfs, is located by projecting bionto the surface of the particle in the direction of n. This is required to correct the Riemann invariant for hydrostatic pressure. The 1D Riemann problem is then 11 solved for the normal velocity, giving ufsn=uin+pi+ig(xxfs) ici: (13) This provides the surface-normal component of free-surface particle transport velocity, while the tangential component is set to the tangential component of the particle material velocity. The nal form of the particle transport velocity is then _xi= (ufsn)n+ui(uin)n: (14) This procedure requires identi cation of free-surface particles. In FVPM, this is unambiguous and comes at no extra cost. SinceP j ij= 0 for interior particles, free-surface particles are simply those with non-zero sum of interpar- ticle areas (allowing a tolerance for rounding errors). It should be noted that this formulation is distinct from the enforcement of the free-surface boundary condition. As discussed in section 2.4.1, the ux across the free surface is zero, independently of the velocity of the free-surface particles. 3.2. Hydrostatic balance For accurate computation of hydrostatic equilibrium, or ows where hydro- static forces overwhelm dynamics, the gravity force must be treated with care. The sum of pressure and weight force on a particle can be computed in an obvious way asZ irpdxgVi: (15) In the context of mesh-based methods, Botta et al. [24] stated the principle that for a hydrostatically well-balanced scheme, weight and pressure gradient should be discretised in the same manner so that their numerical errors cancel exactly. The expression above does not satisfy this rule, and in practice, large errors have been observed in numerical tests. Instead, in the extended method, the force per unit volume is now written asr(pgx). This leads to the FVPM discretisation of force on particle i 12 as Fi=X j(p ijigxij) ij (16) wherep ijis the interface pressure between particles iandjdetermined by a Riemann solver, and xijis thei-jparticle interface location. Local variation in density is neglected. In an exact hydrostatic pressure eld, Eq. (16) evaluates as follows. The pressurep(x) is given by gx. The present FVPM implementation uses a rst- order consistent SPH scheme to estimate gradients at the particle barycentre for use in a MUSCL-type scheme. The pressure and velocity gradients in this linear eld are therefore recovered exactly as gand zero, respectively, and the pressure and (zero) velocity are reconstructed exactly at the interface point xij, from both the iandjsides. The Riemann problem at each interface, with identical left and right states, also returns the exact pressure, p ij=pi+ig (xijbi). Substituting this into (16), the force reduces to Fi=X j[pi+ig(xijbi)igxij] ij = [piigbi]X j ij = 0; sinceP j ij= 0. Thus, equilibrium is exact in a hydrostatic pressure distribu- tion with uniform density. For a particle at the free surface, Eq. (16) is applied, using xfs(as de ned in section 3.1 and Fig. 2) as the interface point xij, andP j ijas the e ective area. The pressure p ijis zero. For clarity, the above development omits limiters. Following standard prac- tice in mesh-based nite volume methods, slope limiters are employed in FVPM to prevent the creation of new extrema in the reconstruction process. In the present work, the Barth-Jespersen slope limiter [20] is used. The limiter function 13 for particle iand a neighbour jis 'ij=8 >>>>>< >>>>>:min 1; max i i r i(xijbi) ;ifr i(xijbi)>0 min 1; min i i r i(xijbi) ;ifr i(xijbi)<0 1; ifr i(xijbi) = 0;(17) where denotes pressure or a velocity component and xij=1 2(xi+xj) is the interface point. Minima and maxima are taken over all neighbours of i. The limiter for particle i,'i, is then taken as the minimum of all pairwise limiters 'ijand is applied to the gradient of in the reconstruction from barycentre to interface points: ij= i+'ir i(xijbi); (18) wherer iis determined by a rst-order consistent SPH gradient operator in this case. Limiter functions are non-linear, and when activated, they compromise the rst-order consistency required for equilibrium. To avoid this, reconstruction and limiting operations are applied to the deviation from local hydrostatic pres- sure distribution, ^ pi(x) = (p(x)pi)ig(xxi), rather than to pressure p(x) directly. The limiter is calculated for this variable according to Eq.(17). It is then reconstructed to xijby ^pi ij= ^pi i+'i~r^pi i(xijbi) (19) and the result is converted back to absolute pressure pij=pi+ig(xxi) + ^pi ij: (20) This is used as input to the Riemann problem at the ijinterface. 4. Numerical tests 4.1. Hydrostatics 4.1.1. Test case Two geometries were used to test the hydrostatic balance model. One is a simple rectangular tank of width 2 H, lled to depth H. The second, referred 14 to as the irregular tank, has oblique surfaces. In both cases, particles are dis- tributed with  x=H = 1=10 andh=x= 0:501, and initialised with zero veloc- ity and a hydrostatic pressure distribution. In a third test, the uid in irregular tank is discretised with non-unifrom particles, ranging from  x=H = 1=10 in most of the tank to a small region of  x=H = 1=40 at the free surface, with h=x= 0:501 throughout. All results presented here were computed with the new hydrostatic balance formulation. Based on the reference velocity ( gH)1=2, the Reynolds number is 1 :1106and Mach number is 0.04. 4.1.2. Results Results are shown in Fig. 4 for the evolution of the total kinetic energy of the system. In the simple rectangular tank, kinetic energy, normalised with respect to the total gravitational potential energy, decays rapidly from 1013to 1018, on the order of machine zero. Energy is slightly higher in the irregular tank. Addition of the enhanced free-surface boundary condition has a minimal e ect on the global energy. With particles of non-uniform size near the free surface, there is higher initial kinetic energy and slower decay. This is because the barycentre is not symmetrically located in particles with neighbours of dif- ferent size, and some recon guration is necessary to reach exact equilibrium. Nonetheless, spurious velocity is at a negligible level, with dimensionless kinetic energy falling below 1013byt50. Velocity elds are shown in Fig. 5 at the instant of maximum velocity for two cases. These indicate the maximum velocity that can be expected as the system approaches numerical equilibrium, on the order of (104)pgHwith non-uniform particles. To examine convergence, simulations were carried out for a simple tank with regularly distributed particles 5 to 100 deep. Computed steady-state pressure distributions are shown in Figure 6. Results display second-order convergence to the analytical solution (21). p(y) =0c2 0 ( 1 +( 1)g c2 0(Hy) 1 1) (21) 15 0 0.5 1 1.5 2 2.5 3 3.5 4 t(g/H)1/210-1910-1810-1710-1610-1510-1410-1310-12KE/(ρgH3)regular tank irregular tank irregular tank, enhanced FS PTV(a) 0 10 20 30 40 50 t(g/H)1/210-1810-1610-1410-1210-1010-8KE/(ρgH3)regular tank irregular tank irregular tank, enhanced FS PTV irregular tank, enhanced FS PTV, non-uniform particles (b) Figure 4: Dimensionless kinetic energy in the regular and irregular hydrostatic tanks, with the hydrostatic balance formulation, with and without enhanced free-surface particle transport velocity (FS-PTV). The same data is shown on di erent scales in the two plots. 16 Figure 5: Non-dimensionalised velocity elds U=juj=pgHin the irregular tank at the instant of maximum velocity in each case. (a) Lagrangian particle transport velocity and uniform particles, t= 0:096. (b) Enhanced free-surface particle transport velocity and non- uniform particles, t= 0:596. Note that particle positions are marked at their geometric centres, and support radius is proportional to spacing ( h=x= 0:501). The tops of the free-surface particles are approximately aligned on a horizontal line, although the centres are not. 17 0 0.2 0.4 0.6 0.8 1 p/(ρ0gH)00.20.40.60.81y/HH/∆x=5 H/∆x=100 analytical(a) Pressure distribution computed for two resolutions, with the analytical solution. 10-210-1 ∆x/H10-710-610-510-4errorrms pressure error (∆x/H)2 (b) Root mean square error in non-dimensionalised pressure. Figure 6: Non-dimensionalised pressure p=(0gH) in a regular tank at t= 10, computed with the hydrostatic balance formulation. 18 4.2. Simple one-dimensional impact 4.2.1. Test case In weakly compressible methods the sound speed is a numerical parameter, usually set below the physical value. Le Touz e et al. [25] and Marrone et al. [26] investigated the role of compressibility in SPH modelling of nearly incompress- ible uids. Numerically and physically, transient impact pressure is dependent on the speed of sound. Waves resulting from impact may be dissipated by phys- ical viscosity or numerical dissipation to yield a nearly incompressible ow at longer time scales. A simple 1D impact ow has been simulated to evaluate compressibility and dissipation e ects in FVPM, preliminary to simulations of an experimental dambreak ow in section 4.3. The problem is similar to the impact of two symmetric jets considered by Marrone et al. [26]. A uid column of length L, at initial pressure 0, impacts on a wall at speed Uand corresponding Mach number M=U=c0. Theoretically, a shock wave is expected to propagate from the wall to the free surface, where it is re ected as a rarefaction. Alternating rarefaction and shock waves impact on the wall. For weak waves, which propagate at the reference speed of sound, the pressure ranges between 0 and Uc 0. Simulations were carried out with L=x= 100 and h=x= 0:48, for the FVPM variants with zero-order and linear reconstruction at M= 0:001;0:01;0:1, without the new free-surface or hydrostatic balance treatments (there is no gravity in this problem). 4.2.2. Results With zero-order reconstruction (Fig. 7), waves rapidly become smeared to an approximately sinusoidal signal, and peak pressure decays by 50% by t= tc0=L'50. With the linear reconstruction (MUSCL) scheme (Fig. 8), for M= 0:01 andM= 0:001, pressure remains essentially constant from cycle to cycle up tot= 100 and beyond, aside from a small under/overshoot which grows with time. At the highest Mach number, M= 0:1, non-linear compressibilty e ects are signi cant, leading to more complex interactions and a decaying wave. At 19 0 2 4 6 8 10 tc0/L-1-0.500.511.5p/(ρc0U)M=0.1 M=0.01 M=0.001 0 20 40 60 80 100 tc0/L-1-0.500.511.5p/(ρc0U)M=0.1 M=0.01 M=0.001Figure 7: Wall pressure in the 1D impact test, computed using FVPM with zero-order recon- struction, shown on two time scales. Pressure is non-dimensionalised by theoretical acoustic impact pressure c0Uand time by the acoustic transit time L=c 0. The same data is shown on di erent scales in both plots. (Computed without the new formulations.) the lowest Mach numbers, normalised pressure p=(Uc 0) is independent of Mach number when plotted as a function of normalised time tc0=L. This con rms that the absolute magnitude and frequency of impact pressure scale with the speed of sound or inversely with Mach number. Although there is negative pressure in these tests, there is no tensile instabil- ity. FVPM can su er a discrete analogue of SPH tensile instability, if negative pressure occurs where there is a gap between particle supports. In this case, however, the expansion is not large enough to open a void. Neighbour overlap is de ned by h=x0= 0:48 in this case, giving 2 neighbours per particle in 1D. When the uid is brought to rest by the shock wave, kinetic energy is stored by compression workR 0p2d. Energy alternates between kinetic and com- pression forms, as illustrated in Fig. 9 for FVPM with linear reconstruction at M= 0:01. Total energy should remain constant in the absence of dissipation, which may occur in strong shock waves (in cases of high Mach number) or as a result of numerical dissipation. In Fig. 10 it is evident that the total en- ergy is dissipated rapidly in the scheme with zero-order reconstruction, whereas the higher-order variant approximately conserves the sum of kinetic and com- 20 0 2 4 6 8 10 tc0/L-1-0.500.511.5p/(ρc0U)M=0.1 M=0.01 M=0.001 0 20 40 60 80 100 tc0/L-1-0.500.511.5p/(ρc0U)M=0.1 M=0.01 M=0.001Figure 8: Wall pressure in the 1D impact test, computed using FVPM with linear recon- struction, shown on two time scales. Pressure is non-dimensionalised by theoretical acoustic impact pressure c0Uand time by the acoustic transit time L=c 0. The same data is shown on di erent scales in both plots (computed without the new formulations). pression energy for M= 0:01;0:001, with losses less than 10% at t96 (24 cycles). For M= 0:1 there is much more rapid dissipation because of non-linear compressibility e ects, as expected. These benchmarks illustrate the e ects of the choice of arti cial sound speed, as discussed by previous researchers [25, 26]. They also highlight the contrast between the dissipative behaviour of FVPM with zero-order reconstruction and the very low dissipation in the higher order scheme. 4.3. Dambreak 4.3.1. Test case and experimental data This simulation models the dambreak experiments of Lobovsk y et al. [3]. As shown in Fig. 11, a volume of still water of height His initially retained behind a vertical gate, which is removed to initiate ow. The experiment was conducted 100 times. In the present work, the experimental case H= 300 mm was simulated with inviscid uid and free-slip walls, and h=x= 0:48, giving 8 neighbours per particle in the initial con guration. In most simulations, the speed of sound is set for a reference Mach number M=pgH=c 0= 0:05. 21 0 2 4 6 8 10 tc0/L00.20.40.60.81E* kinetic compression totalFigure 9: Total energy in the 1D impact test at M=0.01 in FVPM with linear reconstruc- tion, non-dimensionalised with respect to initial kinetic energy (computed without the new formulations). 0 20 40 60 80 100 tc0/L00.20.40.60.81E* M=0.1 M=0.01 M=0.001 (a) First order 0 20 40 60 80 100 tc0/L00.20.40.60.81E* M=0.1 M=0.01 M=0.001 (b) Second order Figure 10: Sum of kinetic and compression energy in the 1D impact test, non-dimensionalised with respect to initial kinetic energy (computed without the new formulations). 22 0.883 H 2HH2H 5.366 Hg p40.267 H1.713 H2.541 H h2h3h4Figure 11: Schematic diagram of the model of the dambreak experiment of Lobovsk y et al.[3], showing the location of pressure sensor p4and water height measurements h2,h3andh4. With velocity up to 2pgHoccurring in the ow, the maximum Mach number occurring is approximately 0.1. The new hydrostatic balance formulation is used in all dambreak simulations presented. Without it, signi cant errors have been observed as the tank approaches equilibrium in the long term, but e ects on the initial transient are negligible. Computational results are compared with experimental measurements of water height at three locations, and pressure at one location on the downstream wall. A video recording of one experiment [27] was analysed to record gate position as a function of time with uncertainty estimated as 5 mm. The median over all 100 experiments of the gate removal mean velocity was reported as 3.46 m/s (2pgH). An analytical function was tted to the video data, and scaled to 3.46 m/s, as shown in Fig. 12. This function was used to prescribe simulated gate motion. 4.3.2. Results Sample pressure elds are shown in Fig. 13. Pressure distribution is smooth in space. Run-up on the downstream wall and the plunging secondary wave are captured. The e ect of the enhanced free-surface particle transport velocity formula- tion is illustrated by Fig. 14, showing wall pressure p4. In comparison with a model based on a simple Lagrangian free surface, the new boundary condition 23 0 20 40 60 80 100 time (ms)050100150200250300350gate position (mm)Data (Lobovský et al., 2014) Fit Fit, scaled to 3.46 m/sFigure 12: Data from one experiment of Lobovsk y et al.[3] for gate removal, with the tted function, and the same function scaled to give the reported median gate velocity. results in lower pressure oscillations. It also provides a more robust simulation. Without any additional treatment such as particle velocity correction (shift- ing), the simulation with a simple Lagrangian free surface fails at t4:5 due to locally poor particle distributions. The motivation for the new method is to provide correct motion for the free surface, rather than the geometric centre, of each particle. Wall pressure p4, computed with various resolutions, is compared quanti- tatively with experimental data [3] in Fig. 15. High-frequency oscillation is evident in the computational results at low resolution. The amplitude of this oscillation decreases rapidly as  xis decreased and the pressure converges. Numerical experiments have been conducted to investigate the relationship between high-frequency oscillations in the dambreak and the acoustic waves discussed in section 4.2 above. Results for various Mach number in Fig. 16 show that amplitude of oscillations increases as Mach number is decreased. Results of a computation with zero-order reconstruction are shown in Fig. 17. In comparison with Fig. 15, oscillations are greatly reduced by the increased dissipation in the lower-order scheme. Together, these results con rm that the oscillations observed in Fig. 15 are due to acoustic waves, which are preserved as a result of the inherent low dissipation of FVPM with linear reconstruction. To suppress unwanted acoustic noise in weakly compressible SPH, Kiara et 24 (a)t= 0:43, gate lifting. (b)t= 6:15, gate removed. (c)t= 7:15, gate removed. Figure 13: Dimensionless pressure p=p=(gH) in the dambreak ow with H=x= 120, M= 0:05, linear reconstruction, and new formulations for hydrostatic balance and enhanced free surface particle transport velocity. Note change of scale in (c). 25 2 3 4 5 6 7 t(g/H)1/2-0.200.20.40.60.811.21.4p/(ρgH)Lagrangian FS enhanced FS-PTVFigure 14: Computed pressure p4on the downstream wall for H=x= 60 andM= 0:05 and linear reconstruction, with and without the enhanced free surface particle transport velocity, with the new formulation for hydrostatic balance. 2 3 4 5 6 7 t(g/H)1/2-0.200.20.40.60.811.21.4p/(ρgH)H/∆x=15 H/∆x=30 H/∆x=60 H/∆x=120 experiment (Lobovsky et al.) Figure 15: Computed pressure p4on the downstream wall for various H=x, with Mach num- berM=pgH=c 0= 0:05, linear reconstruction, and enhanced free surface particle transport velocity. Also shown are the 2.5, 50 and 97.5 percentiles of the experimental data [3]. 26 2.5 3 3.5 4 4.5 500.20.40.60.8p* M=0.0500 2.5 3 3.5 4 4.5 500.20.40.60.8p* M=0.0250 2.5 3 3.5 4 4.5 5 t*00.20.40.60.8p* M=0.0125Figure 16: Computed pressure p4on the downstream wall for various Mach number M= pgH=c 0, withH=x= 60, linear reconstruction, and the new formulations for free surface particle transport velocity and hydrostatic balance. 2 3 4 5 6 7 t(g/H)1/2-0.200.20.40.60.811.21.4p/(ρgH)H/∆x=30 H/∆x=60 H/∆x=120 experiment (Lobovsky et al.) Figure 17: Computed pressure p4on the downstream wall for various H=x, with Mach numberM=pgH=c 0= 0:05, zero-order reconstruction, and the new formulations for free surface particle transport velocity and hydrostatic balance. Also shown are the 2.5, 50 and 97.5 percentiles of the experimental data [3]. 27 al. [28] and Meringolo et al. [29] have proposed ltering approaches, based on the observation that the acoustic signal can often be distinguished from the underlying incompressible ow [25]. Here, a relatively crude post-processing approach is taken to smooth the FVPM pressure data of Fig. 15. The charac- teristic frequency for acoustic signals in the initial water column of height His c0=(4H). A second-order low-pass Butterworth lter is applied with cut-o at half this frequency, i.e. c0=(8H). Filtered data in Fig. 18, along with the raw data in Fig. 15, allow easier examination of convergence and comparison with experiment. Before plunging of the secondary wave at t6:6, the underlying low- frequency pressure signal is converged for all but the coarsest particles, H=x 30. For the pressure peak at t6:6, there is good convergence at H=x= 60, and reasonable agreement on the magnitude of the peak with the higher end of the range of experimental pressure. In comparison with experiment, the pressure levels are generally lower. Similar reduced pressure was observed by other authors simulating the same experiment with single-phase SPH [30] and two-phase particle nite element method [31]. The peak pressure is delayed with respect to the experimental results; this has also been observed in other numerical simulations and has been shown to be due to the absence of the air phase [32, 30]. All dambreak results presented in Figs. 13 to 18 were obtained with the dynamic gate velocity function shown in Fig. 12 (scaled to 3.46 m/s mean velocity) prescribed for the removal of the gate. To investigate the sensitivity of results to the gate removal model, two other models for gate removal were investigated. The no-gate case represents instantaneous gate removal with an incompressible uid. The instantaneous pressure eld at t= 0 was determined by a nite-di erence solution of the Laplace equation, with boundary conditions of zero pressure at the free surface and zero normal pressure gradient at the walls, following the approach of Stansby et al. [33]. The unsteady simulation then proceeded as usual with the weakly compressible model. In the third case, the gate was started impulsively and lifted at constant velocity of 3.46 m/s. 28 2 3 4 5 6 7 t(g/H)1/2-0.200.20.40.60.811.21.4p/(ρgH)H/∆x=15 H/∆x=30 H/∆x=60 H/∆x=120 experiment (Lobovsky et al.)Figure 18: Computed pressure p4on the downstream wall for Mach number M=pgH=c 0= 0:05, withH=x= 60, linear reconstruction, new formulations for free surface particle trans- port velocity and hydrostatic balance, and low-pass lter applied at c0=(8H) (14.2 Hz). See Fig. 15 for raw data. Also shown are the 2.5, 50 and 97.5 percentiles of the experimental data [3]. Results are shown in Fig. 19. In all three cases, as in the experiments, the pressure rise on the downstream wall occurs in two stages, both approximately linear in time, with a much higher rate in the rst stage. With either a fully dynamic or constant-velocity gate, the rst impact is later than without a gate, and pressure reaches a higher level in the initial rapid stage. This behaviour is closer to the experimental results. These moving-gate results suggest that modelling the nite gate removal time is necessary for accurate prediction of the early transient stages of ow far downstream, but the results are not sensitive to the detailed kinematics of the gate. Results for water height are shown in Fig. 20, with experimental measure- ments from an arbitrarily selected run. Transient peaks in water height for the ner resolutions are due to particles splashing high above the bulk of liquid. Apparent noise at low resolution is because the free surface in FVPM consists of circular arc segments of the particles' support boundaries (see Fig. 2), and the numerical water height sensor detects the highest point of the particles cov- 29 2.5 3 3.5 t(g/H)1/2-0.100.10.20.30.40.5p/(ρgH)no gate constant-speed gate dynamic gate experiment (Lobovsky et al.)Figure 19: Computed pressure p4on the downstream wall for Mach number M=pgH=c 0= 0:05, withH=x= 60, linear reconstruction, and the new formulations for free surface particle transport velocity and hydrostatic balance, for no gate (instantaneous gate removal), gate removal at constant speed of 3.46 m/s, and a dynamic gate as shown in Fig. 12. Also shown are the 2.5, 50 and 97.5 percentiles of the experimental data [3]. ering itsxlocation. Before impact on the downstream wall, FVPM results agree closely with the experiment for all resolutions, even H=x= 15. After break- ing of the secondary wave ( t6:6) there is satisfactory convergence at h4(the transducer closest to the downstream wall). As the re ected waves propagate leftwards to h3andh2, with more complex and violent dynamics observed in the experiments, results become increasingly sensitive to resolution. 4.4. Standing wave 4.4.1. Test case Results are presented here for a standing wave in an inviscid uid, following Antuono et al.[34]. A potential solution gives the velocity eld as u(x;y;t ) =r Hg 2!cosh[k(y+H)] cosh(kH)cos(kx) cos(!t) : (22) For an inviscid uid, the wave should oscillate without loss of energy. The circular frequency !is [gktanh(kH)]1=2,His depth, and = 2A=H where A is 30 0 2 4 6 8 10 t(g/H)1/200.20.40.60.81h2/HH/∆x=15 H/∆x=30 H/∆x=60 H/∆x=120 experiment (Lobovsky et al.) 0 2 4 6 8 10 t(g/H)1/200.20.40.60.81h3/HH/∆x=15 H/∆x=30 H/∆x=60 H/∆x=120 experiment (Lobovsky et al.) 0 2 4 6 8 10 t(g/H)1/200.20.40.60.81h4/HH/∆x=15 H/∆x=30 H/∆x=60 H/∆x=120 experiment (Lobovsky et al.)Figure 20: Computed water depth in the dambreak problem for various H=x, compuoted with the new formulations for free surface particle transport velocity and hydrostatic balance, along with data from a single experiment [3], 31 amplitude. The domain has periodic boundaries at x= 0;Hand a wall at y= H. The uid is initialised with a level free surface at y= 0 and the analytical velocity eld (22) evaluated at t= 0, and hydrostatic pressure. The amplitude parameteris 0.05 and H=xranges from 20 to 160, with h=x= 0:48 in all cases. The new hydrostatic balance formulation is applied, and results are presented with and without the new free-surface particle transport formulation. The test case is used here to assess the dissipative behaviour of the method. In the viscous case, the kinetic energy decays as Ek(t) =2gH3 32e4k2t[1 + cos(2!t)]; (23) whereis kinematic viscosity [35]. By tting numerical results to this function, an e ective viscosity can be calculated. This in turn is used to determine an e ective Reynolds number based on Hand the characteristic velocity Hgk= (2!) of Eq. (22), which quanti es the numerical dissipation of the nominally inviscid solution. 4.4.2. Results Sample particle pressure distributions at the fourth peak of the wave are shown in Fig. 21. In this view there is no discernible di erence between results for the pure Lagrangian mode and enhanced free-surface particle transport ve- locity. Time histories of kinetic energy are shown in Fig. 22. Dissipation is strongly dependent on resolution, with e ective Reynolds number ranging from 31 for H=x= 20 to 959 for H=x= 160 with fully Lagrangian particle transport. With enhanced free-surface particle transport velocity, dissipation is slightly less, withRe= 1253 at the nest resolution. 5. Conclusions Two new developments have been presented that enhance FVPM simulation of free-surface ows. The rst is a new method to determine the ALE particle transport velocity for particles at the free surface, based on solution of a partial 32 (a) (b) Figure 21: Dimensionless pressure p=p0g(Hy) 1 20U2 in the standing wave test at t=t=T= 4:25 withH=x= 160 (a) with fully Lagrangian particle motion and (b) with enhanced free-surface particle transport velocity. 33 0 2 4 6 8 10 t/T00.20.40.60.81KE/KE0 Re*=31Re*=90Re*=315Re*=959L/∆x=20 L/∆x=40 L/∆x=80 L/∆x=160(a) 0 2 4 6 8 10 t/T00.20.40.60.81KE/KE0 Re*=31Re*=95Re*=348Re*=1253L/∆x=20 L/∆x=40 L/∆x=80 L/∆x=160 (b) Figure 22: Kinetic energy with tted decay curves e ective Reynolds number Rein the standing wave test for various H=x(a) with fully Lagrangian particle motion and (b) with enhanced free-surface particle transport velocity.34 Riemann problem to determine the free-surface velocity consistent with zero pressure at the surface. These may be considered a form of ALE particle ve- locity correction (shifting), speci cally adapted for the free surface, and based on an analysis of the local ow rather than a measure of uniformity or sym- metry of the particle distribution. Aside from this treatment, all simulations presented here have used fully Lagrangian particle motion. The new free-surface method has been shown to improve robustness and reduce high-frequency pres- sure oscillations in dambreak simulations. Dambreak results show convergence and agreement against experimental data from the literature. In particular, the simulations show that the nite time taken for gate removal has a signi cant e ect on the impact ow downstream. In simulations of a standing wave, the enhanced free-surface transport velocity formulation slightly reduces dissipation. The second enhancement is a well-balanced formulation which enables exact equilibrium between gravity and pressure gradient. In hydrostatic tests with this formulation, normalised kinetic energy tends to machine zero. This per- formance is maintained with irregular geometry, and the enhanced free-surface particle transport velocity formulation, and non-uniform particles. The hydro- static pressure distribution approaches the analytical solution with second-order convergence. Furthermore, the dissipation and compressibility properties of FVPM have been assessed in numerical experiments on one-dimensional impact tests. FVPM with linear reconstruction shows very low dissipation, allowing acoustic waves to be preserved through many re ection cycles. Consequently, high-frequency pressure oscillations are evident in dambreak simulation, particularly with low spatial resolution. Oscillations can be suppressed by reverting to the scheme with zero-order reconstruction (introducing strong arti cial dissipation), or as a preliminary test has shown, by ltering. The enhanced free-surface par- ticle transport velocity formulation reduces high-frequency oscillations in the dambreak test. Results of standing wave simulations suggest that the new ap- proach reduces dissipation, and its e ect of reducing pressure uctuations in the dambreak is therefore associated with better tracking of the free surface, 35 as intended. Collectively, this work improves the accuracy and robustness of FVPM for free surface ows.