Application of background pressure with kinematic criterion for free surface extension to suppress non-physical voids in the finite volume particle method

Lagrangian particle methods such as smoothed particle hydrodynamics (SPH) and the finite volume particle method (FVPM) can suffer from non-physical voids in the spatial discretisation, due to the inability of numerical particles to deform as continuum fluid elements can. It is known that the situation can be improved for wall-bounded flows in SPH by adding a uniform background pressure to ensure positive absolute pressure everywhere. In this article, we investigate the application of background pressure in FVPM, and show that numerical voids grow under negative pressure and collapse under positive pressure. To use this technique in free-surface flow, however, the background pressure must be applied as an atmospheric pressure at the free surface. A kinematic criterion for free surface extension (KCFSE) to differentiate physical free surfaces from new numerical voids has been developed, supplementing the inherent capability of FVPM to identify free-surface particles robustly. The novel method enables background pressure to be applied at physical free surfaces and throughout the fluid, but not in non-physical voids, facilitating the suppression of such spurious voids. The KCFSE is validated for a translating square cylinder inside a rectangular numerical domain, with and without a free surface, and liquid in an oscillating rectangular tank.


INTRODUCTION
In this article we describe an innovative method to differentiate physical free surfaces from unphysical voids in the finite volume particle method (FVPM), with the aim of improving accuracy. We show validation of the KCFSE in 2D dynamic problems.
FVPM is a numerical technique in computational fluid dynamics, introduced by Hietel et al. [1]. The basic idea is to preserve the conservative flux-based formulation of the classical finite volume method (FVM), along with the Lagrangian nature of meshless particle methods such as smoothed particle hydrodynamics (SPH), which is advantageous for moving-boundary and free-surface problems. FVPM has been extended to higher order of accuracy, [2,3] viscous flow [4], and moving-boundary problems [2,5,6]. Top-hat weight functions were introduced to eliminate quadrature error and reduce computational time [7], and these were subsequently extended to 3D [8,9,10].
Like SPH, FVPM may suffer from poor particle distribution as a result of fully Lagrangian particle motion, in which computational particles track material particle trajectories [11]. Unlike Lagrangian continuum elements, SPH and FVPM particles usually cannot change their shape to cover space, and voids may appear. In SPH, this problem is reduced by the tendency of particles to rearrange towards a uniform distribution [12,13]. The problem is more severe in higherorder schemes, since particles follow the exact trajectories more closely [11]. In addition, a void in FVPM may grow unphysically in the presence of negative pressure, as we will show in section 2. This phenomenon is related to tensile instability in SPH [14].
Dilts [15] discussed that conservative form of classical SPH formulations lead to tensile instability since the pressure gradient is not zero-order consistent (spatially constant pressure field produces acceleration), proposed use of pressure differences instead of pressure sums in the momentum equation to satisfy the consistency which has more protection against the tensile instability.
Although the proposed method in [15] is consistent, but it does not preserve the conservation. To use conservative form of SPH equations, use of particle transport methods is important to avoid unphysical acceleration which is closely related to tensile instability. Nestor et al. [2,4,5] showed that favorable particle distribution can be maintained by adding a small correction to the Lagrangian particle transport velocity, taking advantage of the Arbitrary Lagrangian-Eulerian (ALE) nature of FVPM. Lind et al. [16] suggested a particle transport correction based on Fick's diffusion law to improve the particle distribution. Colagrossi et al. [12] proposed an algorithm for particle redistribution before initialisation of flow. Sun et al. [17,18] introduced a tensile instability control (TIC) mechanism to Delta-plus-SPH scheme which can effectively prevent the initiation of the tensile instability. Khayyer and Gotoh [19] investigated the performance and stability of flows with tensile stresses with Moving Particle Semi-implicit (MPS) method. Despite the success of these methods, they do not ensure a homogeneous and isotropic particle distribution in all cases.
Another approach in SPH is to avoid unfavourable particle distribution by applying background pressure throughout the numerical domain. Positive pressure provides a favourable particle redistribution force [12,13] and avoids the tensile instability [20,21]. Morris et al. [20] proposed adding a background pressure to the field to maintain positive pressure everywhere, but also showed that excessive pressure results in a long-wavelength instability. Colagrossi et al. [12] identified a numerical error term which scales with absolute pressure (but also provides the force that drives favorable particle redistribution). Therefore, Marrone et al. [22], simulating flow over a cylinder, employed the lowest possible background pressure to avoid negative values. Violeau and Leroy [23] showed that high background pressure reduces the maximum time step for stability. Korzani et al. [24] used background pressure to avoid negative pressure in flow over a cylinder, and found that fluctuations increase with increasing value of background pressure. In summary, positive background pressure can stabilise SPH by suppressing unphysical voids, but it also introduces error, and the optimum trade-off between these effects is not clear. To date there has been no study on the application of background pressure in single phase flows with a free surface [17,24,25,26], since in the absence of surface tension the pressure at the free surface must be zero in weakly-compressible SPH [27,28].
In this article, the application of background pressure is investigated and developed for FVPM to prevent the growth of unphysical voids in flows with physical free surfaces. In the Navier Stokes equations an arbitrary constant can always be added to the pressure without changing the solution, while this requires proper treatment of background pressure at the free surface, which is facilitated in FVPM by the straightforward identification of free surface. However, it is still necessary to differentiate between physical free surfaces and non-physical ones (i.e. numerical voids). In the present research, a new kinematic criterion for free surface extension (KCFSE) is proposed for this purpose. This paper is organized as follows. Section 2 summarizes the basic FVPM formulation and the role of absolute and background pressure in the method. In section 3, the application of background pressure in free-surface flows is described, with the novel KCFSE to differentiate between physical and non-physical free surfaces. Numerical results are presented in section 4.

