A Dual-Scale Approach for Finite Weber Simulations

Direct Numerical Simulation remains an expensive task for atomization simulations. To de-crease the burden of DNS, a Dual-Scale modeling approach (Gorokhovski and Herrmann, 2008) that describes turbulent phase interface dynamics in a Large Eddy Simulation spatial ﬁltering context is proposed. Spatial Filtering of the equations of ﬂuid motion introduce several sub-ﬁlter terms that require modeling. Instead of developing individual closure models for the interface associated terms, the Dual-Scale approach uses an exact closure by explicitly ﬁltering a fully resolved realization of the phase interface. The resolved realization is maintained using a Reﬁned Local Surface Grid approach (Herrmann, 2008) employing an unsplit geometric Volume-of-Fluid method (Owkes and Desjardins, 2014). Advection of the phase interface on this DNS scale requires reconstruction of the fully resolved interface velocity. In this work results from the Dual-Scale LES model employing a sub-ﬁlter turbulent eddy reconstruction by combined approximate deconvolution and non-linear spectral enrichment (Bassenne et al. 2019) and sub-ﬁlter surface tension model (Herrmann 2013) are compared to DNS results for a phase interface in homogeneous isotropic turbulent ﬂow at two different Weber numbers.


Introduction
Atomization in turbulent environments involves a vast range of length and time scales. Predictive simulations aiming to resolve all relevant scales thus require enormous computational resources, taxing even the most powerful computers available today [1]. Since primary atomization is governed by the dynamics of the interface, a need therefore exists for appropriate interface models that make the computational cost of predicting the atomization outcome more tractable. A wide range of phenomenological models aiming to represent statistically the essential features of atomization have been proposed in the past. Although these aim to introduce the dominant mechanisms for breakup, they often use round blobs injected from the nozzle exit and hence neglect all details of the interface dynamics. Other modeling approaches to atomization include stochastic models [2,3] representing the interface by constituent stochastic particles and the mean interface density transport equation model for Reynolds-Averaged Navier-Stokes (RANS) approaches [4,5]. The former treats the interface dynamics in a stochastic sense but requires the a priori knowledge of the breakup mechanism, whereas the latter is affected by the drawbacks of the RANS approach: the transport of the mean interface density is modeled by a diffusion-like hypothesis, thereby neglecting the spatial grouping effects of liquid elements [1]. In the context of Large Eddy Simulations (LES), [6,7,8,9] have proposed models to close the unclosed terms arising from the introduction of spatial filtering into the governing equations. However, these models typically neglect the contribution of the sub-filter surface tension term and are based on a cascade process hypothesis that may be questionable in the context of surface tension-driven atomization. An exception is the model for the sub-filter surface tension term proposed in [10]. In [11,12], a Dual-Scale approach for LES of interface dynamics was proposed and a model for the sub-filter surface tension induced motion of phase interfaces was developed. The purpose of this contribution is to develop a model for the sub-filter phase interface motion induced by sub-filter turbulent velocity fluctuations. Combining such a model with the surface tension model proposed in [11,12] will result in a LES model applicable to atomizing flows.

Governing Equations
The equations governing the fully resolved motion of an unsteady, incompressible, immiscible, two-fluid system including surface tension are the Navier-Stokes equations, where u is the velocity, ρ the density, p the pressure, and µ the dynamic viscosity. Pertaining to the surface tension source term, σ is the surface tension coefficient, κ is the interfacial curvature, n the interface normal direction, δ σ is the delta function. Furthermore, the continuity equation results in a divergence-free constraint on the velocity field Assuming ρ and µ are constant within each fluid, density and viscosity can be calculated from where indices l and g denote values in liquid and gas, respectively, and ψ is the volume-of-fluid scalar that is defined ψ = 0 in the gas and ψ = 1 in the liquid. The transport equation for the volume-of-fluid scalar is defined as, Here, the last term on the right-hand side is zero for incompressible flows, see Eq. (2).

