Comparison of Commercial and Open-Source CFD Solvers on Surface Tension Dominated Flows

Problems involving multiphase flows require a physical understanding of how the phenomena develop and the specific interactions they manifest. For surface tension relevant flows, such as bubbles and droplets, the importance of modelling and predicting small-scale behaviour is crucial in accurately defining the liquid-gas interface and complex interactions that may take place. Axisymmetric numerical simulations of single droplets impacting onto thin liquid films are performed using commercial and open-source CFD codes. ANSYS Fluent® is the commercial software employed, whereas Basilisk is the open-source CFD solver adopted. The incompressible Navier-Stokes equations are coupled and handled differently throughout each software. A solution-adaptive mesh refinement tool is adopted to reduce computational cost. Software comparison is based on quantitative and qualitative analysis, namely crown height and outer diameter measurements, and the crown curvature and occurrence of splashing, respectively. Results show that Basilisk simulations are in good agreement with the experimental data. Fluent follows the tendency for the crown outer diameter however, in terms of height, the model under-predicts its growth and collapses at later stages of the impact for lower thicknesses.


Introduction
The occurrence of multiphase flows is characterised by the presence of two or more distinct phases, namely solid, liquid and vapour. The interaction between these phases allows for a multitude of flows, such as gas-liquid, gas-particle, liquid-solid and three-phase flows. These flows are described as extremely complex due to the underlying dynamics and interactions at the gas-liquid interface, as well as phenomena of heat transfer and phase change. The phenomenon of droplet impact is a common nature occurrence and found in various industrial applications, including drop impingement in internal combustion engines with direct fuel injection, sprinkling irrigation, spray cooling, ink-jet printing, pesticide spraying of crops, etc. [1,2]. As a multiphase flow, defining how to tackle the specific dynamics of the phenomena and determining the optimal approach is crucial in numerically modelling and solving these types of flows. According to Popinet [3], Eulerian fields are particularly suited for fluid mechanics due to the necessity of dealing with large deformations. These can either be coupled with a Lagrangian or Eulerian representation of the interface. The Lagrangian approach includes the Marker-and-Cell [4] and front-tracking methods, whereas an Eulerian approach consists of the Volume of Fluid (VOF) [5] and the Level-Set (LS) methods [6], among others. The major differences of these interface representations relate to the transport phenomena and topology changes. In the Lagrangian formulation, transport is rather simple, however it presents difficulties in large deformations such as ligament breakup and coalescence. On the contrary, the Eulerian formulation can handle topology changes accurately, however, leading to more complex transport schemes [3]. The difficulty in simulating two-phase flows is accurately defining a sharp interface. The VOF method [5] is based on a fractional volume of fluid for tracking the liquid-gas interface. This interface lies within a control volume with a volume fraction that varies between 0 and 1. This method is particularly suited for mass conservation, however inaccurate in calculating the normal vector and curvature, among other smooth properties, leading to a low resolution interface. On the contrary, the LS method [6] treats the interface as a zero level-set of a continuous function. As a smooth geometric function, this provides increased accuracy in the normal vector and curvature calculation at the expense of mass conservation. Due to this, several methods have been arising in order to find solutions to these disadvantages. The Coupled Level-Set and Volume of Fluid Method (CLSVOF) [7] combines the advantages of the VOF and LS methods, using the VOF method to calculate the volume fraction field and the LS to compute the geometrical properties at the interface. The conservative level-set method, developed by Desjardins et al. [8], employs a hyperbolic tangent level-set function. The transport and re-initialisation of the hyperbolic tangent function occur by means of high order fully conservative schemes, reducing mass conservation issues in level-set methods. The normal vector and curvature calculation can also be estimated using Height-Functions [9]. These geometric functions compute the curvature from volume fractions/level-set function directly on the interface, where the VOF/LS methods estimate the curvature on the cell centres, leading to a better approximation. The selection of the most appropriate models regarding applications involving two-phase flows is crucial in accurately modelling and simulating the specific phenomena and interactions they manifest. The main objective of this work is to numerically simulate droplet impact onto a liquid film using commercial and open-source CFD codes, namely ANSYS Fluent and Basilisk, respectively. Quantitative and qualitative analysis regarding crown height and outer diameter measurements, and crown evolution and splashing occurrence, respectively, are presented.

