A novel inlet boundary condition for VOF-DNS simulations of the primary breakup of prefilming airblast atomizers

Prefilming airblast atomizers are becoming widely used in state-of-the-art aero engines. The planar configuration experimentally studied by the Karlsruhe Institute of Technology (KIT) has been computationally replicated by different research groups through both Large Eddy Simulation (LES) and Direct Numerical Simulation (DNS). The present investigation reports the use of a novel inlet boundary condition for the DNS study of prefilming airblast atomization that allows accounting not only for the inflow turbulence but also for the temporal evolution of the liquid film thickness at the DNS inlet. The methodology implies the resolution of subsequent singlephase LES, two-phase flow LES and two-phase flow DNS. The former simulations are solved by using the interFoam solver from the OpenFoam toolbox, whereas the latter is computed with the PARIS open-source code. Results for a widely known operating condition show a different degree of atomization when the liquid film thickness evolution is accounted for than the one obtained through the use of a constant liquid film thickness at the domain inlet.


Introduction
In prefilming airblast atomizers, fuel is deposited on a prefilmer, onto which it flows in the shape of a film driven by a surrounding high-velocity airstream. When the film reaches the atomizer edge, it atomizes. These atomizers were born in the aviation industry [1] and are commonly used in aircraft engines alone [2] or in piloted configurations together with simplex nozzles [3]. Their atomizing features are then of key importance to the engine performance, fuel consumption and generation of emissions, thus raising interest by the research community. Different groups have performed experimental studies on annular configurations [4,5,6]. The multiscale nature of the problem and their geometrical complexity has hindered the numerical works on 3D realistic configurations [7]. In this sense, planar atomizer configurations have been widely adopted after comparison with their annular counter-part [8]. Among the developed planar experiments, the test rig adopted at KIT-ITS [9] stands out due to the wide amount of operating conditions and fuels tested through a variety of techniques at atmospheric conditions [9,10,11,12] and even at higher ambient pressures [13], constituting a large database. From these works, the liquid accumulation behind the prefilmer trailing edge is acknowledged as a key phenomenon driving atomization. Given the simplicity of the planar configuration while still being representative of the breakup taking place in annular atomizers, most numerical studies have focused on the former geometry. DNS studies include Volume of Fluid (VOF) [14,15,16] frameworks; coupled VOF-Level Set [17,18,19]; or Smoothed Particle Hydrodynamics (SPH) in 2D [20,21] and 3D [22]. Studying the full atomizer geometry through DNS is unfeasible. Hence, the studied domains are restricted to the last part of the prefilmer and the first mm downstream of its edge. Inflow boundary conditions are thus key to the simulation development and the predicted atomization features. While some authors acknowledge using constant velocity profiles for the liquid and gas phases at the inlet [19,22], Sauer et al. [14,15,16] introduced the embedded DNS concept in which they prescribe turbulent fluctuations at the inlet thanks to the use of precursor LES. Warncke et al. [23] compared two DNS using a mean gaseous velocity profile or a turbulent time-varying profile mapped from the precursor LES. They found the imposed turbulence on the gaseous phase has a major impact on breakup behaviour, predicting finer droplets. The main objective of this work is to present a numerical methodology to account not only for gaseous velocity fluctuations at the VOF-DNS inlet of prefilming airblast atomizers, but also for fluctuations in the liquid film height and velocity. This allows accounting for temporal and spatial liquid film waves on the prefilmer that could influence breakup. The impact of using this inlet boundary condition (BC) in the DNS is assessed by comparison to a DNS in which a turbulent gaseous inflow is prescribed but the liquid film height remains constant at the domain inlet.

Case of study
The geometry studied is the generic configuration from the KIT test rig [9], shown in Figure 1. An airfoil shape works as prefilmer with a length L p = 70.9 mm, span b = 50 mm and edge thickness h e = 0.23 mm, with a channel height h c = 8 mm. Fuel is injected through several holes drilled along the span, yielding an effective film length L f = 47.6 mm [13]. A single operating condition is investigated. It corresponds to the test rig operated at atmospheric conditions (T = 298 K) with a gas bulk velocityū g = 50 m/s. Shellsol D70 is introduced with a volumetric flow rateV /b = 50 mm 2 /s as the working fluid replicating kerosene. Details on the fluid properties are summarized in Table 1, corresponding to a moderate Reynolds (Re = ρū g h c µ ≈ 13500) and Weber (W e l = ρū 2 g h l σ ≈ 8). This condition has been chosen since it has been extensively studied via experiments [9,10,12,21] and simulation through different methods [14,15,16,19,23]. Thus, a large amount of data is readily available for validation.