A. Governing Equations
The Navier-Stokes equations without source terms can be written in conservation form as where U t , x , F, G, t, and x represent the vector of conserved quantities, inviscid flux, viscous flux, time, and position, respectively. In the conserved quantities and flux functions, ρ, u, p, I, and τ are density, velocity, pressure, identity tensor, and viscous stress tensor of the fluid, respectively. The Murnaghan-Tait equation is used to calculate the pressure, where ρ 0 , c 0 , and γ are the reference density of the fluid, reference sound speed, and a constant, respectively. The value of γ is set to 7 in the present work [4].

B. Numerical Method
To derive the FVPM formulation of Eq. (1) the kernel and test function for each particle i with smoothing length h are defined as the following, respectively: where σ t , x is defined as ∑ W j t , x j . In the present work non-deformable particles are used with circular compact support of radius R c = 2h. In this paper, a top-hat weighting function W is used [7], defined as The FVPM formulation for particle i is then determined by applying the test function to the integral form of Eq. (1) for inviscid flow which is written as where N, V i , β Ij , U ij t , x ij , and dS are number of neighbours of particle i, particles volume, interparticle area vector, vector of conserved quantities at the interface of particles, velocity of the interface, and the boundary surface area vector, respectively. Introducing F(U i , U j ) to denote a numerical approximation to the inviscid numerical flux F(U ij )(t , x) - where β i b and F i b are the boundary interaction vector and numerical flux approximation, respectively. In the present work, a purely Lagrangian model x i =u i is used, so that the effect of new algorithms can be clearly distinguished from any ALE particle velocity correction. The interparticle area vector is defined as Hietel et al. [1] have shown that ∑ β ij j = 0 is a sufficient condition for interparticle interactions to enforce consistency in the FVPM. By defining the particle's volume V i t as Ψ i t , x dx, the temporal variation of the particle's volume is calculated using d dt where The AUSM + -up scheme [29] is used to compute the spatial inviscid fluxes of Eq. (3). In the current work, higher-order FVPM [29], based on the linear MUSCL finite volume approach, is used. In this approach, linear reconstruction of particle data from the barycentre to the particle interface is used to compute the fluxes. To compute the viscous stresses (Eq. (4)), the required velocity gradients are calculated by using the consistency-corrected SPH approximation [30]. To impose boundary conditions, the integral term at the RHS of Eq. (9) is used and the numerical flux over the boundary segment is prescribed and then boundary interaction vector β i b is calculated. The second-order Runge-Kutta scheme is used for temporal discretisation of Eqs. (4) and (10). The readers are referred to [4] and [7] for detailed derivation of the FVPM formulation and extension of FVPM to viscous flow, respectively.