Numerical Methods
The numerical model consists of solving the Navier-Stokes equations coupled with an interface tracking method in an axisymmetric assumption. The general form of the governing equations for the mass and momentum conservation are represented by equations (1) and (2), respectively: where ⃗ U is the velocity vector, t is time, ρ is the density, p is the static pressure, µ is the viscosity, ⃗ g is the gravitational acceleration constant, and F σ are the surface tension forces. The VOF [5] and LS [6] methods are defined by the volume fraction, α, and the level-set function,ϕ, and can be expressed by equations (3a) and (3b), respectively: for the gas-phase 0 < α < 1 for the interface 1 for the liquid-phase where d is the distance from the interface. The general transport equations for the volume fraction and level-set function are displayed by equations (4) and (5), respectively: The continuum surface force (CSF) model [10] is adopted for the different solvers for the surface tension calculation. For Basilisk, the surface tension forces are based on the volume fraction, as presented by equation (6a), while in Fluent, these forces are based on the level-set function, as displayed by equation (6b): where σ is the surface tension, κ is the interface curvature, H ϕ is a Heaviside function used to smooth the surface tension calculation, δ (ϕ) is a delta function, and ⃗ n is the unit normal vector to the interface. The Heaviside and delta functions are defined by equations (7) and (8), respectively: |ϕ| > a for the liquid-phase |ϕ| ≤ a for the interface where a is the thickness of the interface. In terms of unit normal vector and curvature, each solver has a unique approach to the approximation of these geometrical properties. Fluent calculates these values based on the level-set function, as given by equations (9a) and (10a). Basilisk formulates these properties in a distinct manner, making use of height functions as presented by equations (9b) and (10b), respectively: where h x and h xx are the first and second derivatives of the height function, respectively. The physical properties of density and viscosity, similarly to the unit normal vector and curvature, also differ throughout the solvers. For Fluent, these properties are calculated as an arithmetic mean considering the level-set function, as shown by equations (11a) and (12a). For the Basilisk solver, the density is calculated similarly to Fluent, however employing the volume fraction scalar, as presented by equation (11b). In terms of viscosity, a harmonic mean of the viscosity is applied for increased stability, visualised by equation (12b): where the subscripts l and g refer to the liquid and gas phase, respectively. Figure 1 displays a schematic of the physical phenomenon. A single droplet of diameter D 0 and impact velocity U 0 impacts onto a thin liquid film of thickness h. The droplet and the liquid film are of the same fluid and surrounded by gas. The domain consists of a quadrilateral with four different boundary conditions. The bottom boundary is the impact surface defined as a no-slip wall with a static contact angle, θ w . The left boundary is specified as the axis of symmetry and it is used as the centerline for the axisymmetric numerical simulations. The top and right boundaries are defined as pressure-outlets to represent the far field. Gravity acts on the negative direction of the y-axis, g y = −9.81 m/s. The height and radial distance of the domain are specific to each solver. The fluid adopted for the numerical simulations is a 75%/25% jet fuel and biofuel mixture, more specifically The mesh adopted for the different solvers is an orthogonal cartesian mesh with adaptive mesh refinement (AMR). The base mesh has an initial refinement equal to 3 cells per diameter (cpd) of the droplet. For Fluent, the AMR is based on a gradient adaption approach based on the volume fraction scalar, whereas Basilisk employs an adaptive wavelet algorithm based on the estimation of numerical discretisation errors. The time-step is changed according to the Courant-Friedrichs-Lewy (CFL) condition. The solution methods differ for each of the solvers and were specifically defined for improving the numerical simulations. In terms of pressure-velocity coupling, Fluent adopts the PISO scheme and Basilisk employs a projection method. For the advection scheme, Fluent uses the QUICK scheme, whereas Basilisk employs a momentum-conserving advection scheme. Both solvers adopt PLIC-based reconstruction schemes for the liquid-gas interface. Table 1  For the domain, Fluent is more flexible in defining the height and radial distance of the domain in comparison with Basilisk, which does not allow the combination of non-cubic domains and MPI parallelism with tree grids [12]. Due to this, the dimensions may be analysed differently. For Fluent, the height is defined as 5D 0 , allowing for the phenomena visualisation. However, in terms of radial distance, the influence of the radial pressure outlet on the liquid film must be considered. Therefore, distances of 4D 0 , 5D 0 and 6D 0 were considered, achieving radial independence for a distance of 5D 0 . For Basilisk, due to the cubic domain and quadtree grid restrain, the height and radial distance must be equal. Therefore, following a similar approach, distances of 2.7D 0 , 5.3D 0 and 10.7D 0 defined the domain study, reaching independence for 5.3D 0 . For the mesh study, and for both solvers, the measurements were performed for a mesh resolution of 48, 96 and 192 cpd, obtaining mesh independence for a 96 cpd resolution.
Analogously, in terms of time, for Fluent, the measurements were performed for CFL numbers of 0.5, 0.2 and 0.1, reaching time independence for CFL = 0.2. For Basilisk, in order to ensure numerical stability, the CFL number is defined as CFL = 0.5 [12].

