A Dual Scale Model for Reconstructing Sub-Filter Shear-Induced Instabilities Using a Vortex Sheet Method

A method to predict sub-ﬁlter shear-induced velocities on a liquid-gas phase interface for use in a dual scale LES model is presented. The method reconstructs the sub-ﬁlter velocity ﬁeld in the vicinity of the interface by introducing a vortex sheet at the interface. The vortex sheet is transported by an unsplit geometric volume and surface area advection scheme with a Piece-wise Linear Interface Construction (PLIC) representation of the material interface. At each step and desired location the shear-induced velocities can be calculated by integrating the vortex sheet and other relevant quantities over the liquid-gas surface with the sub-grid velocity reconstruction limited to a small number of cells near the phase interface. The vortex sheet method is tested and compared against prior literature.


Introduction
Liquid atomization is an important process occurring in many engineering applications. Internal combustion engines depend on the rapid atomization and evaporation of fuels to quickly and efficiently mix with air before combustion can occur. Evaporation of a liquid fuel is a slow process, but can be greatly enhanced by increasing the surface area of the liquid fuel. Since the residence time in combustion engines is small, the liquid fuel must be rapidly atomized into many small drops to provide a large surface area and increase the rate of evaporation. Many modern engines rely on the fuel injector design to generate the required turbulence to rapidly atomize the liquid fuel and obtain an optimal gaseous fuel/air mixture. Engine performance, efficiency and pollutant production strongly depend on the quality of the fuel/air mixture prior to combustion, and thus the details of the liquid-gas interface dynamics and atomization are of great importance. Predicting the turbulent interface dynamics remains a challenging task for numerical simulations. Direct numerical simulations (DNS) have provided great insight for studying many aspects of turbulent immiscible interfaces. However, DNS must resolve all the relevant scales of motion, requiring enormous computational resources to simulate even simple geometries [1]. This requirement severely limits the range of resolvable time and spatial scales available to DNS and restricts DNS from being a viable tool for engineering design. A need therefore exists for alternative modeling approaches to predicting turbulent interface dynamics. A number of models have been introduced to predict breakup, including stochastic models [2,3] and interface transport equations for Reynolds-averaged Navier-Stokes (RANS) equations [4]. The stochastic model requires a priori knowledge of the break-up mechanism for accurate predictions. Meanwhile the RANS approach models the mean interface density with a gradient diffusion-like hypothesis, which ignores the spatial grouping effects of liquid elements [1]. Many engineering applications of atomization, including aircraft engine combustors and diesel engine injectors can exhibit swirling flows, recirculation regions and jets in cross-flow or co-flow that are hard to predict using a RANS approach. Large-Eddy Simulations (LES) are often preferred in these applications, and therefore an atomization model consistent with the LES methodology is a desirable engineering design tool. Several LES models for turbulent immiscible interfaces have been proposed in the past [5,6]. These LES models, however, require the existence of a cascade process to predict the unresolved scales, and furthermore require the dynamics of the unresolved scales to be inferred from the dynamics of the resolved scales. The LES methodology has proven to be remarkably successful in single phase turbulent flows due to the existence of the energy cascade. However, It remains unknown whether a similar cascade process can be taken advantage of to model the dynamics of turbulent immiscible interfaces and atomization. What's more, high resolution simulations of turbulent liquid jets show that small droplets can be ripped out from large ligaments in areas of high shear, circumventing any cascading process that a traditional LES approach would use [7]. Surface tension forces tend to increase at increasingly small length scales due to the smaller local radius of curvature, similar to the viscous forces responsible for the energy cascade. Although viscosity acts to dissipate kinetic energy at small scales, surface tension can either reduce surface corrugations or amplify them via an instability mechanism like Kelvin-Helmholtz, Rayleigh-Taylor or Rayleigh-Plateau. These instabilities all rely on sub-filter interface geometry to predict sub-filter corrugation growth, and thus require knowledge of the sub-filter interface geometry. Details of the sub-filter interface geometry are unavailable in the traditional LES approach, and therefore a dual-scale approach was proposed to provide a fully resolved realization of the sub-filter interface geometry [8] and properly handle the sub-filter effects. In this work a model is presented capable of predicting the effects of sub-filter shear-driven dynamics on the resolved interface geometry.

Methods
The governing equations for the fully resolved motion of an unsteady, incompressible, immiscible, two-fluid system in the absence of surface tension and gravitational acceleration are the Navier-Stokes equations, where u is the fluid velocity, ρ is the density, p is the pressure, and µ is the dynamic viscosity. Surface tension and gravitational acceleration are neglected to focus on the shear driven instabilities of the interface. In addition to the momentum equation, the conservation of mass constrains the velocity field to be divergence-free. In the incompressible regime, fluid properties are taken to be uniform throughout each fluid. Therefore, ρ and µ are evaluated with a volume-of-fluid scalar, ψ as, where the l and g subscripts denote properties in liquid and gas respectively. The volume-offluid scalar ψ is evaluated as ψ = 0 in the gas and ψ = 1 in the liquid. In addition ψ must also be transported with the flow field as, where the last term on the right-hand-side is zero for incompressible flows due to Eq. (1).