C. Free-Surface Identification in FVPM
FVPM provides an unambiguous criterion for identification of particles that lie on a fee surface. In SPH, a particle is considered a neighbour of particle i if its centre is inside the compact support of particle i [31,32]. In FVPM, in contrast, if the compact supports of two particles overlap each other they are considered as neighbours. As illustrated in Figs. 1(a) and (b), particles j, k, m, n, and p are neighbours of i. If the boundary of a particle's support is fully covered by neighbouring particles ( Fig. 1(a)), its area vectors sum to zero (∑ β ij j = 0) [1,2]. In contrast, if ∑ β ij j ≠ 0 it means that the particle is at a free surface, adjacent to a void ( Fig. 1(b)). In practice, for fully immersed particles, ∑ β ij j = 0 is satisfied within some rounding error. Therefore, a particle is identified as a free-surface particle if it has magnitude of ∑ β ij j greater than a small tolerance.

D. Absolute, background and astmosphierc pressure in FVPM
In the following discussion, for illustration, we consider a field of particles, far from boundaries, with zero velocity at uniform pressure p. For a consistent numerical flux function F U i , U j , the right hand side of the momentum equation (10) then reduces to -∑ pβ ij j . If there are no voids between particle supports, ∑ β ij j = 0 [1,2], and hence the pressure forces acting on particle i are in balance. Fig. 2(a) indicates the pressure forces on particle i due to interacting particles k, m, n, p, and j in a region where pressure is negative, and the pairwise force -pβ ij is directed from i towards j. Thus, for any particle distribution with complete cover, FVPM ensures equilibrium independently of the absolute pressure. Fig. 2(b) shows an example with a void in the particle cover. In this case, the condition ∑ β ij j = 0 is not satisfied for any particle i adjacent to the void. In the presence of negative pressure, the resultant force on particle i acts away from the void. The imbalance of pressure forces leads to the growth of the void. Thus, unphysical voids tend to grow in the presence of negative pressure, and collapse in the presence of positive pressure. Superficially, this phenomenon resembles the tensile instability in SPH. It is also consistent with findings in SPH that addition of background pressure stabilizes the simulation, and it suggests that the same benefit may result in FVPM.
(a) Schematic view of pressure forces act on particle i in negative pressure zone while ∑ pβ ij j = 0 for particle i, (b) tracking of streamlines by fixed shape particles in fully Lagrangian motion leads to the creation of unphysical voids (blue outline). If the unphysical voids occur in negative pressure zone the act of pressure forces leads to the growth of the voids because ∑ pβ ij j ≠ 0 for particle i.

A. Background pressure at a free surface
To ensure positive pressure throughout the domain, background pressure is added to the momentum equation, so that it adds to the pressure forces between particles, as depicted in Fig. 3. A free-surface particle (like particle i in Fig.   3) has an effective surface area ∑ β ij j (directed outwards) not covered by neighbouring particles, and therefore considered to be exposed to the void. Therefore, a force ∑ β ij j is added to ensure equilibrium in the case p = 0. Thus ∑ β ij j can be interpreted as the area of the exposed free surface, and p b acts as atmospheric pressure. However, the application of background pressure force p b ∑ β ij j to all particles in this manner cannot lead to favourable results, because p b ∑ β ij j is applied not only at physical free surfaces, but also at unphysical free surfaces (exposed surface patches of particles k, m, n, o, and p in Fig. 4(a)). This would negate the purpose of background pressure to suppress the growth of voids, since particles once again experience forces away from the void wherever pressure is negative. Therefore, it is necessary to differentiate the physical free surface (exposed surface patches of particles i and j in Fig. 4(b)) from unphysical interior voids (exposed surface patches of particles k, m, n, o, and p in Fig. 4(b)), and impose background or atmospheric pressure only at the physical free surface.

