On Simulating Primary Atomization Using the Refined Level Set Grid Method

To simulate primary atomization, one has to track the position of the phase interface accurately, handle large numbers of topology changes and drops, treat the singular force of surface tension in an accurate and stable manner, and ensure grid-independent numerical results. To address all of these challenges we present a balanced force Refined Level Set Grid (RLSG) method for collocated, unstructured finite volume flow solver grids that can be coupled to a Lagrangian spray model. Special emphasis is placed on the accurate treatment of surface tension forces, since during the atomization of liquid jets by coaxial fast-moving gas streams, the details of the formation of small-scale drops from aerodynamically stretched out ligaments are governed by capillary forces [1]. Several different generic verification examples are presented, discussing the accuracy, volume preservation, and grid-convergence properties of the balanced force RLSG method. ∗Corresponding Author Introduction The simulation of the primary atomization process of liquid jets and sheets is a challenging numerical problem. One not only has to track the position of the liquid/gas interface and handle a large number of topology changes, but one also has to account for the fact that on scales typically associated with atomizers, the phase interface is a discontinuity and the surface tension force represents a singular force. Treating the surface tension force numerically in a stable and accurate manner is of crucial importance, since, for example, during the atomization of liquid jets by coaxial fast-moving gas streams, the details of the formation of small-scale drops from aerodynamically stretched out ligaments are governed by capillary forces [1]. From a numerical point of view, surface tension poses a unique challenge since it is a singular force, active only at the location of the phase interface where material properties, like density and viscosity, change discontinuously. One of the prerequisites for correctly treating surface tension forces is therefore the ability to locate the position of the phase interface accurately. To this end, several phase interface tracking schemes exist for fixed grid flow solvers, among them the marker method [2], the Volume-of-Fluid (VoF) method [3], and the level set method [4]. Here we will use a variant of the level set method, termed Refined Level Set Grid (RLSG) method [5, 6]. It handles all topology changes automatically, allows for easy grid convergence studies of the phase interface representation, and ensures good fluid volume conservation properties. Furthermore, the RLSG approach can provide a straightforward interface to couple the tracked representation of the phase interface during primary atomization to a Lagrangian point particle description during secondary atomization. Such a coupling is a prerequisite to simulate the atomization process as a whole and to handle the vast number of generated drops in an efficient manner. Different strategies exist to discretize the surface tension force once the location of the phase interface is known. The most commonly used method is due to Brackbill et al. [7] called Continuum Surface Force (CSF). Here, the ideally singular surface tension force is spread into a narrow band surrounding the phase interface by the use of regularized delta functions. These can take the form of a discrete derivative of a Heaviside scalar, i.e., the volume fraction, in VoF methods, or spread out delta functions, like the popular cosine approximation due to Peskin [8] in level set methods. Especially in level set methods, the use of spread out delta functions can be problematic, since convergence under grid refinement is only guaranteed for certain, not commonly employed delta function approximations [9]. An alternative to the CSF method is the Ghost Fluid Method (GFM) [10] that aims to apply the jump conditions and surface tension force as singular source terms within the context of finite difference schemes. Both the CSF and the GFM method, however, are prone to generating unphysical flows, so-called spurious currents, near the location of the phase interface when surface tension forces are present. In the canonical test cases of an equilibrium column and an equilibrium sphere, these velocity errors can grow unbounded very fast, unless they are artificially damped by introducing viscosity. The amplitude of the spurious currents when damped by viscosity is of the order of u ∼ 0.01σ/μ for VoF and level set methods and u ∼ 10−5σ/μ for marker methods [11], where σ is the surface tension coefficient and μ is the viscosity. Thus, numerical simulations are limited by a critical Laplace number, La = σρR/μ, where ρ is the density and R is a characteristic phase interface radius of curvature, since for large La, i.e., large σ, spurious currents start to dominate the physical flow [11]. The reason for the occurance of spurious currents is twofold. The major reason is a discrete imbalance between the surface tension force and the associated pressure jump across the phase interface [12]. The second source of error is due to errors associated with evaluating phase interface curvature. To address the former source of error, Young et al. [13] proposed a modification to the procedure of Kim and Choi [14] to regain discrete consistency. However, they were using the CSF method with smeared out delta functions in a level set context and, hence, the exact discrete balance was not achieved. Francois et al. [12] proposed a socalled “balanced force algorithm” for VoF schemes on structured Cartesian meshes that discretely balances the surface tension force and the associated pressure jump across the interface. In that paper, the discrete evaluation of the delta function as the derivative of the volume fraction scalar naturally results in the discrete balance when following a similar approach to the one proposed in Young et al. [13]. The approach by Francois et al. [12] eliminates spurious currents up to machine precision zero, if the interface curvature is prescribed exactly. Different strategies exist to increase the accuracy of curvature evaluation. For VoF methods, the height-function approach [15] allows second-order converging curvature calculation. However, the required stencil sizes are large and thus problematic for interfaces close to each other. For level set methods, curvature at the node location can be calculated with high-order accuracy, however, the phase interface curvature is approximated to first order at most, due to the fact that nodal location and phase interface position typically do not coincide. In this paper, we will extend the balanced force algorithm of Francois et al. [12] and Young et al. [13] to unstructured flow solver grids using the RLSG level set method to track the phase interface. To achieve second-order converging curvature evaluation, an interface projected curvature evaluation method is proposed. The performance of the resulting balanced force RLSG method is demonstrated analyzing equilibrium columns and spheres on structured and unstructured flow solver grids. Then, coupling of the RLSG method to a Lagrangian point particle representation of small scale drops is briefly discussed. Finally, to demonstrate the capability of the new method in complex flows and to analyze grid convergence behavior, a Rayleigh-Taylor instability is presented. Governing equations The equations governing the motion of an unsteady, incompressible, immiscible, two-fluid system are the Navier-Stokes equations,