A surface tension model has been developed in the finite volume particle method (FVPM). FVPM is a conservative, consistent, meshless particle method that incorporates properties of both smoothed particle hydrodynamics and the mesh-based finite volume method. Surface tension force is applied only on free-surface particles, which are inexpensively and robustly detected using the FVPM definition of interparticle area, analogous to cell face area in the finite volume method. We present a model in which the direction of the pairwise surface tension force is approximated by the common tangent of free-surface particle supports. The new surface tension model is implemented in 2D. The method is validated for formation of an equilibrium viscous drop from square and elliptical initial states, drops on hydrophobic and hydrophilic walls, droplet collision, and impact of a small cylinder on a liquid surface. Results are practically free from parasitic cu rrent associated with inaccurate curvature determination in some methods. Key Words: Finite volume particle method, meshless , surface tension, coalescence, wetting, cylinder impact on liquid . 1. INTRODU CTION In this article we describe a simple model for surface tension in the finite volume particle method (FVPM), with the aim of reducing or eliminating the spurious or parasitic velocity which is a feature of many surface tension models in both mesh-based and meshless methods. At present, the model is implemented in 2D , with validation in 2D static and dynamic problems. Many problems in fluid dynamics are dominated by surface tension , including spray coating, aerosolisation , and fuel injection. In mesh -based computational fluid dynamics, the accurate determination of interfaces, with or without surface tension, is not straigh tforward because interfaces do not generally conform to the mesh. In particle methods, in contrast, the computational nodes or particles are free to move, making it trivia l in principle to track the surface of the fluid. In addition, Page 1 of 23 M. H. Moghimi & N. J. Quinlan Computers & Fluids, March 2018 FVPM in particular allows definitive identification of particles at the free surface, greatly simplifying the implementation of surface tension models. FVPM is a relatively novel numerical tec hnique in computational fluid dynamics, introduced by Hietel et al. in 2000 [1]. It preserve s the conservative flux- based formulation of the classical finite volume method (FVM). FVPM particles are closely analogous to cells in the FVM as they have well- defined volume and interparticle area, and interact through numerical fluxes. However, they are allowed to overlap each other and move freely or arbitrarily, taking advantage of the Arbitrary Lagrangian -Eulerian (ALE) nature of FVPM . In the present work, a purely Lagrangian model of particle motion is strictly used. Particle methods to simulate surface tension may be categorized as microscopic or macroscopic. Microscopic methods in smoothed particle hydrodynamics (SPH) are based on application of van der Waals (vdW) -like potentials or other attractive -repulsive forces between particles [2,3 ,4]. This approach exploits the particle framework to mimic the molecular physics that underlie surface tension. The surface tension coefficient does not appear explicitly, and consequently the model requires calibration of coefficients . Van der Waals- based models are also sensitive to temperature . Yang et al. [5] proposed a hyperbolic kernel function with central discontinuity to improve particle distribution in this approach. Intermolecular forces manifest at macroscopic level as a net tensile stress in the plane of a free surface. In a planar free surface, these forces are in equilibrium. T he macroscopic continuum surface force (CSF) model [ 6] is therefore based on explicit calculation of the interface curvature in a two -phase system. This enables the macroscopically measured surface tension coefficient to be used in the model. Morris [7], Muller et al. [8], Liu et al. [9], Adami et al. [10], Breinlinger et al. [11], and Schnabel et al. [12] applied this technique in SPH. The method is sensitive to the numerical calculation of gradients of the color function . Maertens et al. [13] introduced surface tension modelling in FVPM. They simulated 3D drop s in isolation and interacting with solid surfaces, using the physical surface tension coefficient, without explicitly estimating the curvature of the interface or tuning any parameters. This is possible in FVPM because free -surface particles can be identified unambiguously and the free surface itself is well- defined in terms of the particle geometry. In this technique, the surface tension force is applied at the intersection curves of the free surfaces of spherical particles , tangent to the surface . All the met hods cited above display some spurious or parasitic velocity. Page 2 of 23 M. H. Moghimi & N. J. Quinlan Computers & Fluids, March 2018 In this article, a variant model for surface tension in FVPM is presented . Section 2 summarizes the FVPM formulation for fluid flow computations. In section 3, the new surface tension model in FVPM is described. In section 4, numerical results are presented and discussed. 2. FVPM METHOD A. Governing Equations The Navier -Stokes equations without source terms can be written in conservation form as ∂U(t , x) ∂t+∇ . F �U(t , x)�=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 . ρ, u, p, I, and τ are density, velocity, pressure, identity tensor, and viscous stress tensor of the fluid, respectively. Equation (5 ) is used for pressure calculation : p = ρ0c02 γ��ρ ρ0�γ -1� , (5) where ρ0, c0, and γ are reference density of the fluid, reference sound speed, and a constant, respectively. The value of γ is set to 7 in the present work [14]. B. Numerical Model The FVPM formulation of Eq. (1) for particle i is d dt�Ui(t)Vi(t)� = -∑βij . �F �Uij�(t , x) - Uij(t , x) ẋij(t)�N j=1 -∫Ψi(t , x)F �Uij�(t , x) . dS , (6) where N, Vi, βij, Uij(t), ẋij, Ψi, and dS are number of neighbo urs of particle i, particle volume, interparticle area vector between particles i and j, vector of conserved quantities at the interface of particles i and j, velocity of the interface, test function of the particle, and the surface area vector of the support of particle i, respectively. Introducing F(Ui,Uj) to denote a numerical approximation to the inviscid numerical flux F(Uij)(t , x) - Uij(t , x) ẋij(t), Eq. ( 6) is written as Page 3 of 23 M. H. Moghimi & N. J. Quinlan Computers & Fluids, March 2018 d dt�Ui(t)Vi(t)� = -∑βij . �F �Ui , Uj��N j=1 -βib.Fib , (7) where βib and Fib are boundary interaction vector and numerical flux approximation respectively . Both the flux and the interface velocity depend on the p article transport velocity ẋi, which can be chosen arbitrarily. In the present work, we set ẋi=ui, i.e. a purely Lagrangian model. The interparticle area is defined as: βij = γij - γji , (8) where γij=∫Ψi(t , x)∇ Wj(t , x) σ(t , x)dx . (9) In the above, Wj , σ and Ψi are the kernel function , kernel summation and test function, respectively, defined as follow s: Wi(t , x) = mi W(x - xi(t) , h) , (10) σ(t , x) =∑Wj(t , x)N j=1 , (11) Ψi(t , x) = Wi(t , x) σ(t , x) . (12) In Eq. (10), mi and h are the particle’s mass and the smoothing length of the kernel function , with compact support radius 𝑅𝑅𝑐𝑐 equals to 2h. In the present work, W is defined as a top -hat function as the following: Wi(t , x)=�1 x ∈ Ωi 0 x ∉ Ωi , (13) where Ωi is the support of particle i . To close the system of equations, the evolution of the particle volume is calculated using d dtVi(t) =∑�γij . ẋj(t) - γji . ẋi(t)�N i=1 , (14) where Vi(t) is defined as Vi(t) =∫Ψi(t , x)dx. In the present work, particle supports are circles of constant radius . Interparticle inviscid fluxes of Eq. (3) are computed using the AUSM+-up scheme [ 15]. In first -order FVPM, based on the Godunov method, the flux between a pair of particles is computed using a zero - order consistent reconstruction of particle data from the barycentre to the particle -particle interface. In the nominally second -order version [15 ], based on the MUSCL finite volume approach, linear reconstruction is used to compute the fluxes. As discussed by Nestor et al. [16], the viscous terms (Eq. Page 4 of 23 M. H. Moghimi & N. J. Quinlan Computers & Fluids, March 2018 (4)) are computed at each particle -particle and particle -boundary interface . Velocity gradients are required to compute the viscous stresses and the linear reconstruction . Following Bonet and Lok [ 17], velocity gradients are calculated by using the consistency -corrected SPH approximation. Temporal discritisation of Eqs. (4) and (7) is based on the second -order Runge -Kutta scheme . The time step Δt determined by a Courant -Friedrichs -Lewy (CFL) criterion for FVPM is ∆t = CVi ∑��ui+uj 2�. βij+�ci+cj 2��βij�� j , (15) where C is the Courant number [18]. For detailed derivation of the FVPM formulation and further discussion of the extension of FVPM to viscous flow , readers are referred to [ 16] and [18] respectively. To impose boundary condition s, the in tegral term in the RHS of Eq. ( 6) is used . A numerical flux over the boundary segment is prescribed and then the boundary interaction vector βib is calculated by integration (Eq. ( 15)) between intersection points Pa and Pb as illustrated in Fig. 1. In Eq. ( 15), nib and η are the boundary normal vector and coordinate along the boundary segment , respectively. βib =∫Ψinibdη∂Ω . (16) Fig. 1 . Boundary normal vector nbi of truncated support of particle i. C. Free -Surface Detection in FVPM Neighbouring particles and their weight function s determine the interparticle area vectors βij, computed using Eq. (8) and (9). As illustrated in Fig. 2, particles j, k, m, n, and p are considered as the neighbours of particle i, because their compact support s overlap the compact support of particle i. If the boundary of a particle’s support is fully covered by neighbouring particles (Fig. 2 (a)), it can be shown from Eqs. (8) and (9) that its area vectors sum to zero ( ∑βij j =0). In contrast, if the summation of the interaction vectors between the particle of interest and its neighbours is not equal to zero (∑βij j≠0) it means that the particl e is on the free surface (Fig. 2 (b)). Page 5 of 23 M. H. Moghimi & N. J. Quinlan Computers & Fluids, March 2018 By using the top -hat kernel function βij can be calculated exactly [ 18]. In practice, for fully immersed particles, ∑βij j =0 is satisfied within some rounding error. Free -surface particles are identified when the magnitude of ∑βij j is greater than a small tolerance. (a) (b) Fig. 2. Summation of interparticle area vectors between particle i and its neighbours is (a) zero and (b) non -zero. 3. SURFACE TENSION MODEL A. Surface Tension (Cohesion) To calculate the surface tension force, the free -surface particles must be detected , and the appropriate force applied between particles. A free -surface particle pair is show n in Fig. 3. The cohesive surface tension force may be written as Fij c= σLVlijn , (17) where σLV and lij, are the surface tension coefficient between liquid and gas and the length of the intersection curve of the particle surfaces (or unit out-of -plane length in 2D) . The unit vector n is the direction of the surface tension force , which is to be determined . Particle s i and j exert equal and opposite cohesion force F ij on each other. At macroscale, surface tension acts in a plane tangent to a free surface . In FVPM , the free surface is composed of the exposed surfaces of the particle support volumes, giving a bumpy approximation to the true free surface. In the present work , the common tangent of particles i and j is taken as an approximation for the free surface tangent. If the particles are of equal size, this is parallel to the compact supports’ centre -to-centre line. The general expression for cohesive force between particle i and its neighbors is then Page 6 of 23 M. H. Moghimi & N. J. Quinlan Computers & Fluids, March 2018 Fi c= ∑σLVlijnij N j=1 , (18) where nij is defined as 𝐱𝐱𝑗𝑗−𝐱𝐱𝑖𝑖 �𝐱𝐱𝑗𝑗−𝐱𝐱𝑖𝑖� , and the sum is taken over all neighbours on the free surface . In contrast to other macroscopic methods, in FVPM the curvature need not be calculated explicitly. Fig. 3 . Schematic view of surface tension force vectors between particle i and j (cohesive force) at point Pa, and the free - surface common tangent (dashed line) of the particles . Maertens et al . [13] derived the interparticle force in 3D by a different approach, considering changes in surface energy of the particles’ exposed spherical surface patches due to motion of the particles. This yields forces tangent to the particle support boundaries. For a pair of equally sized particles, the result ing force reduces to a centre -to-centre vector between the particles, as in the present work, with a geometric correction factor. Maertens et al. employed an ALE particle transport velocity correction , while in the present work the particles are advected with purely Lagrangian velocity . In the present 2D implementation, the interface length lij is 1. To extend it to 3D, spherical particles should be considered , following Jahanbakhsh et al . [19], with lij computed as the length of arcs of intersection of the particles’ free surfaces. B. Surface Tension (Adhesion) In the presence of a solid -liquid interface, the adhesive force is applied at the intersection point on the wall , point Pf. In Fig. 4, σLV, σSV, and σSL are the surface tension coefficients of the liquid -gas interface, solid -gas interface, and solid -liquid interface respectively. Young’s relation links these coefficients with the equilibrium contact angle as follow s [20]: σSV - σSL= σLV cos(θe) . (19) The normal components of the total cohesive and adhesive forces are in balance. Only the tangential components are responsible for the movement of the contact line. Defining t he tangent vector of the wall as nt, the corresponding tangential component of adhesive force in 2D is written as Fi,t a= �|Fi c|cos(θd)- σLV cos(θe)lij�nt , (20) Page 7 of 23 M. H. Moghimi & N. J. Quinlan Computers & Fluids, March 2018 where θd is defined as the instantaneous dynamic contact angle between the cohesive force vector and the wall . In equilibrium, the dynamic and equilibrium contact angle s are equal . (a) (b) Fig. 4 . Schematic diagram of surface tension force between particle i and a solid wall ( Fia) at point Pf in (a) non -equilibrium (b) equilibrium state (Fi c is the total cohesive force on particle i here only due to particle k). 4. NUMERICAL RESUTLS The method is applied for initially square 2D droplet s in free space and on a solid wall , as problems expected to result in static equilibrium . To test the method in strongly dynamic cases, we simulated oscillation of an elliptical drop, collision of two similar drops , and impact of a solid cylinder on a liquid surface . A. Oscillation of a 2D Liquid Square Droplet The first test case is formation of a circular drop from an initial square distribution of liquid particles (Fig. 5) with surface tension coefficient σ, density ρ, and dynamic viscosity µ. The Reynolds number based on the square ’s side length L and characteristic velocity σµ⁄ is 3633. The ratio L∆x⁄ is 25, where ∆x is the initial particle spacing. The radius of all particles is 𝑅𝑅𝑐𝑐 = (1.2)∆x. The reference sound speed c0 is 80 m/s. The maximum Mach number umaxc0⁄ is 0.034. In the absence of external forces , the initial square arrangement of the particles becom es a circular droplet through a sequence of oscillation s. The evolution of the 2D viscous droplet is illustrated in Fig. 6. For an infinite, invis cid cylindrical drop the period of small- amplitude oscillation is given by Rayleigh [2,21] as τ = 2π �ρ R3 6 σ . (21) Page 8 of 23 M. H. Moghimi & N. J. Quinlan Computers & Fluids, March 2018 The dimensionless time t* is t τ . In steady state (t*=0.049) , gaps appear in the particle distribution near the initial corners of the drop . However, there is no void inside the drop , i.e. the particles fully cover the volume of the drop . The method is reaching equilibrium despite this non -uniform particle distribution. The dimensionless radius r* is defined as and |xi - x�| R, where xi, x�, and R are the particle ’s compact support centre , the centre of mass of the steady -state drop let, and the droplet radius based on volume , respectively . The centre of mass of steady state circular droplet x� is calculated as follow s: 𝑥𝑥̅ =∑ xiρiViN i=1 ∑ρiViN i=1 , (22) 𝑦𝑦� =∑ yiρiViN i=1 ∑ρiViNi=1 . (23) Pressure error εp* is defined as pNUM - pAN pAN where pNUM and pAN are the steady -state numerical and analytical pressure. The analytical internal pressure for an infinite ly long cylindrical liquid drop under surface tension in zero ambient pressure is given by pAN =σ R . (24) Fig. 5. Initial arrangement particles to form a liquid droplet , L∆x⁄=25 (red points are free -surface particles). t*=0.000 49 t*=0.000 86 t*=0.0016 t*=0.0028 t*=0.00 32 t*=0.0047 t*=0.0053 t*=0.006 3 t*=0.009 t*=0.049 Fig. 6. Snapshots of 2D viscous liquid droplet evolution from an initial square shape, L∆x⁄=25 (red points are free -surface particles). Fig. 7 shows the distribution of dimensionless pressure for L∆x⁄=25 across drops with various Ohnesorge numbers µ �ρσR⁄ based on the drops’ radii R defined in terms of 2D volume as Page 9 of 23 M. H. Moghimi & N. J. Quinlan Computers & Fluids, March 2018 R =�VDrop π , (25) where VDrop is ∑Vi i=1 . It is worth noting that the dimensionless radius of the outermost particle is less than 1. This is because the plot shows 𝑟𝑟∗ at particle compact support centres, which lie inside the radius R. Dimensionless pressure error is approximately 4.6×10–3 or less. Pressure within the drop is practically uniform , with normalized pressure range �pNUMmax-pNUMmin�pAN� less than 4.8×10–11. These results show an improvement on prior SPH results. In a 2D two-phase model with a CSF method , Adami et al. [10] reported dimensionless pressure error up to 0.07 and with a range of 0.03 for L∆x⁄ = 24 and a density ratio of 1000. Zhang et al. [22] reported error of 0.0 4 with range of 0.03, for similar density ratio of 798 and resolution L∆x⁄ = 30. In FVPM simulations of a 3D drop with L∆x⁄= 10, Maertens et al. [13] found error of 0.023. For the same physical and numerical conditions, we have obtained a similar level of error. In comparison with other particle methods, the present method results in similar or lower error with a practically uniform and steady pressure field. Results of a convergence study are shown in Fig s. 8 and 9. Pressure error εp* decreases linearly as a function of dimensionless particle spacing ∆x*=∆xL⁄. This is expected, since there is a first -order error between the surface of the numerical drop (defined by free- surface arc segments) and the circular circum ference of the exact solution. Fig. 7. Normalized error εp* between the numerical pressure and the Laplace analytical pressure , for all particles , as a function of r*(radial position of particle s’ compact support centres , normalized to final droplet radius) for 2D viscous droplet s for various Ohnesorge numbers µ�ρσR⁄ evolved from square configuration s, L∆x⁄=25. Fig. 8. Pressure error εp* (discrepancy between the numerical pressure and the Laplace analytical pressure, normalized to the analytical pressure) as a function of particle spacing. Page 10 of 23 M. H. Moghimi & N. J. Quinlan Computers & Fluids, March 2018 ∆x L = 0.1 ∆x L = 0.05 ∆x L = 0.025 ∆x L = 0.0125 Fig. 9. Final states of initially square drop for various ∆xL⁄ in a co nvergence study (red points are free -surface particles). Fig. 10 show s the decay of dimension less kinetic energy KE∗ for viscous fluids with nominally first - order (zero -order reco nstruction) and second -order (linear MUSCL reconstruction) schemes. In the first - order scheme (Fig. 10 (a)), numerical dissipation dominates physical viscosity. In second -order scheme (Fig. 10(b)), decay is slower . Oscillation of total kinetic energy is less with the zero -order reconstruction than with the MUSCL reconstruction . To avoid excessive numerical dissipation , linear MUSCL reconstruction is used in the remainder of this work . Figs. 10(a) and (b) show practical elimination of spurious kinetic energy . The normalized kinetic energy KE*=KE 2πRσ⁄ (kinetic energy normalized to equilibrium surface energy) at steady state is on the order of 10–13 (Fig. 10(b)). Final k inetic energy normalized to maximum kinetic energy has the order of 10-11 (Fig. 10(b)). The order of residual dimensionless velocity u*=|u|�σµ⁄�⁄ for first -order and second - order schemes are 10–14 and 10–9 respectively . Detailed results show that this is not due to classical parasitic velocity, but rather a rigid -body motion of the drop . Adami et al. [10] proposed an approach to obtain a stable and accurate surface -curvature calculation based on CSF in SPH . The order of kinetic energy ( normalized to the maximum) was 10–3. Liu et al. [23], using a hybrid particle -mesh method for 2D drops, found final kinetic energy (normalized to initial surface energy) on the order of 10–3 to 10–2, depending on the initial configuration. In the 3D FVPM results of Maertens et al. [13], the final kinetic energy ( normalized to maximum kinetic energy ) ranges from 10–2 to 10–3 for a 3D drop, depending on sound speed. For the same fluid properties and particle spacing, we have found normalized residual kinetic energy of 10-9 for a 2D drop. Page 11 of 23 M. H. Moghimi & N. J. Quinlan Computers & Fluids, March 2018 (a) (b) Fig. 10. Decay of KE*=KE 2πRσ⁄ (kinetic energy normalized to final surface energy) as a function of t* (time normalized to period of droplet oscillation) for 2D droplet evolution from square configuration for various methods: (a) viscous fluid and first -order FVPM and (b ) viscous fluid and second- order FVPM. B. Oscillation of a 2D Elliptical Droplet The droplet was deformed from its circular final state in the last section to an elliptic configuration as the initial condition for a new test, following Nugent and Posch [2] and Yang et al. [5]. The deformation is applied using �x' y'�=�2 sinφ�xsin�ϕ 2� ycos�ϕ 2�� , (32) where ϕ is 𝑒𝑒𝑒𝑒,and x'and y'are new value s of x and y. The value of 𝑒𝑒 (eccentricity) is 0.55. Figs. 11 and 12(a) show the oscillating deformation of the surface of drop . In these figures , R* is the dimensionless deformation parameter defined as ∑�xi,FS - x��NFS i=1(RN FS)� where xi,FS and NFS are the position s of free - surface particles and their number . In the deformation of the droplet surface the effect of viscosity is evident , lead ing to faster damping of oscillations . Page 12 of 23 M. H. Moghimi & N. J. Quinlan Computers & Fluids, March 2018 (a) (b) Fig. 11. Surface deformation of 2D droplet evolved from elliptic configuration versus time : (a) viscous and (b) inviscid. (a) (b) Fig. 12. (a) Three cycles of surface deformation of 2D inviscid droplet evolved from elliptic configuration as a function of time and (b ) computed and analytical dimensionless oscillation period τ*as a function of surface tension coefficient σ*, where both are normaliz ed by the values in the first case τ*=σ*=1. To determine oscillation period, sinusoidal curves with decaying amplitude were fitted to the numerical results. The results, in Fig. 12(b), show good agreement of the FVPM computed oscillation period with the analytical solution ( Eq. (21)) for a range of surface tension coefficients. The evolution of the 2D viscous droplet from the initial elliptical shape is illustrated in Fig. 13. Kinetic energy of the drop is shown in Fig 14. The strong surface tension at the poles of the perturbed droplet leads to intense oscillation in kinetic energy. In the viscous drop the energy dissipation is faster than the inviscid dro p. Page 13 of 23 M. H. Moghimi & N. J. Quinlan Computers & Fluids, March 2018 t*=0 t*=0.0022 t*=0.0064 t*=0.011 t*=0.015 t*=0.054 Fig. 13 . The snapshots of 2D viscous liquid droplet evolution with initial ellipse shape with basic Lagrangian FVPM (red points are free -surface particles). (a) (b) Fig. 14. Decay of KE*=KE 2πRσ⁄ (kinetic energy normalized to initial surface energy) as a function of t* (time normalized to period of droplet oscillation) for 2D droplet evolved from elliptic configuration: (a) viscous and (b) inviscid . C. Collision of Similar Droplets The results of binary collisions are presented here. Fig. 15 shows a schematic view of the collision modes. The ratio s L∆x⁄ and 𝑅𝑅𝑐𝑐∆x⁄ for each drop are 50 and 1.2 respectively . The reference sound speed c0 is 150 m/s. In this figure, VL, VR, and R are the velocit ies and radius of the droplet s. The magnitude of relative velocity |Vrel| is |VL-VR|, and the dimension less relative velocity |Vrel|�σ µ�� is 0.18. The dimensionless time t* used here is t |Vrel|(2R)⁄ . The Mach number |Vrel|c0⁄ is 0.043. The Weber and Reynolds numbers based on the relative velocity and the diameter of the droplets are 10 and 56 respectively. In off -centre collision mode the spacing between the droplet centre s, normal to the relative velocity , is equal to the radius , as depicted in Fig. 1 5. Results for the head -on collision are shown in Fig. 16. The process of merging starts by forming a liquid bridge at the interaction interface . This section elongates in the vertical direction as the droplets Page 14 of 23 M. H. Moghimi & N. J. Quinlan Computers & Fluids, March 2018 continue their motion towards each other. The vertical expansion and horizontal contraction continue until the merged liquid reaches its maximum elongation in vertical direction and forms a quasi - lemniscate droplet . Due to the high surface tension force at the top and end of the formed quasi - lemniscate , it starts to contract from the poles. The merged droplet then oscillates until it forms a circular droplet. Fig. 17 shows a binary collision of drops in off -centre mode. The same process is seen in this collision as the liquid bridge forms at the initiation of the collision , and then expands in the direction normal to the initial centreline of the droplets. In contrast to the head -on collision , rotation of the liquid can be seen. The merged liquid forms a quasi -lemniscate droplet , whic h oscillates and rotates until it reaches a circular shape. The same behaviours of these collision modes have been reported by Melean et al. [3] and Yang et al. [5]. Fig. 15. Schematic view of droplets collision for head -on and off -centre mode s. t*=0 t*=0.0028 t*=0.0039 t*=0.0092 t*=0.012 t*=0.0166 t*=0.024 t*=0.0314 t*=0.037 t*=0.0406 t*=0.0535 t*=0.365 Fig. 16. Evolution of 2D head -on binary collision of two similar viscous droplets with |Vrel|�σ µ��=0.18, We=10, and Re=56 (red points are free -surface particles). Page 15 of 23 M. H. Moghimi & N. J. Quinlan Computers & Fluids, March 2018 t*=0 t*=0.0042 t*=0.0052 t*=0.00 81 t*=0.0129 t*=0.0148 t*=0.0185 t*=0.024 t*=0.0369 t*=0.0462 t*=0.148 t*=0.338 Fig. 17. Evolution of 2D off-centre binary collision of two similar viscous droplets with |Vrel|�σ µ��=0.18, We=10, and Re=56 (red points are free -surface particles). D. Liquid Drop on a Solid Wall In this section, t he results of droplet formation from a n initial truncated circular arrangement of particles on hydrophilic ( θe=40o, 60o, and 80o) and hydrophobic ( θe=120o and 140o) solid surfaces are presented (Fig. 18). The ratio s 2R∆x⁄ for equivalent circular drop and 𝑅𝑅𝑐𝑐∆x⁄ are 50 and 1.2 respectively . The reference sound speed c0 is 160 m/s. The test case with a contact angle of 40o resulted in the maximum Mach number umaxc0⁄ of 0.015. On the hydrop hilic surface (Fig. 18(b), (c), and (d)) , the base of the droplet extends on the solid surface, increase s the contact area between the wall and the liquid , while in the hydrophobic case (Fig. 18(e) and (f) ), the contact surface is decreased . The contact angles for all cases have been checked in the final condition using the coordinates of the particles near the wall (Fig. 4), and are in close agreement with the prescribed equilibrium contact angles . Fig. 19 shows the decay of the kinetic energy. The results in Fig. 19 show practical elimination of spurious residual kinetic energy at contact angles 40o to 140o. The kinetic energy normalized to maximum kinetic energy for contact angles from 40o to 140o ranges from 10-7 to 10-11. Breinlinger et al. [11] presented SPH simulations of a 2D drop on a flat surface , in which normalized kinetic energy normalized to maximum kinetic energy is 10-3 for contact angles 60o and 150o. The SPH method of Liu et al. [23], in a similar test case for contact angles 60o and 120o, resulted in final kinetic energy , normalized by initial surface tension energy , on the order of 10-3. Page 16 of 23 M. H. Moghimi & N. J. Quinlan Computers & Fluids, March 2018 (a) (b) (c) (d) (e) (f) Fig. 18. Evolution of v iscous liquid drop with initial truncated circular configuration (a) on hydrophilic (b, c, and d) and hydrophobic (e and f) surfaces , 2R∆x⁄=50 (red points are free -surface particles) . Fig. 19 . Decay of KE*=KE 2πRσ⁄ (kinetic energy normalized to initial surface energy) as a function of t* (time normalized to period of droplet oscillation) for a 2D viscous droplet on hydrophilic ( θe=40o, 60o, and 80o) and hydrophobic surfaces (θe=120o and 140o). Page 17 of 23 M. H. Moghimi & N. J. Quinlan Computers & Fluids, March 2018 The diameter D of contact area and height H of a two dimensional drop on a flat surface (Fig. 20) in the absence of external body forces are given by Breinlinger et al. [11] as follows : D= 2δ sin(θe) , (33) H= δ (1-cos(θe)) , (34) δ= �A θe-sin(θe)cos(θe) , (35) where A and δ are the cross -section of circular drop and the radius of curvature . Fig. 21 shows the variation of contact diameter and height of the drop with contact angles. The FVPM results agree closely with the analytical solution . Fig. 20 . Schematic diagram of a liquid drop on a solid wall at equilibrium state with height H, contact area diameter D, and contact angle θe. Fig. 21 . Normalized c ontact area diameter D and height H of liquid drop on a solid wall at equilibrium state as a function of contact angle θe. E. Impact of Small Cylinder on a Liquid Surface In this section, the new method is used to simulate a strongly dynamic, approximately two - dimensional experiment with strong surface tension effects . Vella et al. [24] provide experimental data and an approximate analytical solution for impact of a horizontal steel cylinder dropped onto a water free surface. In the experiment, the terminal velocity of the cylinder in air is much higher than the impact speed of the cylinder , so that the cylinder is in the early stages of free fall and aerodynamic forces can be neglected . The cylinder is released from rest and impacts with velocity U. The impact Froude number is defined as F=U (glc)12⁄ , (36) Page 18 of 23 M. H. Moghimi & N. J. Quinlan Computers & Fluids, March 2018 where lc=�σρg⁄ is the capillary length. The dimensionless radius R and weight W are defined as 𝑟𝑟lc⁄ and πr2ρsgσ⁄ respectively . The parameters r and ρs are the radius and density of cylinder respectively. The experiment was carried out with water and isopropanol- water mixtures. The Reynolds number based on the d iameter of the cylinder and the impact velocity, for cases on the threshold between floating and sinking, varies between 151 and 322. The ratio 2r∆x⁄ is 16. The reference sound speed c0 is 20 m/s. The Mach number based on the impact velocity, for cases on the threshold between floating and sinking, varies between 0.016 and 0.027. The density ratio ρsρ⁄ is 7.85. The equilibrium contact angle θe between the cylinder and the fluid is 80o. Test conditions and results are shown in Table 1 and Fig. 22 for the experimental and analytical results of Vella et al. and the present FVPM method. FVPM results shown are the maximum Froude number for which floating resulted and the minimum that resulted in sinking, at each value of weight W. The gap between floating Froude number FF and sinking Froude number FS is due to the interval between tested impact velocities. FVPM results for the boundary between sinking and floating regimes are close to both analytical calculations and experiments in the lower range of weight. At higher weight, FVPM and analytical models are in good agreement, but both predict lower critical Froude number for sinking than observed in experiment . Vella et al. [24] suggested that the differences between the experiment and the 2D analytical model may be due to surface tension force at the ends of the body , which can increase decelerat ion and allow the body to float for higher impact velocit y. The end effects are also absent from the present 2D computational model, and may explain the discrepancy from experiment. Snapshots of the cylinder motion are presented in Fig. 23. It is seen that after downward motion the cylinder rebounds upward until it floats on the surface of the liquid. Table 1. Comparison of experimental data [24], FVPM calculations , and analytical results [24]. FF and FS denote the maximum Froude number observed for floating and the minimum Froude number for sinking, respectively, at the corresponding W, R and Re. % isopropanol (vol.) W R Re Experiment [24] FVPM Analytical [ 24] FF FS FF FS F 0 0.25 0.10 299 4.7 4.9 3.303 3.548 3.448 0 0.3 0.11 306 4.3 4.5 3.120 3.365 3.228 6.25 0.36 0.12 220 3.5 4.0 3.075 3.342 2.991 0 0.41 0.13 322 2.7 3.0 2.814 3.059 2.826 10 0.66 0.16 151 1.7 1.9 2.127 2.401 2.191 Page 19 of 23 M. H. Moghimi & N. J. Quinlan Computers & Fluids, March 2018 Fig. 22. Floating -sinking regime diagram of cylinder impact on a liquid surface . t*=0 t*=0.0042 t*=0.014 t*=0.018 t*=0.024 t*=0.032 t*=0.052 Fig. 23. Snapshots of the motion of a circular cylinder with Re=220 , FF=3.075 , FS=3.342 , and W=0.36 (red points are free - surface particles). Page 20 of 23 M. H. Moghimi & N. J. Quinlan Computers & Fluids, March 2018 5. CONCLUSION A simple model for surface tension in the finite volume particle method (FVPM) has been presented in this article . It has been tested for 2D static and dynamic liquid flows. The results obtained are in good agreement with existing data. In stati c tests, drops in free space and on a wall reach equilibrium states practically free of spurious velocity, with uniform pressure and final kinetic energy on the order of 10–7 to 10–14 of initial surface energy. In simulations of cylinder impact on a liquid surface, predictions of floating or sinking outcome are in agreement with experimental data. All simulations were carried out with purely Lagrangian particle transport, and equilibrium was achieved despite notably non-uniform particle distribution in some cases. The only input parameters required are surface tension coefficient and equilibrium contact angle. FVPM has the advantage that the free surface can be detected definitively at no extra cost. In principle, t his surface tension approach could be implemented in smoothed particle hydrodynamics (SPH) using a method for free surface detection such as that of Marrone et al. [25]. In further work, the extension of the method to 3D is intend ed, using spherical particles. ACKNOWLEDGE MENT This publication has emanated from research conducted with the financial support of Science Foundation Ireland (SFI) , Aerogen®, and the European Regional Development Fund under Grant Number 13/RC/2073.