B. Kinematic Criterion for Free Surface Extension (KCFSE)
In this section, we describe a novel test to distinguish between the appearance of spurious voids in particle cover and extensions of the physical free surface. Since the timestep is limited by a Courant criterion, a particle's displacement in each timestep is less than its diameter. (With typical Mach number 0.1 or less as in most free-surface applications of weakly compressible particle methods, the displacement is an order of magnitude less than particle size.) Therefore, a particle far from a free surface cannot join the free surface within one timestep. More precisely, we define a criterion that a particle can join the free surface only if at least one of its neighbours was on the free surface at the previous timestep.
The application of this novel kinematic criterion for free surface extension (KCFSE) is illustrated schematically in Fig. 5. At time t (Fig. 5 (a)), there are no interior voids, and a subset of particles are identified as free-surface particles.
At time t+∆t, due to the motion of particles, surface patches (shown in blue in Fig. 5(b)) of some particles become exposed to the void. All particles with non-zero ∑ β ij j , that do not already belong to a free surface, are candidates for addition to the free surface. Now the KCFSE is applied. Particles i and j had free-surface neighbours at time t.
Therefore, the free surface is extended to include these particles and they are labelled as such. Particles k, m, n, o, p did not have neighbours on the free surface. Therefore, their non-zero ∑ β ij j is considered to be the result of a nonphysical void, and they are not labelled as free-surface particles (Fig. 5(c)). Later, when the momentum equation is updated, the free-surface background pressure force p b ∑ β ij j (atmospheric pressure) is applied only on the labelled free-surface particles. For those particles, the net force due to background pressure is zero. Particles at the void (k,m,n,o,p), not labelled as part of the physical free surface, experience an unbalanced force due to background pressure, tending to close the void.

NUMERICAL RESULTS
The novel KCFSE algorithm has been applied for a translating square cylinder in an enclosed rectangular domain with and without a free surface, and oscillating liquid in a rectangular tank. Fully Lagrangian particle motion is used in all simulations.

A. Trasnlating Square Cylinder in an Enclosed Domain
The first test case is SPHERIC benchmark 6 [33], flow over a moving square cylinder in an enclosed rectangular domain without a free surface (Fig. 6). Here we use this test to investigate the effect of background pressure p b in FVPM. The Reynolds number based on the square's side length L and maximum velocity u s,max is 100. The ratio L ∆x ⁄ is 20, where ∆x is the initial particle spacing. The ratio h ∆x ⁄ is 0.7. The reference dynamic pressure p d and dimensionless background pressure p b * are defined as respectively, where u s (t) is the square velocity. The velocity of the square (Fig. 7) is prescribed according to the benchmark [33], with smooth acceleration from rest to constant velocity. In the following, the dimensionless flow speed u * is defined as |u| u s,max ⁄ . Results are compared against solution from an incompressible finite difference level set method [33]. are in good agreement with the benchmark (Fig. 8(a)). In the test case with zero background pressure ( Fig. 8(b)), voids appear near the top and bottom surfaces of the square cylinder and in the wake of the cylinder. This effect is removed by increasing the background pressure, as shown in Figs. 8(c-e). Qualitatively, the structure of the velocity field and the particle distribution is independent of background pressure for p b * from 4 to 200.  Fig. 8. Velocity fields at Re=100, t * =5 of (a) benchmark solution [33] at M=0.05 with L ∆x Results are shown in Fig. 9 for the pressure drag coefficient of the square cylinder, C Dp , defined as where