Numerical methodology
As stated, the large cost of DNS able to account for atomization implies only a small domain can be simulated. In this case, the DNS domain (covered in another subsection) will be restricted to the last part of the prefilmer and a few mm downstream of its edge. To account for the liquid film thickness evolution at the inlet of such VOF-DNS, precursor simulations need to be performed. Figure 2 sketches the workflow of the simulations to allow imposing a realistic time-varying BC at the DNS inlet in terms of velocity and liquid shape. A single-phase channel flow LES is first performed. Velocity results from this simulation are directly mapped onto the bottom side of the DNS inlet, below the prefilmer. They are also mapped onto a two-phase flow LES that allows accounting for the film evolution on the prefilmer. Velocity and liquid volume fraction results from the latter simulation are mapped onto the top side of the DNS inlet, above the prefilmer. The DNS computed according to the previous workflow (hereinafter referred to as Time-varying h l DNS or h l = f (t)) will be compared with results from a DNS in which only the velocity field from the single-phase flow LES is mapped at the inlet BC (different LES time ranges are mapped onto both sides of the prefilmer to avoid coupling frequencies at both sides). In this latter simulation (hereinafter referred to as Constant h l DNS or h l = f (t)) the liquid film thickness at the inlet BC is set as the mean value from the two-phase flow LES. Please note that both simulations then account for a turbulent gaseous inflow, the difference only residing in the liquid film height and the initial liquid-gas interface velocity profile.

Single-phase flow LES
The first step of the methodology is to perform a one-phase channel flow LES. The standard PISO solver pisoFoam in OpenFOAM [24] was used to solve the filtered governing equations: Closure on Eq. 2 is given by the subgrid scale viscosity (ν t ), which has been modelled by the Wall Adapting Local Eddy-viscosity (WALE) model introduced by Nicoud and Ducros [25]. The computational domain is shown in Figure 3(left) together with the mesh. The domain height includes the full channel height h c , its width chosen as h c · π/2· according to Moser et al. [26]. A zonal mesh is used following the strategy by Sauer et al. [15], with a cell size grading in the y-direction leading to a y + < 1 at the walls (no wall functions are then used).
The final cell count is 20.1 M elements. As boundary conditions, the bulk velocityū g = 50 m/s is prescribed at the inlet. A no-slip condition is used at the physical walls, whereas a periodic condition is used in the spanwise boundaries. Turbulence is initially triggered in the domain using the boxTurb tool from OpenFOAM. Once a steady-state is achieved, both the inlet and outlet boundaries are converted into periodic boundary conditions. Thus, the domain acts as an infinitely long prefilmer in order for turbulence to be developed and passed on to the subsequent simulations. 2 nd order discretization schemes are chosen in time (backward) and space (Gauss linear, central differencing). Variable time-stepping is used, the maximum CFL number limited to 0.2. This results in time steps in the order of 1 × 10 −7 s. Velocity data for the subsequent simulations were sampled in the central YZ plane of the domain every 1 × 10 −6 s. To validate this simulation, the mean velocity profile (temporally and spatially averaged along z) and turbulent fluctuations are successfully compared to DNS data from Iwamoto [27] (Figure 4).

Two-phase flow LES
The two-phase flow channel LES is performed with the interFoam solver in OpenFOAM [24]. This solver makes use of a VOF method, solving the advection of the cell liquid volume fraction α l = V l /V cell and adding a source term f σ for the surface tension force to Eq. 2: The density and viscosity at a cell are linear expressions of the fluid properties weighted by α l . The computational domain for this simulation is shown in Figure 3(right) with its mesh. The domain keeps the full channel height h c , its length (26.2 mm) includes half the effective film length L f and enough distance to separate the air and fuel inlets. A zonal mesh is used with a cell size grading in the y-direction leading to y + < 0.9 at the walls (no wall functions are used) and a slight grading also in the x-direction. The final cell count of this mesh is 21.4 M. A time-varying condition with data mapped from the single-phase LES is used for the air inlet (linear interpolation is used among samples). The fuel inlet velocity is prescribed through the mass flow rate from Table 1. A no-slip condition is set at the walls, whereas a periodic condition is used in the spanwise bounds. Turbulence is modelled through WALE. 2 nd order discretization schemes are used in time and space. A fixed time step (5 × 10 −8 s) is used keeping CFL < 0.25. Velocity and α l data for the DNS were sampled in a YZ plane close to the outlet every 2 × 10 −6 s. The appearance of this simulation is shown in Figure 2, where the wavy appearance of the film can be observed. After a certain distance, the mean film thickness is stabilized along the x-direction. A temporal averaging of the film thickness, also spatially averaged among the stable range in x, yielded h l = 74.5 ± 9.1 µm, matching the experimentally observed value [9] (h l ≈ 75 µm) and thus validating the simulation results.