Filtered Governing Equations
Following the methodology of LES modeling, a spatial filter is applied to Eq. (1), where the overbar ( * ) implies spatial filtering, and where τ 1 , τ 2 and τ 3 represent the sub-filter effects due to acceleration, advection and viscosity respectively [5]. Applying the spatial filter to Eq. (2) yields The spatially filtered volume-of-fluid can be evaluated by solving where τ ψ is the sub-filter liquid flux and is obtained by applying the spatial filter to Eq. (3) and making use of Eq. (1).

The Dual-Scale Approach to Modeling Sub-Filter Shear-Induced Velocities
The classical LES modeling approach in single-phase flows assume the existence of a cascading process where there is a net transfer of energy to small scales. Applying a cascade process to the atomization process for a liquid jet for example, would imply that the jet first breaks up into large scale structures and then continues to break up into increasingly small-scale structures. However as mentioned previously, evidence from high-resolution simulations of atomizing liquid turbulent jets suggest that the atomization process does not follow a cascade process. These simulations show that small-scale drops can be ejected during the ligament-formation process, circumventing any cascade process for the phase interface geometry [7]. Instead of relying on a cascade process for the sub-filter motion, the dual-scale approach aims to maintain a fully resolved realization of the interface geometry at all times [8] as shown in Fig. (1). The dynamics of this interface are governed by Eq. (3), where u is the fully resolved fluid velocity. The fully resolved velocity is decomposed into its filtered u and sub-filter u sg components, which can then be substituted in Eq. (3) as Finally, ψ in Eq. (9) can be evaluated by a direct explicit filter, where G(x) is a spatial filter function. With a model for u sg and making use of Eq. (14), there is no need to construct a model for τ ψ , and Eq. (9) can be evaluated directly. What's more the unclosed terms in the Navier-Stokes equations, τ 1 , τ 2 and τ 3 , can be calculated directly by first evaluating the fully resolved realization of ρ, µ, ψ and u using Eqs. (2), (12) and (13), and then taking an explicit filter of the terms in Eqs. (6), (7) and (8). Finally, it is worth noting that because ρ and µ are uniform throughout each fluid, the terms τ 1 and τ 3 reduce to zero when the spatial filter does not contain an interface, and τ 2 reduces the standard sub-grid stress term that can be modeled by any classical single-phase LES technique. Therefore the dual-scale procedure needs only be applied in the vicinity of the interface. The dual-scale method does present an exact closure of the sub-grid terms, however the modeling task is now shifted to maintaining a fully resolved realization of the interface geometry by solving Eq. (13) and modeling the sub-filter velocity u sg . A model for u sg is proposed consisting of four parts, where u , δu, u σ and u g are the sub-filter velocities due to turbulent fluctuations, shear-induced instabilities, surface tension and acceleration instabilities respectively. Models for u and u σ are presented in [9] and [10] respectively, a model for δu is presented here and u g is the subject of future work.

Sub-Filter Velocity due to Shear-Induced Instabilities
To model the shear-induced sub-filter velocities, we consider the motion of the phase interface between two fluids that are two-dimensional, inviscid and incompressible. This motion is governed by the Euler equations presented here in dimensionless form as where the velocity u and density ρ are defined on either side i of the interface Γ, and p is the pressure. Additionally, the boundary conditions at the interface are given by and the velocities are constrained to U ±∞ far from the interface. In these boundary conditions, reference density, velocity and length respectively, and κ is the local curvature of Γ. The evolution of the vortex sheet strength can be derived by introducing velocity potentials into the Euler equations and its accompanying boundary conditions to produce a vortex sheet transport equation [11,12,13].
The terms on the left-hand side of Eq. (19) represent the temporal changes and convective transport of the vortex sheet strength. The terms on the right-hand side describe the stretching of the vortex sheet and surface tension effects. Note that η is a surface quantity so Eq. (19) only needs to be solved at the location of the interface, so in addition to tracking the vortex sheet the interfacial area A in each cell is transported via the following equation [14].
Additionally the circulation C around the interface is defined as the line integral around the interface, taking the form of Finally the velocities due to the presence of the vortex sheet can be evaluated via a line integral of the form where the interface Γ is parameterized by the arc length of the interface s. Milne-Thomson showed that Eq. (22) can be transformed into v(x, y) = 1 2 L 0 η(s) sin (2π(x − x(s))) cosh (2π(y − y(s))) − cos (2π(x − x(s))) + µ 2 ds, for vortex sheets of length L that are periodic in the x-direction [15]. Note that the to avoid singularities in the integrals in Eqs. (22), (23) and (24) a desingularization parameter µ 2 is added to the denominator [16].