B. Translating Square Cylinder in a Domain with Physical Free Surface
In this test case, the upper wall of the previous case is replaced with a free surface, as shown in Fig. 10(a). The square is submerged at a depth of (7.5)L, deep enough that the effect of the free surface on the local flow is expected to be small. However, background pressure cannot be applied straightforwardly to suppress voids as in the previous fully enclosed case. Here, the new KCFSE is investigated. The Reynolds number based on the square's side length L and maximum square velocity u s,max is 100. The ratio L ∆x ⁄ is 20, where ∆x is the initial particle spacing. The ratio h ∆x ⁄ is 0.7.
Results with zero background pressure are shown in Figs. 10(b) and (c). Large voids appear in the wake, making the solution meaningless. The simulation failed shortly after this point. In Fig. 10(b) the KCFSE is not applied, and particles at the unphysical void are identified as free-surface particles, along with the true free surface. In contrast, in  Pressure fields are shown in Fig. 12 with background pressure subtracted, in the dimensionless form p * = p-p b p d .With p b acting as the effective atmospheric pressure, * is dimensionless gauge pressure. With KCFSE off (Figs. 12(a-c)), there are large unphysical voids surrounded by regions of negative gauge pressure. This is consistent with the expectation that numerical voids grow in the presence of negative pressure. With KCFSE on and positive background pressure, * ranges between -2 and +2, confirming that absolute pressure is positive everywhere, since the lowest p b * here is 4.

C. Oscillating Rectangular Tank with Physical Free Surface
In this test case, a rectangular tank contains liquid with a free surface, and oscillates vertically with displacement y where A, f, and t are amplitude, frequency, and time, respectively. Instead of moving the tank walls, the problem is solved in the reference frame of the tank, and the necessary fictitious force

D. Dam-break Problem
In this test case [34], the effectiveness of KCFSE in a problem with large deformation of free surface is assessed.
As shown in Fig. 17(a), the initial height and length of liquid are H and 2H respectively. The Reynolds number based on the liquid length and characteristic velocity gH is 1.03 × 10 6 . The ratio H ∆x ⁄ is 60, where ∆x is the initial particle spacing. The Mach number M = gH c 0 is 0.05. The ratio h ∆x ⁄ is 0.95. The dimensionless background pressure p b * is defined as p b ⁄ and dimensionless time t * is defined as t g H ⁄ .
In Fig. 17(b) the KCFSE is not applied, and particles at unphysical voids are identified as free-surface particles, along with the true free surface. In contrast, in Fig. 17(c) when the novel KCFSE is applied, only the top surface of the fluid and the void due to the wave breaking and creation of secondary wave are identified as physical free surfaces, as expected. In this test, with p b * = 1, the result confirms the capability of the KCFSE to differentiate between physical and nonphysical free surfaces.

CONCLUSIONS
The application of background pressure has been investigated in the finite volume particle method (FVPM). It has been shown that in FVPM, numerical voids (regions not covered by any particle) may grow in the presence of negative pressure, a phenomenon similar to the SPH tensile instability. It has been found that by ensuring positive absolute pressure, non-physical voids can be suppressed on FVPM, as previous researchers have shown in SPH. Furthermore, in contrast with SPH, results are independent of the magnitude of positive pressure, which therefore does not need tuning.
A new kinematic criterion for free surface extension (KCFSE) has been presented, enabling FVPM to differentiate between physical free surfaces and unphysical voids. The KCFSE is based on the principle that new isolated free surfaces should not appear in the course of a simulation, although an existing free surface may extend. This in turn allows an accurate atmospheric pressure boundary condition to apply on physical free surfaces only, which in turn enables background pressure to be applied in a flow with free surface. In principle, this approach could be implemented in any particle method where the free surface can be identified reliably.
The method has been tested for 2D flows prone to numerical void formation, with and without free surfaces, and results are in good agreement with benchmark data where appropriate.

ACKNOWLEDGEMENT
This publication has emanated from research conducted with the financial support of Science Foundation Ireland (SFI) and Aerogen ® , co-funded under the European Regional Development Fund, under Grant Number 13/RC/2073.