Two-phase flow DNS
The final DNS is carried out with the PARIS (PArallel, Robust, Interface Simulator) code [28].
Since the code also uses the VOF method, the governing equations solved are versions of Eqs. 3 and 4 discretized on a cartesian grid, with the particular use of the Continuous Surface Force (CSF) method and height functions to estimate κ to compute surface tension [29].
The DNS domain is shown in Figure 5(left) together with its dimensions. The studied Re yields a Kolmogorov length scale of η = 12 µm, implying a minimum grid spacing ∆x min = 25 µm according to Pope [30]. A cartesian grid with ∆x = 10 µm is chosen for a total of 140 M cells. The inflow BC for velocity and α l has already been described. An outflow condition is used both at the top and bottom bounds. Outflow is also chosen for the streamwise outlet, with a special convective (time-varying) treatment to reduce reflections [31]. A periodic condition is kept in the spanwise boundaries. A liquid contact angle of 60 • is set according to [22]. 2 nd order discretization schemes are used in time and space, with a modified Piecewise Linear Interface Calculation (PLIC) used for the α l advection. Variable time-stepping is used to keep CFL < 0.2, leading to time steps in the order of 1 × 10 −8 s. 4 ms of physical time have been simulated for each version of the DNS inlet, allowing the observation of 2 breakup events in each case. The first step to post-process the DNS data is to identify all continuous liquid structures. This is done scanning the whole domain looking for free surfaces according to a α l threshold. Once this is done, the intact core (including any ligament attached to the film) is extracted in order to process it separately from the droplet cloud. The result is illustrated in Figure 5(right). The projection of the core to the XZ plane allows identifying the spray contour. As far as the droplet cloud is concerned, each individual droplet is assigned its volumetric diameter: Droplet velocities are averages of the velocity components from each cell composing the droplet. Any other droplet characteristics (local Re, local W e, etc.) can be computed from these magnitudes. This algorithm is applied for the whole duration of the second breakup event computed in each simulation, sampling data every 2 × 10 −6 s.

Results and discussion
Qualitative evolution of the spray The temporal evolution of the flow for the second breakup event of the h l = f (z, t) case is shown in Figure 6 (the first breakup event is discarded as it constitutes a transient not showing the same breakup mechanism depicted for this second event). Qualitative results show the expected mechanism including liquid accumulation behind the prefilmer edge, bag breakup leading to fine droplets, ligament formation and ligament breakup leading to relatively larger droplets. This is justified considering the torn sheet break-up regime is expected according to [32] given the momentum flux ratio M = (ρ g U g )/(ρ l U l ) = 15.7. The representation of the h l = f (z, t) has been omitted for the sake of brevity. It shows a clearer wrinkling of the film in the streamwise and spanwise directions before reaching the atomizer edge, but this does not significantly influence the droplet cloud appearance from a qualitative standpoint. Based on observations with relatively thick prefilmers, Holz et al. [21] stated that primary atomization is governed by liquid accumulation at the atomizer edge, thus decoupling the breakup frequency from the frequency of the film waves. This explains the small influence observed among cases. However, this may not hold for thin edges (in the order of the film height) commonly used in aero engines, for which there is no room for important accumulation. A mean breakup length has been determined from the XZ liquid core projection as the global liquid surface divided by the span, yielding L bu = 1.6 mm in the h l = f (z, t) case and L bu = 1.9 mm in the h l = f (z, t) case, far from the experimental L bu = 3.2 mm [10] but close to the value computed by Mukundan et al. [19] with the same procedure. Future works will include the individual detection of the ligaments so as to only average ligament lengths.

Droplet cloud quantitative analysis
The post-processing technique allowed detecting 1320 and 1655 distinct droplets along the second breakup event of the h l = f (z, t) and the h l = f (z, t) DNS, respectively. The droplet size and velocity distributions are depicted in Figure 7 for both DNS, comparing them against experimental data [16]. Results in both simulations show a bimodal size distribution (a distinct peak is seen around 65 µm for h l = f (z, t) and less clearly around 60 µm for h l = f (z, t)), thus catching the distinct atomization produced by bag breakup and ligament breakup. The experiments did not allow detecting droplets with a diameter d < 20 µm, similar to the ones allowed by the computational grid. The h l = f (z, t) case exhibits a better agreement with experimental results especially in terms of velocity, although the predicted velocity distributions are flatter than the ones gathered at the experiment, with more droplets in the 0 to 5 m/s and 20 to 35 m/s ranges. This will have to be further investigated by reducing the cell size or gathering more distinct breakup events for averaging. Figure 7. Droplet size pdf (left) and velocity pdf (right) predicted by DNS and compared to experimental data [16].

Conclusions and future works
A methodology to account for the liquid film evolution in VOF-DNS of prefilming airblast atomizers through a novel inlet BC has been presented. Both the precursor simulations and the DNS have been validated against experimental data. Even though the gain in predicting the breakup mechanism and the droplet size and velocity distributions when accounting for the liquid film thickness temporal variations at the inlet has been marginal, the authors expect this workflow to show its full potential in future simulations of thinner prefilmer edges. In such configuration, the liquid accumulation phenomenon is expected to lose importance in favor of membrane sheet breakup, for which the liquid film evolution may play a driving role.