Filtered Governing Equations
Introducing spatial filtering into Eqs. (1) and (2) and assuming that the filter commutes with both the time and spatial derivatives, the filtered governing equations can be derived [7], where τ 1 , τ 2 , and τ 3 are sub-filter terms associated, with acceleration, advection, and viscous effects, respectively [7]. These sub-filter terms are defined as, Using Eqs. (3) and (4), the filtered density and viscosity in Eq. (6) are Finally, filtering of Eq. (5) leads to where the term uψ − u ψ must be modeled. The spatial filtering operator * is defined as, where G is a normalized spatial filter function and integration is over the filter volume. The surface filtering operator * s is defined as, The Dual-Scale Approach to Modeling Sub-filter Interface Dynamics Instead of relying on a cascade process by which dynamics on a sub-filter scale can be inferred from the dynamics on the resolved scale, the Dual-Scale approach proposed in [12] aims to maintain a fully resolved realization of the immiscible interface geometry at all times, expressed, for example, in terms of a volume-of-fluid scalar ψ. Then ψ can be calculated exactly by explicit filtering using Eq. (14). Although this is an exact closure, the problem of modeling is of course simply shifted to the problem of maintaining a fully resolved realization of the interface geometry, i.e., describing the fully resolved motion of the interface, Eq. (5). Since the fully resolved velocity is the sum of the filtered velocity and the sub-grid velocity, u = u + u sg , this results in where the only term requiring modeling is u sg . In [12], a model for u sg is proposed consisting of three contributions, where u is due to sub-filter turbulent eddies, u σ is due to sub-filter velocities induced by subfilter surface tension forces, and δu is attributed to the interface velocity increment due to relative sub-filter motion between the two immiscible fluids. The focus of the current contribution is on the first and second terms; for a modeling outline of the last term, the reader is referred to [12,13].

Sub-filter Turbulent Fluctuation Velocity Model
In this work we propose to reconstruct the sub-filter turbulent fluctuation velocity, u , using a Spectrally-enriched Differential Filter (SDF) model [14]. This model uses combined approximate deconvolution on an interpolated velocity field along with non-linear spectral enrichment. We first consider the coarse velocity field on grid size (2∆). The SDF-modeled velocity on the refined grid is The first term u ADI , is an approximately-deconvolved velocity field [15,16] that has been interpolated to a finer grid (∆). The second term u (∆) SGS , is the spectrally-enriched component of the sub-filter scale velocity. The purpose of this component is to add spectral energy to the sub-filter scales ranging from (∆) to (2∆) [17]. Each application of this method results in a velocity field that exists on a grid twice refined in all directions. For additional details of the SDF model, the reader is referred to [14].

Sub-filter Surface Tension Model
On the sub-filter scale, we propose to use the resolved interface geometry to obtain the subfilter surface tension velocities, u σ . The model for this velocity is that of a Taylor analogy where the surface tension forces are modeled as a spring and damper system. This model was first proposed in [11], and has been modified for this work. The model for the sub-filter surface tension velocity is where c σ = 16 and c µ = 1. For this work, c σ has been tuned such that it reproduces correct results for the damped surface waves test case [18]. It is important to note the new parameter κ which provides information about the scale of the curvature within a filter cell. The parameter κ is defined as, As pointed out by [19], for cases such as the damped surface waves test case [18] the spring component of Eq. (19) would be unable to capture acceleration of the interface. This is because it is possible that κ ≈ 0 inside of a filter cell even though the interface is curved.

Numerical Methods
Equation (16) is solved using an unsplit geometric transport scheme for volume-of-fluid scalars that ensures both discrete volume conservation of each fluid and boundedness of the volumeof-fluid scalar ψ [20]. Geometric reconstruction of the interface within each computational cell is done using PLIC reconstruction, employing analytical formulas [21] using ELVIRA estimated normals [22,23].
To efficiently solve Eq. (16) for the fully resolved immiscible interface, the RLSG method [24] is employed. By design, it solves the interface capturing advection equation on a separate, highly resolved Cartesian overset grid of mesh spacing h G , independent of the underlying LES flow solver grid of mesh spacing h. In the Dual-Scale LES approach, h G needs to be chosen sufficiently small to maintain a fully resolved realization of the phase interface.
The velocity u at RLSG scale h G is calculated from the two components discussed previously. The velocity contribution due to turbulent fluctuation is calculated by K-times application of the SDF model, with The sub-filter surface tension velocity component is computed by maintaining a cell centered, collocated velocity component to solve Eq. (19). Since κ and n are defined only near the interface, a problem arrises for u σ calculation away from the interface location. In order to overcome this, a Fast-Marching Method [25] is employed to extend the surface velocities outward to fill an entire filter cell containing an interface. Finally, the cell center, collocated velocity is calculated at the face centers using an arithmetic average. The unsplit, geometric advection scheme of [20] requires face-centered velocities that are discretely divergence-free to ensure both conservation and boundedness. While discretely divergence free filtered velocities u are available on the flow solver mesh in a standard fractional step method, i.e., ∇ h · u = 0, the RLSG velocities u are not necessarily divergence free on the fine overset mesh.
To ensure ∇ h G · u = 0, u is projected into the subspace of solenoidal velocity fields using the projection/correction step of a standard fractional step method applied to the overset mesh contained within each LES cell separately. Note that the resulting Poisson systems that need to be solved to determine the Lagrange multiplier for projection of u are defined on a per LES cell basis, necessitating an additional step to make the LES cell face average of u equal to the LES scale filtered velocity u. This is achieved by computing and applying a per cell face constant correction velocity to u, similar to the approach used to determine an outlet correction velocity in the fractional step method to make the approach satisfy the continuity equation on the entire domain. Here the correction velocity applied to each cell face ensures satisfying the continuity equation for u on the LES filter scale h. The SDF model is applied at every LES time step such that the resolved scale velocity u is actively evolving with the LES velocity u on the LES time scale. However, the sub-filter surface tension velocity evolves with the interface at the RLSG time scale. Due to this, the projection/correction of the RLSG velocities must happen throughout the RLSG time steps. During these sub-cycles, the SDF modeled velocity is frozen constant. Finally, to calculate ψ, Eq. (14) is evaluated by setting the filter size to the local flow solver mesh spacing h and evaluating the integral by explicitly summing the volume-of-fluid scalar ψ of those overset-mesh cells that are contained within a given LES flow solver cell.

