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 discr etisation, due to the inability of numerical particles to deform as continuum fluid elements can. It is known th at the situation can be improved for wall-bounded flows in SPH by adding a uniform background pressu re to ensure positive absolute pressure everywhere. In this article, we investigate the application of background pressure in FV PM, and show that numerical voids grow under negative pressure and collapse under positive pressure. To use this t echnique in free-surface flow, however, the background pressure must be applied as an atmospheric pressure at th e 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 partic les robustly. The novel method enables background pressure to be applied at physical free surfa ces 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. Key Words: Finite volume particle method, m eshless method, particle method, free surface 1. INTRODUCTION In this article we describe an innovative method to differ entiate physical free surfaces from unphysical voids in the finite volume particle method (FVPM), with the aim of im proving accuracy. We show validation of the KCFSE in 2D dynamic problems. FVPM is a numerical technique in comput ational fluid dynamics, introduced by Hietel et al. [1]. The basic idea is to preserve the conservative flux-base d formulation of the classical finite volume method (FVM), along with the Lagrangian nature of meshless part icle methods such as smoothed partic le hydrodynamics (SPH), which is advantageous for moving-boundary and free- surface problems. FVPM has been ex tended to higher order of accuracy, [2,3] viscous flow [4], and moving-boundary problems [2,5,6]. Top-hat weight functions we re introduced to eliminate quadrature error and reduce computational time [7], and these were subsequently ex tended to 3D [8,9,10]. M. H. Moghimi & N. J. Quinlan Engineering Anal ysis with Boundary Elements, December 2018 Page 2 of 21 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 La grangian continuum elements, SPH and FVPM particles usually cannot change their shape to cover spa ce, and voids may appear. In SP H, this problem is reduced by the tendency of particles to rearrange towards a uniform distribution [12,13]. The problem is more severe in higher- order 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 w ill show in section 2. This phe nomenon 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 e quation to satisfy the consiste ncy which has more protection against the tensile instability. Although the proposed method in [15] is consistent, but it does not preserve th e conservation. To use conservative form of SPH equations, use of particle transport methods is importa nt to avoid unphysical acce leration 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 part icle transport velocity, ta king advantage of the Arb itrary Lagrangian-Eulerian (ALE) nature of FVPM. Lind et al. [16] suggested a particle transport co rrection 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 c ontrol (TIC) mechanism to Delta-plus-SPH scheme which can effectively prevent th e initiation of the tensile instability. Khayyer and Gotoh [19] investigated the performance and stability of flows with te nsile stresses with Moving Particle Semi-implicit (MPS) method. Despite the success of these methods, they do not ensure a homogene ous and isotropic particle distribution in all cases. Another approach in SPH is to avoi d 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-wavelengt h instability. Colagrossi et al. [12] identified a numerical error term which scales with absolute pressure (but also provides th e force that drives favorable particle redistribution) . Therefore, Marrone et al. [22], simulating flow over a cyli nder, employed the lowest possible background pressure to avoid negative va lues. 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 fluctuatio ns 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]. M. H. Moghimi & N. J. Quinlan Engineering Anal ysis with Boundary Elements, December 2018 Page 3 of 21 In this article, th e 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 equa tions an arbitrary constant can always be added to the pressure without changing the solution, wh ile this requires proper treatme nt of background pressure at the free surface, which is facilitated in FVPM by the straightforward iden tification of free surface. However, it is still necessary to differentiate between physi cal free surfaces and non-physical ones (i.e. numerical voids). In the present research, a new kinematic criteri on 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 ro le of absolute and background pressure in the method. In section 3, the app lication of background pressure in free-surface flows is described, with the novel KCFSE to differentiate betwee n physical and non-physical free surfaces. Numerical results are presented in section 4. 2. THE FINITE VOLUME PARTICLE METHOD A. Governing Equations The Navier-Stokes equations without source te rms can be written in conservation form as ∂Uሺt , xሻ ∂t+׏ . ൫F ൫Uሺt , xሻ൯-G൯=0 , (1) U = ቄρ ρuቅ , (2) F = ቄρu ρuu+pIቅ , (3) G = ቄ0 τቅ , (4) 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 conserve d quantities and flux functions, ρ, u, p, I, and τ are density, velocity, pressure, identity tensor, and viscous stre ss tensor of the fluid, respec tively. The Murnaghan-Tait equation p= ρ0c02 γ൤൬ρ ρ0൰γ -1൨ (5) is used to calculate the pressure, where ρ0, c0, 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) th e kernel and test function for each particle i with smoothing length h are defined as the following, respectively: Wiሺt , xሻ =൜Wሺx - xiሺtሻ , h ሻ x∈Ωi , 0 x∉Ωi , (6) M. H. Moghimi & N. J. Quinlan Engineering Anal ysis with Boundary Elements, December 2018 Page 4 of 21 Ψiሺt , xሻ =ቐWiሺt , xሻ σሺt , xሻ x∈Ωi , 0 x∉Ωi , (7) where σሺt , xሻ is defined as ∑Wjሺt , xሻj . In the present work non-deformable par ticles are used with circular compact support of radius Rc = 2h. In this paper, a to p-hat weighting function W is used [7], defined as Wiሺt , xሻ =൜1 x∈Ωi 0 x∉Ωi . (8) 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 d dt൫UiሺtሻViሺtሻ൯ = -෍βij . ൣF ൫Uij൯ሺt , xሻ - Uijሺt , xሻ xሶijሺtሻ൧ j-නΨiሺt , xሻF ൫Uij൯ሺt , xሻ . dS , (9) where N, Vi, βIj, Uijሺtሻ, xሶij, and dS are number of neighbours of particle i, particles volume, interparticle area vector, vector of conserved quantities at the inte rface of particles, velocity of the inte rface, and the boundary surface area vector, respectively. Introducing F(Ui , Uj) to denote a numerical approximati on to the inviscid numerical flux F(Uij)(t , x) - Uij(t , x) xሶij(t), Eq. (6) is written as d dt൫UiሺtሻViሺtሻ൯ = -∑βij . ൣF ൫Ui , Uj൯൧ j -βib.Fib , (10) where βib and Fib are the boundary interaction vector and numerical flux approximation, respectively. In the present work, a purely Lagrangian model xሶi=ui 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 βij = ׬Ψiሺt , xሻ׏ Wjሺt , xሻ σሺt , xሻdx - ׬Ψjሺt , xሻ׏ Wiሺt , xሻ σሺt , xሻdx . (11) Hietel et al. [1] have shown that ∑βij j= 0 is a sufficient condition for interpartic le interactions to enforce consistency in the FVPM. By defini ng the particle’s volume Viሺtሻ as ׬Ψiሺt , xሻdx, the temporal variation of the particle’s volume is calculated using d dtViሺtሻ =෍൬γij . xሶjሺtሻ - γji . xሶiሺtሻ൰ j. (12) where γij is ׬Ψiሺt , xሻ׏ Wjሺt , xሻ σሺt , xሻdx . The AUSM+-up scheme [29] is used to compute the spatial invisc id fluxes of Eq. (3). In the current work, higher-order FVPM [29], based on the linear MUSCL fi nite volume approach, is used. In th is 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 M. H. Moghimi & N. J. Quinlan Engineering Anal ysis with Boundary Elements, December 2018 Page 5 of 21 segment is prescribed and th en boundary interaction vector βib is calculated. The second- order Runge-Kutta scheme is used for temporal discretisation of Eqs. (4) and (10). The re aders 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 voi d (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. (a) (b) Fig. 1. Summation of interparticle area vectors between particle i and its neighbours is (a) zero and (b) non-zero. D. Absolute, background and astmosphierc pressure in FVPM In the following discussion, for illustration, we consider a field of particle s, far from boundaries, with zero velocity at uniform pressure p. For a consistent numerical flux function F ൫Ui , Uj൯, 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) indicat es the pressure forc es 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 en sures equilibrium independently of the absolute pressure. Fig. 2(b) shows an example with a void in th e particle cover. In this case, the condition ∑βij j= 0 is not satisfied for any M. H. Moghimi & N. J. Quinlan Engineering Anal ysis with Boundary Elements, December 2018 Page 6 of 21 particle i adjacent to the void. In the presence of nega tive 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 pres ence of positive pressure. S uperficially, this phenomenon resembles the tensile instability in SPH. It is also consistent with findings in SPH that addition of background pressure stabilizes the simulation, a nd it suggests that the same benefit may result in FVPM. (a) (b) Fig. 2. (a) Schematic view of pre ssure 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 for ces leads to the growth of the voids because ∑pβij j≠ 0 for particle i. 3. NEW FREE -SURFACE EXTENSION METHOD AND VOID SUPPRESSION A. Background pressure at a free surface To ensure positive pressure throughout the domain, backgr ound pressure is added to the momentum equation, so that it adds to the pressure forces between particles, as de picted in Fig. 3. A free-surface particle (like particle i in Fig. 3) has an effective surface area െ∑βijj (directed outwards) not covered by neighbouring particles, and therefore considered to be exposed to the void. Therefore, a force ݌௕∑βijj is added to ensure equilibrium in the case p = 0. Thus െ∑βijj can be interpreted as the area of the exposed free surface, and pb acts as atmospheric pressure. M. H. Moghimi & N. J. Quinlan Engineering Anal ysis with Boundary Elements, December 2018 Page 7 of 21 Fig. 3. Schematic view of the forces due to uniform pressure act on the free-surface particle i in FVPM with zero velocity along imposing of background pressure. However, the application of background pressure force pb∑βijj to all particles in this manner cannot lead to favourable results, because pb∑βijj 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 differen tiate the physical free surface (exposed surface patches of particles i and j in Fig. 4(b)) from unphysical interior vo ids (exposed surface patches of particles k, m, n, o, and p in Fig. 4(b)), and impose backgrou nd or atmospheric pressure only at the physical free surface. (a) (b) Fig. 4. Schematic view of im posing background pressure pb on free-surface particles (red particles) and suppression of unphysical void: (a) an unphysical void (blue lines created by particles k, m, n, o, and p) is inside the domain of flow and the background pressure is imposed on it incorrectly, (b) imposing of background pressure along us ing the KCFSE prevents the creation of unphysical void. B. Kinematic Criterion for Free Surface Extension (KCFSE) In this section, we describe a novel test to distinguish be tween the appearance of spurious voids in particle cover and extensions of the physical free surface. Since the time step is limited by a Courant criterion, a particle’s M. H. Moghimi & N. J. Quinlan Engineering Anal ysis with Boundary Elements, December 2018 Page 8 of 21 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 (show n in blue in Fig. 5(b)) of some particles become exposed to the void. All particles with non-zero ∑βijj, 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 th ese 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 ∑βijj is considered to be the result of a non- physical void, and they are not labelle d as free-surface particles (Fig. 5(c)). (a) (b) (c) Fig. 5. Schematic view of novel kinematic criterion for free su rface extension (KCFSE): (a) particles at time t (red particles are free-surface particles); (b) updated part icle distribution at time t+∆t (candidate free-surface patc hes shown in blue); (c) after implementing the KCFSE the physical free-surface particles (red) are differentiated from the internal void. Later, when the momentum equation is update d, the free-surface background pressure force pb∑βijj (atmospheric pressure) is applied only on the labelled free-surface particl es. 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. 4. NUMERICAL RESULTS The novel KCFSE algorithm has been app lied for a translating square cylinder in an enclosed rectangular domain with and without a free surface, and oscillating liquid in a re ctangular tank. Fully Lagrangian particle motion is used in all simulations. M. H. Moghimi & N. J. Quinlan Engineering Anal ysis with Boundary Elements, December 2018 Page 9 of 21 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 e ffect of background pressure pb in FVPM. The Reynolds number base d on the square’s side length L and maximum velocity us,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 pd and dimensionless background pressure pb* are defined as 1 2ρus,max2 and pbpd⁄ respectively. The Mach number M is defined as us,max c0⁄. The dimensionless time t* and non-dimensional square velocity us*ሺݐሻ are defined as tus,max L and usሺtሻ us,max respectively, where us(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|us,max⁄ . Results are compared against solution from an incompressible finite difference level set method [33]. Fig. 6. Translating square cylinder in side a rectangular numerical domain. Fig. 7. Normalized velocity us* as a function of t*. The velocity fields of the FVPM solutions for various dimensionless background pressure s, shown in Figs. 8(b–e), are in good agreement with the benchmark (F ig. 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 sh own in Figs. 8(c-e). Qualitatively, the st ructure of the velocity field and the particle distribution is indepe ndent of background pressure for pb* from 4 to 200. (a) t*us* 0 0.5 1 1.5 2 2.500.20.40.60.811.2 M. H. Moghimi & N. J. Quinlan Engineering Anal ysis with Boundary Elements, December 2018 Page 10 of 21 (b) pb*=0 (c) pb*=4 (d) pb*=10 (e) pb*=200 Fig. 8. Velocity fields at Re=100, t*=5 of (a) benchmark solution [33] at M=0.05 with L∆x⁄ =20 for (b) pb*=0, (c) pb*=4, (d) pb*=10, and (e) pb*=200. Results are shown in Fig. 9 for the pressure drag coefficient of the square cylinder, CDp, defined as CDp = Fxp 1 2ρus2L , (14) where Fxp is the x-component of the pressure force on the cylinder. For all background pressure from pb*ൌ0 to pb*ൌ 200, results agree with the benchmark solution in the accele rating phase of the motion. FVPM results fluctuate in the constant-velocity phase due to acoustic wa ves in the weakly compress ible solution, in contrast with the incompressible benchmark solution. Results for pb* ൌ 0 differ significantly from all cases of pb* ൐ 0, and display very large oscillations around ݐ∗ൌ 6. For pb* 10 and 200, the FVPM results are not sensitive to the magnitude of the background pressure. At lower Mach number, the amplitude of oscillation in drag is lo wer, confirming that oscillati ons are due to compressibility effects. M. H. Moghimi & N. J. Quinlan Engineering Anal ysis with Boundary Elements, December 2018 Page 11 of 21 (a) M=0.05 (b) M=0.05, zoom view of (a) (c) M=0.02 (d) M=0.02, zoom view of (c) Fig. 9. Pressure drag coefficient of a transla ting square cylinder inside a rectangular numerical domain for various background pressu res and for (a) M=0.05, (b) zoom view of (a) for t*=5 to 7, (c) M=0.02, and (d) zoom view of (c) for t*=5 to 7. B. Translating Square Cylinder in a Do main 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 fr ee surface on the local flow is expected to be small. However, background pressu re cannot be applied straightforwardl y 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 us,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 th is point. In Fig. 10(b) the KCFSE is not applied, and particles at the unphysical void are identified as free-surface pa rticles, along with the true free surface. In contrast, in Fig. 10(c) when the novel KCFSE is applied, only the top su rface of the fluid is identified as a physical free surface, t*CDp 0 2 4 6 80246810 pb*=0 pb*=10 pb*=200 Benchmark[33] t*CDp 5 5.5 6 6.5 701234 t*CDp 0 2 4 6 80246810 pb*=0 pb*=10 pb*=200 Benchmark[33] t*CDp 3 3.5 4 4.5 511.522.53 M. H. Moghimi & N. J. Quinlan Engineering Anal ysis with Boundary Elements, December 2018 Page 12 of 21 as expected. In this test, with pb* ൌ 0, identification of the free surface has no consequences, but the result confirms the capability of the KCFSE to differentiate between physical and nonphysical free surfaces. (a) (b) (c) Fig. 10. (a) Translating square cylinder inside a rectangular numerical domain with physical free surface. Particle distribution with L∆x⁄ =20 at pb*=0, Re=100, M=0.05, t*=4.4, (b) without KCFSE and (c) with KCFSE. Id entified free-surface particles are coloured red. Fig. 11 shows the velocity fields for various positiv e background pressures. Without KCFSE, the simulations developed large voids, and in some cases failed. However, with KCFSE activated and pb* ൒ 4, background pressure leads to suppression of unphysical voids (Figs. 11(d-f)), and the velocity fiel d is practically independent of the magnitude of the background pressure. The results with KCFSE are in good agreement w ith benchmark [33] shown in Fig. 8(a). (a) pb*=4, KCFSE off, t*ൌ૝ .ૠ (b) pb*=10, KCFSE off, t*ൌ૝ .ૠ (c) pb*=200, KCFSE off, t*ൌ૞ M. H. Moghimi & N. J. Quinlan Engineering Anal ysis with Boundary Elements, December 2018 Page 13 of 21 (d) pb*=4, KCFSE on, t*ൌ૞ (e) pb*=10, KCFSE on, t*ൌ૞ (f) pb*=200, KCFSE on, t*ൌ૞ Fig. 11. Velocity fields at Re=100, M=0.05, with L∆x⁄ =20, with and without KCFSE. Where t*൏5, simulation failed shortly after the time shown. Pressure fields are shown in Fig. 12 with backgr ound pressure subtracted, in the dimensionless form p*=൫p-pb൯pdൗ.With pb acting as the effective atmospheric pressure, ݌∗ is dimensionless gauge pressure. With KCFSE off (Figs. 12(a-c)), there are large unphysical voids su rrounded by regions of negativ e 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 pb* here is 4. (a) pb*=4, KCFSE off, t*ൌ૝ .ૠ (b) pb*=10, KCFSE off, t*ൌ૝ .ૠ (c) pb*=200, KCFSE off, t*ൌ૞ M. H. Moghimi & N. J. Quinlan Engineering Anal ysis with Boundary Elements, December 2018 Page 14 of 21 (d) pb*=4, KCFSE on, t*ൌ૞ (e) pb*=10, KCFSE on, t*ൌ૞ (f) pb*=200, KCFSE on, t*ൌ૞ Fig. 12. Pressure distributions at Re=100, M=0.05, with L∆x⁄ =20 with and without KCFSE. Where t*൏5, simulation failed shortly after the time shown. (Note the discontinuity in the colour map at p*=0.) Fig. 13 shows the pressure drag coefficient CDp of the square cylinder for cases with KCFSE. Results are similar to the FVPM results for an enclosed dom ain, and agree with the benchmark [33] to the same extent, with fluctuations in the constant-velocity phase. The solution for zero background pressure fails around t*=4.4 and t*=4.5 (depending on Mach number) as shown in Figs. 13(a) and (c) respectively. For pb*൐0, the results are practically independent of the background pressure. (a) M=0.05 (b) M=0.05, zoom view of (a) t*CDp 0 2 4 6 80246810 pb*=0 pb*=10 pb*=200 Benchmark[33] t*CDp 2 2.5 3 3.5 40.511.522.53 M. H. Moghimi & N. J. Quinlan Engineering Anal ysis with Boundary Elements, December 2018 Page 15 of 21 (c) M=0.02 (d) M=0.02, zoom view of (c) Fig. 13. Pressure drag coefficients at Re=100, with L∆x⁄ =20, with KCFSE. (a) M=0.05, simulation failed for pb*=0, (b) zoom view of (a), (c) M=0.02, simulation failed for pb*=0, and (d) zoom view of (c). The effect of the KCFSE on the calcu lated pressure drag coefficient for pb*=10 and pb*=200 is shown in Fig. 14. In the cases in which the KCFSE is not used, three simulations fail ( M=0.02 for pb*=10 and pb*=200 and M=0.05 for pb*=10) and one solution ( M=0.05 for pb*=200) continues with very large erro rs. For the cases where KCFSE is used, the pressure drag coefficients are in good agreement with benchmark solution [33] and i ndependent of the magnitude of the background pressure. (a) M=0.05 (b) M=0.02 Fig. 14. Comparison of pressure drag coefficients by activating and deactivating the KCFSE. (a) M=0.05, simulation failed for pb*=10 without KCFSE and (b) M=0.02, simulation failed for pb*=10 and pb*=200 without KCFSE. t*CDp 0 2 4 6 80246810 pb*=0 pb*=10 pb*=200 Benchmark[33] t*CDp 3 3.5 4 4.5 511.522.5 t*CDp 0 2 4 6 802468101214 pb*=10 (KCFSE:off) pb*=10 (KCFSE:on) pb*=200 (KCFSE:off) pb*=200 (KCFSE:on) Benchmark[33] t*CDp 0 2 4 6 802468101214 pb*=10 (KCFSE:off) pb*=10 (KCFSE:on) pb*=200 (KCFSE:off) pb*=200 (KCFSE:on) Benchmark[33] M. H. Moghimi & N. J. Quinlan Engineering Anal ysis with Boundary Elements, December 2018 Page 16 of 21 C. Oscillating Rectangular Tank with Physical Free Surface In this test case, a rectangular tank co ntains liquid with a free surface, and o scillates vertically with displacement y = A sinሺ2π f tሻ, where A, f, and t are amplitude, frequency, and time, respec tively. Instead of moving the tank walls, the problem is solved in the reference frame of the tank, and the necessary fictitious force Fo=ρiViቀ-A(2πf)2sinሺ2πftሻቁ . (16) is applied to each particle i. These results were obtained for Reynolds num ber 2011 (based on the liquid depth and the amplitude of velocity oscillations, 2 πfA), LA⁄ = 200, and L∆x⁄ = 36, where L is the liquid depth. The Mach number M = ሺ2πfAሻc0⁄ is 0.02. The ratio h∆x⁄ is 0.485. The reference dynamic pressure pd and dimensionless background pressure pb* are defined as 2 ρሺπfAሻ2 and pbpd⁄ respectively. This is a challenging test case with displacements less than the particle spacing, but strong pressure oscillations which can trigger instability. Result are shown in Figs. 15 and 16 for tf = 21 (the end of the 21st cycle), with and without KCFSE applied, and various values of background pressure. With pb* = 0, large voids form near the floor of the tank (Fig. 15(a)). They are associated with regions of negative pressure (Fig. 16(a)), caused by rarefaction waves generated by downward acceleration of the floor, which interact with waves reflect ed from the free surface. Without KFSCE, as shown in Figs. 15(a-c), increasing background pressure does not improve the situ ation. Particles at the surfaces of internal voids are identified as part of the free surface and subjected to pb as an atmospheric pressure. Results for elevated pb* along with KCFSE are shown in Figs. 15(e-f) and 16(b-c). At pb* = 50, there are no major voids. However, regions of negative pressure remain and the pressure field is not uniform in the x direction. At pb* = 200, there are no voids, pressure is positive everywhere, and there are no variations in the x direction, as expected. (a) pb*=0, KCFSE off (b) pb*=50, KCFSE off (c) pb*=200, KCFSE off M. H. Moghimi & N. J. Quinlan Engineering Anal ysis with Boundary Elements, December 2018 Page 17 of 21 (d) pb*=0, KCFSE on (e) pb*=50, KCFSE on (f) pb*=200, KCFSE on Fig. 15. Particle distributions with LA⁄=200, tf =21, L∆x⁄=36 at Re=2011, M=0.02 by activating and deactivating the KCFSE. Identified free-surface particles are colored red. (a) pb*=0, KCFSE on (b) pb*=50, KCFSE on (c) pb*=200, KCFSE on Fig. 16. Pressure distributions with LA⁄=200, tf =21, L∆x⁄=36 at Re=2011, M=0.02 with KCFSE for various background pressures. D. Dam-break Problem In this test case [34], the effectiven ess 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 2 H respectively. The Reynolds number based on the liquid length and characteristic velocity ඥgH is 1.03 × 106. The ratio H∆x⁄ is 60, where ∆x is the initial particle spacing. The Mach number M = ඥgHc0ൗ is 0.05. The ratio h∆x⁄ is 0.95. The dimensionless background pressure pb* is defined as pb⁄ ܪ݃ߩ and dimensionless time t* is defined as t ඥgH⁄ . In Fig. 17(b) the KCFSE is not applied, and particles at unphysical voids are identif ied 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 pb* = 1, the result confirms the capability of the KCFSE to differentiate between physical and nonphysical free surfaces. M. H. Moghimi & N. J. Quinlan Engineering An alysis with Boundary Elements, December 2018 Page 18 of 21 (a) (b) pb*=1, KCFSE off (c) pb*=1, KCFSE on Fig. 17. (a) Configuration of dam-break numeri cal test. Particle distribution with H∆x⁄ =60 at pb*=1, Re=1.03 × 106, M=0.05, t*=6.5, (b) without KCFSE and (c) with KCFSE. Identified free-surface particles are coloured red. 5. CONCLUSIONS The application of background pressure has been investigat ed in the finite volume particle method (FVPM). It has been shown that in FVPM, numerical voids (regions not cove red 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 pr evious 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 KC FSE 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 M. H. Moghimi & N. J. Quinlan Engineering Anal ysis with Boundary Elements, December 2018 Page 19 of 21 enables background pressure to be applied in a flow with free surface. In principle, this a pproach 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 numeri cal void formation, with and without free surfaces, and results are in good agreement with be nchmark data where appropriate. ACKNOWLEDGEMENT This publication has emanated from research conducted wi th the financial support of Science Foundation Ireland (SFI) and Aerogen®, co-funded under the European Regional Devel opment Fund, under Grant Number 13/RC/2073.