Numerical Approach
The Navier-Stokes equations are solved using NGA, a structured, staggered, finite difference flow solver with a fractional step method [17]. The task of maintaining a fully resolved realization of the phase interface geometry is achieved by solving Eq. (13) on a high resolution auxiliary Cartesian grid independent of the underlying flow solver grid. The Refined Level Set Grid (RLSG) method [18] is used to manage the auxiliary grid and activate it in regions where the spatial filter contains an interface as illustrated in Fig. (1). Eqs. (13) and (20) are advanced using an unsplit geometric transport scheme. The volumeof-fluid scalar transport is executed via the computational framework of Owkes and Desjardins directly [19], and the split transport interfacial area transport scheme of James and Lowengrub [14] is extended to the unsplit framework of Owkes and Desjardins. This method ensures that volume-of-fluid scalars remain bounded and that discrete volume is conserved for each fluid. The geometric interface within each computational cell of the RLSG is built using PLIC reconstruction with analytical formulas [20] and ELVIRA estimated normals [21]. Additionally the circulation C is regarded as a surfactant and transported geometrically with the interfacial area via the scheme of James and Lowengrub [14]. The vortex sheet strength is subsequently updated in each RLSG cell containing an interface upon completion of the advection step with the following equation.
The sub-filter velocities in Eq. (22) are numerically evaluated by first triangulating the PLIC reconstruction in each cell and then numerically integrating Eq. (22) (or Eqs. (23) and (24) in the periodic case) with a 9-point Gaussian Quadrature formula for triangles [22]. To evaluate the vortex sheet strength η at the Gaussian Quadrature points in the triangulated PLIC surfaces, it is linearly interpolated using gradients computed by neighboring cells [14].

Results
The sub-filter velocity generation technique is tested with a desingularized version of the Moore singularity problem following Krasny [16]. The initial condition for the interface is given by where A 0 and B are the amplitude and wavelength of the disturbance respectively. The initial vortex sheet strength is then given by where η is the unnormalized vortex sheet strength. For this test case the initial amplitude and wavelength are A 0 = 0.01 and B = 1, and the unnormalized vortex sheet strength is η = −1.
The initial velocity field can then be evaluated with Eqs. (23) and (24) with the aforementioned Gaussian Quadrature formula for triangles. The simulation will take place in a 1 × 1 domain with an equidistant 256 × 256 cartesian grid and periodic boundary conditions on the left and right walls. As shown in Fig. (2) the results from this study (labeled "VoF-PLIC") are in good qualitative agreement with that of the vortex sheet roll-up study conducted by Krasny [16]. The overall shape between the results is identical for early times with some slight deviation occuring late in the simulation, however the number of turns generated by the vortex sheet is the same between each set of results. Quantitatively this method can be compared to Krasny's study with the arc length of the interface after 1.0 second of simulated time at several desingularization parameters µ. The arc length of the VoF-PLIC vortex sheet method is computed by simply summing up the interfacial area A in each cell. The results of this comparison are presented in Table (1) and they show excellent agreement with this study and that of Krasny [16] with the notable exception of µ = 0.10. This is believed to be caused by the desingularization parameter µ, which acts to numerically diffuse the vortex sheet. In cases where µ is sufficiently large the numerical viscosity will dampen out high wavenumber disturbances. Since the method relies on a discontinuous PLIC geometric representation of the interface, high wavenumber disturbances can exist at cell faces and be amplified. This results in roll-up occurring not only at the low wavenumber of the domain but at high wavenumbers leading to excessive stretching of the interface and a larger arc length. Many atomization models are based on the Kelvin-Helmholtz Instability, it is therefore necessary that the VoF-PLIC model is capable of reproducing growth rates in the linear regime of the KHI. For a pair of fluids with equal density in the linear regime, the growth rate w is given by Eq. (28) [13] where k is the wavenumber of the disturbance. The initial conditions for the interface will be identical to those of the previous numerical experiment with the exception that the initial amplitude A 0 is reduced to A 0 = 1 × 10 −5 , and the experiment is only run until t = 0.5. This ensures that the growth of the instability remains within the linear regime. The growth rate of the disturbance will then be measured by taking a least squares regression of the amplitude A(t) to Eq. (29).
The results of this experiment are given in Table (2), and they show that the VoF-PLIC method is capable of matching linear theory as µ is reduced. These results also show how µ can replicate the effect of having a finite shear layer thickness by suppressing the growth rate at larger values.

Conclusions
In this paper a method to reconstruct sub-filter shear driven velocities for use in a Dual-Scale LES model has been presented. The method generates shear driven velocities by applying a vortex sheet at the interface location of a phase interface between a liquid and gas. The velocities induced by that vortex sheet can then be found by numerically integrating Eq. (22) with a Gaussian Quadrature formula. These velocities are then used to transport the interface and surface quantities with an unsplit geometric transport scheme. Finally, the updated interface geometry can be explicitly filtered and sent back to the underlying Navier-Stokes flow solver. The method has been tested against well-known results and shows excellent agreement in capturing the motion of the interface under reasonable conditions. Future work includes extending the method into 3D, including the effect of surface tension and further testing the method in a dual scale environment.