Comparison Metrics
The auto-covariance spectrum of a passive scalar in a homogeneous isotropic turbulent field is a quantitative measure for how well the turbulent velocity reconstruction performs. [26] reported a range of k −1 scaling for the auto-covariance of passive scalars in turbulent flow that is followed by a k −5/3 scaling. The volume-of-fluid scalar fluctuation about a plane normal to the y-axis can be calculated by where * x,z denotes averaging with respect to the x and z directions. The auto-covariance of the volume-of-fluid scalar fluctuation about the location of the initially flat interface can be calculated as Computing the Fourier transform of R ψ (r x , r z ) using Eq. (24) leads to the volume-of-fluid scalar auto-covariance spectrum The probability density function (PDF) of interface curvature will also be used as a comparison metric. The curvature that is sampled is the mean curvature defined by, where κ 1 and κ 2 are the principle curvatures calculated by height function. Samples are taken from every point on the interface and |κ| is used to create the PDF.

Results and Discussion
In this section results from the case of a material interface in decaying homogeneous isotropic turbulence are presented. An initially flat interface is placed inside a box of fully developed isotropic turbulence. The fully developed isotropic turbulence field was graciously provided by R. Chiodi and O. Desjardins [27]. Both density and viscosity ratio are unity in a turbulent field that has a Reynolds number of Re λ = 156. The Weber numbers in this work are W e λ = ∞ and W e λ = 1.1.
Results were obtained using an LES using the Dual-Scale approach employing a LES mesh resolution of 32 3 and an overset mesh resolution of 512 3 . DNS results were obtained using a 512 3 mesh for both the flow solver and interface transport. The compensated volume-of-fluid auto-covariance spectrum, k R ψ (k), after approximately one tenth of an eddy turnover time is shown in Figure 2. A flat compensated auto-covariance spectrum denotes k −1 scaling. It can be seen that the SDF model is able to accurately reproduce the compensated auto-covariance spectrum for both Weber numbers. The results for the infinite Weber number case are noticeably worse than the W e = 1.1 case. This could be due to the absence of surface tension combined with the absence of the viscous sub-range of turbulence in the SDF velocity reconstruction. Figure 3 compares the Probabilty Density Function of interfacial curvature. It can be seen that the populations of mean curvature are reproduced very well using the Dual-Scale model. In the infinite Weber number case, there is a dropoff in probability at (3h G ) −1 due to numerical surface tension effects related to the PLIC reconstruction. The Hinze scale for the W e = 1.1 case is also accurately reproduced.

Conclusions
A Dual-Scale modeling approach for phase interface dynamics in turbulent flows is presented that is based sub-filter velocity models to describe effects from sub-filter surface tension and sub-filter turbulent velocity fluctuations. The method uses an overset high-resolution mesh to capture a resolved realization of the phase interface geometry that can be explicitly filtered to close the terms that require modeling in the filtered Navier-Stokes equations.
Comparison of DNS and Dual-Scale LES results show great agreement for both the W e = ∞ and W e = 1.1 test cases. The Dual-Scale LES model is able to reproduce the Volume-of- Fluid auto-covariance spectrum which demonstrates accurate transport of the Volume-of-Fluid scalar. The PDF of interfacial curvature is also accurately reproduced and the Hinze scale is correctly captured using the sub-filter surface tension model. Agreement in both metrics are an indicator that the proposed Dual-Scale LES model may be applicable to atomization cases where droplet generation is initiated by turbulent eddies.

Acknowledgements
The support of NASA TTT grant NNX16AB07A is gratefully acknowledged.