A model for surface tension in the meshless ﬁnite volume particle method without spurious velocity

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 current associated with inaccurate curvature determination in some methods.


INTRODUCTION
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 straightforward 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 trivial in principle to track the surface of the fluid.In addition, 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 technique in computational fluid dynamics, introduced by Hietel et al. in 2000 [1].It preserves 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 welldefined 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.The 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 drops 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 methods cited above display some spurious or parasitic velocity.
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.

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.ρ, 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: where ρ 0 , c 0 , 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 where N , V i , β ij , U ij (t) , ẋi j , Ψ i , and dS are number of neighbours 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(U i ,U j ) to denote a numerical approximation to the inviscid numerical flux F(U ij )(t , x) -U ij (t , x) ẋi j (t), Eq. ( 6) is written as where β i b and F i b are boundary interaction vector and numerical flux approximation respectively.
Both the flux and the interface velocity depend on the particle transport velocity ẋi, which can be chosen arbitrarily.In the present work, we set ẋi=u i , i.e. a purely Lagrangian model.The interparticle area is defined as: where In the above, W j , σ and Ψ i are the kernel function, kernel summation and test function, respectively, defined as follows: In Eq. (10), m i 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: where Ω i is the support of particle i.
To close the system of equations, the evolution of the particle volume is calculated using where V i (t) is defined as V i (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 zeroorder 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.( 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 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 conditions, the integral term in the RHS of Eq. ( 6) is used.A numerical flux over the boundary segment is prescribed and then the boundary interaction vector β i b is calculated by integration (Eq.( 15)) between intersection points P a and P b as illustrated in Fig. 1.In Eq. ( 15), n i b and η are the boundary normal vector and coordinate along the boundary segment, respectively.

C. Free-Surface Detection in FVPM
Neighbouring particles and their weight functions 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 supports 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 particle is on the free surface (Fig. 2(b)).
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. 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 shown in Fig. 3.The cohesive surface tension force may be written as where σ LV and l ij , 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.Particles 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 where n ij 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.In the present 2D implementation, the interface length l ij is 1.To extend it to 3D, spherical particles should be considered, following Jahanbakhsh et al. [19], with l ij 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 P f .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 follows [20]: 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 the tangent vector of the wall as n t , the corresponding tangential component of adhesive force in 2D is written as 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 angles are equal.is the total cohesive force on particle i here only due to particle k).

NUMERICAL RESUTLS
The method is applied for initially square 2D droplets 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 c 0 is 80 m/s.The maximum Mach number u max c 0 ⁄ is 0.034.
In the absence of external forces, the initial square arrangement of the particles becomes a circular droplet through a sequence of oscillations.The evolution of the 2D viscous droplet is illustrated in Fig. 6.For an infinite, inviscid cylindrical drop the period of small-amplitude oscillation is given by Rayleigh [2,21] as 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 , where x i , x �, and R are the particle's compact support centre, the centre of mass of the steady-state droplet, and the droplet radius based on volume, respectively.The centre of mass of steady state circular droplet x � is calculated as follows: Pressure error ε p * is defined as where p NUM and p AN are the steady-state numerical and analytical pressure.The analytical internal pressure for an infinitely long cylindrical liquid drop under surface tension in zero ambient pressure is given by where . 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 �p NUM max -p NUM min � p AN � 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.04 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 Figs. 8 and 9. Pressure error ε p * decreases linearly as a function of dimensionless particle spacing ∆x * = ∆x L ⁄ .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 circumference of the exact solution.for first-order and secondorder 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.

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 where ϕ is ,and x '  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 drop.

C. Collision of Similar Droplets
The results of binary collisions are presented here.Fig. 15 shows a schematic view of the collision modes.The ratios L ∆x ⁄ and   ∆x ⁄ for each drop are 50 and 1.2 respectively.The reference sound speed respectively.In off-centre collision mode the spacing between the droplet centres, normal to the relative velocity, is equal to the radius, as depicted in Fig. 15.
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 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 quasilemniscate droplet.Due to the high surface tension force at the top and end of the formed quasilemniscate, 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, which 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].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 40 o to 140 o .The kinetic energy normalized to maximum kinetic energy for contact angles from 40 o to 140 o 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 60 o and 150 o .The SPH method of Liu et al. [23], in a similar test case for contact angles 60 o and 120 o , resulted in final kinetic energy, normalized by initial surface tension energy, on the order of 10 -3 .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: A θ e -sin(θ e )cos(θ e ) , 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.

E. Impact of Small Cylinder on a Liquid Surface
In this section, the new method is used to simulate a strongly dynamic, approximately twodimensional 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 where l c =�σ ρg ⁄ is the capillary length.The dimensionless radius R and weight W are defined as  l c ⁄ and πr 2 ρ s g σ ⁄ 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 diameter of the cylinder and the impact velocity, for cases on the threshold between floating and sinking, varies between 151 and 322.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 deceleration and allow the body to float for higher impact velocity.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.

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 static 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, this 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 intended, using spherical particles.

Fig. 1 .
Fig. 1.Boundary normal vector n b i of truncated support of particle i.

Fig. 2 .
Fig. 2. Summation of interparticle area vectors between particle i and its neighbours is (a) zero and (b) non-zero.

Fig. 3 .
Fig.3.Schematic view of surface tension force vectors between particle i and j (cohesive force) at point P a , and the freesurface 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 resulting 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.

Fig. 4 .
Fig. 4. Schematic diagram of surface tension force between particle i and a solid wall (F i a ) at point P f in (a) non-equilibrium (b) equilibrium state (F i cis the total cohesive force on particle i here only due to particle k).

Fig. 7
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

Fig. 7 .
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 particles' compact support centres, normalized to final droplet radius) for 2D viscous droplets for various Ohnesorge numbers µ �ρσR ⁄ evolved from square configurations, L ∆x ⁄ =25.

Fig. 8 .
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.

Fig. 9 .
Final states of initially square drop for various ∆x L ⁄ in a convergence study (red points are free-surface particles).

Fig. 10
Fig.10shows the decay of dimensionless kinetic energy KE * for viscous fluids with nominally firstorder (zero-order reconstruction) and second-order (linear MUSCL reconstruction) schemes.In the firstorder 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 kinetic 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-

Fig. 10 .
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.

Fig. 14 .
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.
and R are the velocities and radius of the droplets.The magnitude of relative velocity |V rel | is |V L -V R |, and the dimensionless relative velocity |V rel | � dimensionless time t * used here is t |V rel | (2R) ⁄ .The Mach number |V rel | c 0 ⁄ is 0.043.The Weber and Reynolds numbers based on the relative velocity and the diameter of the droplets are 10 and 56

Fig. 15 .
Fig. 15.Schematic view of droplets collision for head-on and off-centre modes.

Fig. 20 .Fig. 21 .
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 contact area diameter D and height H of liquid drop on a solid wall at equilibrium state as a function of contact angle θ e .

Fig. 23 .
Snapshots of the motion of a circular cylinder with Re=220, F F =3.075, F S =3.342, and W=0.36 (red points are freesurface particles).
The ratio 2r ∆x ⁄ is 16.The reference sound speed c 0 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 80 o .Test conditions and results are shown in Table1and Fig.22for 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 F F and sinking Froude number F S 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,

Table 1 .
[24]arison of experimental data[24], FVPM calculations, and analytical results[24].F F and F S denote the maximum Froude number observed for floating and the minimum Froude number for sinking, respectively, at the corresponding W, R and Re.