Results and Discussion
The validation of the axisymmetric numerical simulations was performed using available experimental data [11].  Figure 3a, that the Basilisk and Fluent numerical dimensionless outer crown diameter measurements are slightly lower than the experimental values. Despite that, for both cases, the measurements follow the tendency of the experimental case. In Figure 3b, the Basilisk height values are in good agreement with the experimental case. For Fluent, the height measurements are under-predicted for the earlier phase of the impact and, at τ = 7, the crown disintegrates, leading to a sudden decrease of the crown height and to the formation of secondary droplets. Figure 4 displays the visualisation of the droplet impact phenomena for different impact stages between the experimental and numerical results for the lower thickness case, h * = 0.2. For the Basilisk solver, there is a smooth crown formation with secondary atomisation at τ = 2, as well as bubble entrapment. The crown shape and curvature is similar to the experimental crown. For τ = 5 and τ = 8, the crown continues to develop, increasing its height and outer diameter and exhibiting a thin liquid crown without ligament breakup. On the other hand, Fluent presents a distinct behaviour. For τ = 2, the crown early formation is quite similar, displaying secondary atomisation, and bubble entrapment in the impact region. For τ = 5, the occurrence of secondary atomisation from the crown rim is visible and, in comparison with the Basilisk results, the height presents a lower value, as verified in Figure 3b. At τ = 8, the disintegration of the crown is occurring, relating to the sudden height decrease previously mentioned.    Figure 5b, the Basilisk solver predicts the crown height evolution, displaying a similar rate of change with the experimental data. Fluent also displays a similar behaviour in terms of tendency, however the overall values are still under-predicted. Dissimilar to the previous case, the crown does not collapse for later stages of the impact, maintaining a regular growth. Figure 6 shows the crown outer diameter and height measurements for the higher thickness case, where D 0 = 3.0 mm, U 0 = 3.0 m/s, and h * = 1. Similarly to the previous case, the crown outer diameter numerical measurements in Figure 6a are in good agreement with the experimental results. In terms of growth rate, for later stages of the impact, Fluent shows a smaller discrepancy in which the values tend to stabilise for lower dimensionless time values in comparison to the experimental trend, whereas Basilisk displays a similar tendency to the experimental case throughout the different impact stages. The height measurements in Figure  6b follow the same trend as the previous case, with a smaller under-prediction for Fluent. It is possible to visualise that the numerical results are in a better agreement for higher dimensionless thicknesses, h * = 0.5 and h * = 1, in comparison to the lower thickness value, h * = 0.2, both for height and outer diameter measurements. This implies that further studies are required to understand thin liquid sheets and breakup, and such influence on crown geometrical parameters.  The simulations were performed in parallel on an octa-core local machine with 64GB of memory and hyper-threading enabled. The wall-clock time for the Basilisk and Fluent simulations are 5 − 6 minutes and 1.2 − 1.3 hours, respectively. Due to this, the Basilisk solver is approximately 15 times faster than Fluent, thus computationally more efficient.

Conclusions
Axisymmetric numerical simulations of single droplets impacting onto liquid films were performed using commercial and open-source CFD solvers, namely ANSYS Fluent and Basilisk, respectively. The Basilisk results are in good agreement and follow the tendency of the experimental data, despite being slightly lower regarding crown outer diameter measurements. This is also visualised for Fluent, however, in terms of height, the model under-predicts its growth. Moreover, for h * = 0.2, the crown breaks into secondary atomisation and the results displayed are less accurate in comparison with h * = 0.5 and h * = 1, meaning that further work is required regarding thin liquid ligament formation and breakup. The unit normal vector and curvature calculation are an indicative of the difference of the results obtained for the CFD solvers. Basilisk is superior in computational efficiency in comparison with Fluent.

Acknowledgements
The