Direct Numerical Simulation of Turbulent Spray Combustion in the SpraySyn Burner: Impact of Injector Geometry

Spray combustion is one of the most important applications connected to modern combustion systems. Direct numerical simulations (DNS) of such multiphase flows are complex and computationally very challenging. Ideally, such simulations account for atomization, breakup, dispersion, evaporation, and finally ignition and combustion; phase change, heat and mass transfer should be considered as well. Considering the complexity of all those issues, and to simplify again the problem, virtually all DNS studies published up to now replaced the injector geometry by an approximated, simple configuration, mostly without any walls within the DNS domain. The impact of this simplification step is not completely clear yet. The present work aims at investigating the impact of a realistic injector geometry on flow and flame characteristics in a specific burner (called SpraySyn burner). For this purpose, two cases are directly compared: one DNS takes into account the inner geometry of the injector, including walls of finite thickness; a second one relies on a simplified description, as usually done in the literature. It has been found that considering the details of the geometry has a noticeable impact on the evaporation process and ultimately on the flame structure. This is mostly due to the effect of recirculation zones appearing behind thick injector walls; though quite small, they are sufficient to impact the evolution of the flow and of all connected processes.


Introduction
More and more researchers are aware that experiments, numerical simulations, and modeling should work hand in hand. Regarding spray combustion-the configuration considered in this work, two main workshop series concentrated on designing a proper experiment and inviting modelers for validation studies: (1) the workshops on Turbulent Combustion of Spray (TCS), mostly addressing atmospheric and low-pressure spray burners (TCS; http:// www.tcs-works hop.org/); (2) the Engine Combustion Network (ECN), where high-pressure, dense sprays are considered under conditions that are representative of engine conditions (ECN; https ://ecn.sandi a.gov/ ). These two series follow the pioneering work initiated by the international workshop on measurements and calculations of turbulent flames (TNF) (TNF; http://www.sandi a.gov/TNF/abstr act.html). In TCS, to date, three burners are considered as benchmarks: (1) Sydney's burner for turbulent piloted dilute sprays (Gounder et al. 2012), (2) CORIA spray burner in Rouen (Verdier et al. 2017), and (3) Cambridge swirl spray flame (Yuan et al. 2015). These three burners are laboratory-scale burners that are suitable to trigger model development and to support our fundamental understanding of spray combustion processes. Building on top of previous research (Kammler et al. 2001;Mädler et al. 2002;Weise et al. 2015), another coordinated initiative recently started in Germany to consider nanoparticle formation in spray flames. It will deliver an additional benchmark data based on the SpraySyn burner (Dreier and Schulz 2016;Rittler et al. 2017;Schneider et al. 2019). One of the advantages of this burner is that it was designed from the start while taking into account later numerical requirements. For example, the pilot flame in this burner is injected through a relatively wide slot, facilitating later simulations since the geometry is simpler and a very fine grid is not needed at the inlet of the pilot flame to discretize this region.. The SpraySyn burner is designed as part of a collaborative DFG project (SPP1980) entitled "Nanoparticle Synthesis in Spray Flames SpraySyn: Measurement, Simulation, Processes".
For this SpraySyn burner, most existing simulations rely on Large Eddy Simulations, like for instance (Rittler et al. 2017;Schneider et al. 2019). These two publications considered a computational domain that starts just above the injector lip. Therefore, the inner geometry of the injector is not considered. Schneider et al. (2019) divided their publication into two parts: (1) Large-Eddy simulations (LES), and (2) Experiments, obviously considering the real geometry, but without direct comparisons between both parts. In the studies of (Rittler et al. 2017;Schneider et al. 2019), the numerical models have been tested individually and work fine. The authors are involved as well in this collaborative DFG project, with the ultimate objective of providing DNS results for conditions as close as possible to those found in the SpraySyn burner. For this purpose, one of the first steps is to examine the impact of the injector geometry. This is potentially important, since any change in injection condition might modify the flow and flame characteristics, impacting the evaporation process and, ultimately, nanoparticle synthesis. To the best of the authors' knowledge, this is the first DNS study pertaining to spray combustion and taking into account the inner injector geometry of the burner. Obviously, several interesting past studies considered the impact of walls regarding gaseous or liquid injection, but with a somewhat different focus. For example, Gruber et al. (2010) used DNS to study the direct flame-wall interaction, and compared the results with the free flame structure far away from any wall. Hamzehloo and Aleiferis (2013) employed LES to investigate different injector geometries regarding hydrogen direct injection. Badock et al. (1999) and Yao et al. (2016) examined in purely experimental studies the cavitation, and the impact of different nozzle geometries on spray structure and combustion performance, varying nozzle diameter and orifice length, respectively. Lyra et al. (2015) simulated a jet in cross-flow using DNS, but without considering explicitly the real geometry used in the experiment, replaced by mathematical profiles for all variables at the inlet. None of these studies concentrated on the issue considered in the present work: by a direct comparison, quantify the differences observed in DNS between two identical cases taking into account-or not-the geometry and wall thickness of the injector.
For this comparison, DNS of spray evaporation and combustion is done in a turbulent, spatially-evolving jet (SEJ), either (1) considering the inner injector geometry, including thick walls-Case I, or (2) starting the DNS just above the injector lip-as done in most other DNS and LES studies, leading to Case II. This is a first step toward an even more faithful representation of the exact conditions used for the SpraySyn burner.
This paper is organized as follows: Sect. 2, the employed simulation domains are described. All numerical models are presented in Sect. 3. The obtained results are discussed in Sect. 4, before concluding in Sect. 5.

Employed Simulation Domains
As stated before, the ultimate purpose of this work is to carry out DNS with conditions as close as possible to that of the experimental studies conducted at the University of Duisburg-Essen (Schneider et al. 2019) with the SpraySyn burner. In current experiments, the main solvent is o-xylene (liquid), which is injected together with a dispersion gas ( O 2 ) through the central SpraySyn injector. The liquid is evaporated by a pilot flame ( CH 4 /Air) burning under lean conditions (equivalence ratio of 0.7) as shown in Fig. 1.
In the current DNS simulation, a skeletal kinetic mechanism is used to describe o-xylene oxidation. It was developed at the University of Duisburg-Essen, and consists of 24 species and 107 elementary reactions. A computational domain with dimensions of 8.192 mm (transverse direction) ×16.38 mm (streamwise direction) ×1.024 mm (spanwise-direction) is employed, which is discretized over 33.6 million equidistant grid points (resolution of 16 μm ). Taking into account recommendations from previous DNS studies of spray flames (Wang and Rutland 2007;Neophytou et al. 2011;Borghesi et al. 2013), the selected grid size is small enough to capture all relevant scales, while being still larger than the droplet diameter, as needed for non-resolved spray DNS simulations. The extension of the domain in spanwise direction is just sufficient to allow for the development of three-dimensional vortical structures, but is kept small to obtain acceptable computing times. Previous studies revealed that 2D simulations are already sufficient to get the main features of flame and nanoparticle production in a similar burner (Weise et al. 2015), so that a small spanwise extension should be sufficient for the present study as well, keeping also in mind that the turbulence level is quite low. The numerical domain is combined with inflow/outflow boundary conditions in streamwise direction, and to periodic boundary conditions in the other directions. In published DNS studies, accurate non-uniform profiles for velocity, turbulence, temperature, species are typically not available as boundary conditions. Sometimes, some measurement data or preliminary RANS or LES results might exist. But, in most cases, just the overall information in terms of flow-rates, averaged compositions and temperatures are implemented as uniform boundary conditions. This approach is the one retained as well in this comparison. Hence, the profiles are uniform in both cases at the level of the boundary conditions. It is of course expected that the agreement would become 1 3 better if the boundary conditions of Case II would be set based on (non-uniform, possible time-dependent) results from Case I at the corresponding level. However, if DNS results from Case I have been computed and are available, nobody would be interested in computing Case II any more...
In the present publication, the emphasis is set on the impact of taking into account the injector geometry in the DNS. As a consequence, all the models connected to the modeling of nanoparticles are deactivated to reduce computation requirements; the injected spray droplets do not contain any precursor and consist of a pure single-component liquid (o-xylene); the Population Balance Equations implemented to represent the formation and evolution of nanoparticles are inactivated. Only the turbulent two-phase flow, including evaporation and flame characteristics, are represented in the DNS.
In all simulations, initially liquid droplets with a diameter of 10 μm and a temperature of 300 K are injected through a planar nozzle with a diameter of d j = 0.3 mm into a domain filled with air at a temperature of 500 K.The dispersion gas ( O 2 ) is injected with jet velocity, U j , of 91.34 m/s, leading to a jet Reynolds number of 711 and a maximum Mach number is 0.3. The velocity of each droplet at the inlet is set equal to the velocity of the surrounding gas; this value is computed by multi-linear interpolation of the instantaneous gas velocity from the Eulerian (DNS) grid to that in Lagrangian space (droplet center). The gases supporting the pilot flame ( CH 4 /Air) enter the domain with a speed of 8 m/s through an area with an inner diameter of 1.244 mm and an outer diameter of 3.0 mm . Finally, the external coflow used in the experiments to insulate the process against external perturbations is injected with a low speed of 0.4 m/s through an area with an inner diameter of 3 mm and an outer diameter of 8.192 mm . The simulation has been stopped after two flow-through times based on the jet velocity and streamwise extension of the domain.

Numerical Models
All simulations are performed using the in-house DNS code DINO (Abdelsamie et al. 2016;Chi et al. 2017Chi et al. , 2018. Three essential models are activated to conduct this study: (1) The solver for the low-Mach Navier-Stokes equations for the gas phase, including detailed models for all thermo-chemical processes, (2) The Lagrangian approach for tracking the spray droplets, including evaporation and mass exchange between liquid and gaseous phase, (3) Finally, a second-order, in-house immersed boundary method (IBM) to take into account the injector geometry. A brief review of these different models is given below.
For the gas phase, the low-Mach formulation of the Navier-Stokes equations is discretized in space using a 6th-order finite-difference central stencil. The Poisson equation is solved using fast Fourier transform (FFT) for periodic as well as non-periodic boundary conditions, using in the last case an in-house array pre-and post-processing technique (Abdelsamie et al. 2016). The explicit third-order Runge-Kutta method is used for time integration in the current simulations. The mixture-averaged model is employed to compute species diffusion velocity and all transport properties using the open-source library Cantera 2.4.0. These equations are coded on top of a 2D pencil decomposition to enable large-scale parallel simulations on distributed-memory supercomputers, using the open source library 2DECOMP&FFT (Li and Laizet 2010).
The disperse (droplet) phase is tracked in a Lagrangian frame as non-resolved droplets [following the concept of discrete particles simulation, DPS Thévenin 2017, 2019)]. Hence, the resulting simulations correspond overall to the DNS-DPS approach. In this study, the droplet volume fraction is O(10 −4 ) and the average Weber number is 3.2; therefore, breakup processes can be safely ignored (Irannejad and Jaberi 2014). A two-way coupling between both phases is implemented via the exchange of mass, momentum, and energy. The droplet equations rely on the model first introduced by Abramzon and Sirignano (1989), with later improvements. This model is valid to describe the transient droplet behavior (Miller et al. 1998). It assumes infinite heat conductivity and uniform temperature inside each droplet. The approach would not be appropriate for large droplets (larger than the grid resolution). Many published DNS and LES studies used this model since it is affordable and provides a good level of accuracy (Miller et al. 1998;Neophytou et al. 2011Neophytou et al. , 2012Borghesi et al. 2013;Rittler et al. 2015Rittler et al. , 2017Thévenin 2017, 2019). The implemented equations describing droplet location, momentum, mass transfer, and heat transfer have been detailed in previous publications (Abdelsamie and Thévenin 2017, 2019), and read as follows: In Eqs.
(1)-(4), k , a k , and ∞ are the velocity, diameter of the kth droplet and velocity of the surrounding gas at droplet location k , respectively. Also, T ∞ , T k , L v , W F , C F p,f and B T,k are mixture temperature in far-field, liquid droplet surface temperature, molar latent heat of droplet vaporization, molar mass of the fuel, specific heat of the fuel vapor in the film region and heat transfer number, respectively. The properties and variables in the film region are computed based on the one-third rule (Abramzon and Sirignano 1989;Wang and Rutland 2007;Abdelsamie and Thévenin 2017) and have the subscript f, as mentioned above. In the current work, the density ratio between liquid and gas phase is large ( ∼ 1300 ), and the droplet diameter is small; therefore, the drag force is the only force that must be considered (Maxey et al. 1997;Squires and Eaton 1991), and is found in Eq. 2. Motion and evaporation of the droplets is characterized by three characteristic time scales: momentum relaxation time ( τ v,k ), evaporation delay ( τ a,k ) and heating delay ( τ T,k ) (Neophytou et al. 2011(Neophytou et al. , 2012Borghesi et al. 2013): In these equations, the characteristic time scales are computed as a function of various dimensionless numbers: the droplet Reynolds number, Re k = ∞ | ∞ − k | a k ∕ f , the Spalding mass transfer number ( B m ) and the heat transfer number ( B T ), Where f , Y s,k , Y F,∞ , W O , P ∞ and P sat,k are the gas viscosity, the vapor surface mass fraction (saturated vapor mass fraction), fuel mass fraction in far-field gas mixture, oxidizer molar mass, far-field pressure and saturated vapor pressure computed with the Clausius-Clapeyron equation: In Eq. (12), R f , P ref , and T ref are perfect gas constant, reference pressure and temperature, taken here as atmospheric pressure and boiling temperature of the fuel at this pressure, respectively, while L v is corrected using the Watson equation, Here, L v,s and T cr are the molar latent heat at temperature T ref , and critical temperature of the fuel (here, o-xylene), respectively. As shown in Eqs. (10) and (11) the heat transfer number depends on the fuel vapor to gas mixture specific heats ( C F p,f , C p,f ) at film region, Prandtl number (Pr), Schmidt number (Sc), Sherwood number (Sh), and Nusselt number (Nu) which are computed as in Borghesi et al. (2013): The third important model relies on the immersed boundary method (IBM), which has been used to take into account the injector geometry including non-zero wall thickness on the regular DNS grid. For this purpose, a novel 2nd-order IBM algorithm based on a ghostcell level-set method recently developed in our group has been implemented, as described and validated in detail in Chi et al. (2020). In this method, a local directional extrapolation scheme is employed, leading to an accurate representation of the boundaries on the DNS grid.

3 4 Results
As explained previously, the results of two different DNS will be compared directly in this section: (1) Case I, taking into account the inner details of the injection system; (2) Case II, neglecting completely the details of the injector geometry by simulating only the red rectangular box shown in Fig. 1; in the latter case, the DNS inflow boundary condition is set just above the tip of the injection nozzle. Figure 2 shows a cut plane through the center of the domain with the temperature field at four time instants t = 88 j , 93 j , 97 j , and 103 j , where j = d j ∕U j is the jet time-scale, equal to 3.28 μs here. The first time, t = 88 j , has been chosen so that the head of the starting jet has fully left the computational domain, while the following time intervals have been set arbitrarily. The initial part of the process is not included in the analysis; these four times correspond to a fully established burning jet flame. Looking at this figure, it can be noticed that the onset of turbulence is found earlier in Case I, i.e., closer to the injector. Additionally, the structure of the pilot flame is quite different. While two clearly separated, independent pilot flames appear in the cut-plane in Case II, the pilot flame is strongly bent toward the central injector in Case I, being suddenly interrupted by the colder dispersion gas containing the droplets near the axis. In spite of these significant differences concerning the temperature profile at injection, the pilot flame stabilizes again detached from the injector inlet in a similar way in both cases (look for example at the red rectangles in Fig. 2 at t = 103 j ). However, the length of this detached flame is almost 50% longer in Case II compared to Case I, and it is also slightly wider. Figure 3 shows the time evolution of velocity magnitude for the two cases. It can be seen that the velocity magnitude is noticeably larger in Case I compared to Case II, with peak values typically higher by 20%. Additionally, the turbulence is stronger. It must be kept in mind that the flow injection conditions are exactly the same in both cases regarding velocity and turbulent fluctuations. In Case I the injector walls guide the jet, leading initially to a reduced transversal expansion and, hence, to a sharper jet as it will be discussed later in connection with Fig. 7. Furthermore, the relatively thick walls of the injection system induce the formation of relatively strong vortices, increasing locally the turbulence intensity directly after injection. This explains the faster expansion of the jet in lateral direction and the higher peak velocities found farther from the injector lip, while the jet looks initially wider but evolves with a much slower expansion in Case II compared to Case I. This will be quantified later. Figure 4 depicts the corresponding time evolution of vorticity magnitude for the two DNS simulations. In Case I, stronger and larger vortices are observed compared to those found in Case II. This will be responsible for the local flame extinction events discussed in the last subsection. These strong and large vortices also explain the faster radial expansion of the central jet in Case I. Periodic vortex shedding induced in Case I by the thick injector 1 3 Fig. 2 Contour of instantaneous temperature at different time instants, from top to bottom: t = 88 j , 93 j , 97 j , 103 j . Left: Case I. Right: Case II. Solid spheres represent liquid droplets; the size of these droplets is multiplied by a factor of 5 for visualization purposes. The colour scale is the same for all figures 1 3 Fig. 3 Contour of instantaneous gas velocity magnitude at different time instants, from top to bottom: t = 88 j , 93 j , 97 j , 103 j . Left: Case I. Right: Case II. Solid spheres represent liquid droplets; the size of these droplets is multiplied by a factor of 5 for visualization purposes. The colour scale is the same for all figures 1 3 Fig. 4 Contour of instantaneous gas vorticity magnitude at different time instants, from top to bottom: t = 88 j , 93 j , 97 j , 103 j . Left: Case I. Right: Case II. Solid spheres represent liquid droplets; the size of these droplets is multiplied by a factor of 5 for visualization purposes. The colour scale is the same for all figures 1 3 walls lead to pulsations within the whole DNS domain. This periodic vortex shedding appears in Case II as well, due to the boundary conditions, but with weaker magnitude.

Scatter Plots and Conditional Averaging
In order to quantify more precisely the differences that have been discussed previously, it is helpful to examine scatter plots and to use conditional averaging, as shown in Fig. 5. It is obvious by comparing the scatter plots for Case I (Fig. 5, left) and Case 2 (Fig. 5, center) over all times that the vorticity magnitude is larger in Case I, particularly so in the intermediate temperature range of 800-1200 K. This is confirmed by looking at the conditional average presented in Fig. 5, right. At some time instants (here, t = 88 j ), the vorticity magnitude in Case I is even larger than that of Case II for the whole temperature range, as it can be seen from the conditional average in Fig. 5 Right. However, in most conditions, Case II shows a slightly larger vorticity magnitude for the high temperature range (above 1200 K).
The coupling between vorticity and temperature has of course important consequences regarding evaporation, which is the central process for nanoparticle production in spray flames. It is thus expected that the differences found between Case I and Case II will lead as well to noticeable changes regarding evaporation. In order to quantify this point, the scatter plot of the mass fraction of o-xylene Y F in the gas mixture (i.e., after evaporation) versus temperature is presented in Fig. 6. This quantity, which is a direct measure of evaporation intensity in this initial stage, shows considerably higher values for Case I compared to Case II. This is particularly clear when looking at the conditional average in Fig. 6, Right. This is attributed to the fact that higher vorticity regions enhance the evaporation rate for Case I. The main differences between the two cases are again observed in the intermediate temperature range, the one indeed corresponding to evaporation of the cold jet by the hot pilot flame. As a consequence of these differences regarding evaporation, strong differences regarding nanoparticle generation are expected as well, as discussed in more details in a separate publication ; those aspects are not examined further here.

Transverse Profiles as Function of Height Above Burner
Under real burner operation for nanoparticle production, the exact location of the extraction probe used to collect the nanoparticles is critical. A slight change in this position may lead to completely different properties in terms of particle size distribution (PSD) and amount of nanoparticles. This problem is even more complex when keeping in mind that the extraction probe will perturb the reacting flow in different ways (impacting momentum, temperature, trajectories...), so that its effect should indeed be considered as well in the simulation; this will be the subject of a dedicated publication. In order to check this point, the evolution of key quantities (velocity magnitude, gas temperature, mass fraction of o-xylene in the gas phase) are compared between Case I and Case II along transverse Fig. 7 Transverse profiles of velocity magnitude (top), temperature (center) and mass fraction of o-xylene in the gas mixture (bottom) averaged in time and over the spanwise direction at different heights above burner (from left to right, HAB/d j = 20 , 30 and 40). The black continuous lines correspond to Case I, the dashed blue lines to Case II. All horizontal axes are identical; the vertical axes are kept identical for each quantity cuts (i.e., at a constant distance from injector lip) for three different heights above burner (HAB): HAB=6 mm, 9 mm, and 12 mm, or in non-dimensional units HAB/d j = 20 , 30 and 40. The results shown in Fig. 7 have been obtained by double averaging: in space over the spanwise direction, taking advantage of the periodic boundary condition; and in time over a relevant time-window ( t = 82 j − 103 j ), in order to get rid of spurious instantaneous effects. As it has been already discussed previously (look back at Fig. 2) the velocity in Case I has been observed to be larger than that of Case II. This is found to be true over all HAB (Fig. 7, top). The differences are particularly large close to the center, due to the guiding effect initially induced by the injector walls, and to the differences in vortex structures induced further downstream. Looking now at the average temperature profiles, it can be seen that the initial temperature profile do not differ much initially. However, due to the differences observed in the vorticity magnitude and later on regarding evaporation, the results for Case I and Case II become increasingly different further downstream. Noticeably higher peak temperatures are found close to the exit of the numerical domain in Case II, confirming Fig. 2. A wider temperature profile is also observed in Case II, in line with the much delayed evaporation discussed next.
The higher velocity and vorticity observed in Case I and already discussed in connection with Figs. 2 and 3 has been found to increase the evaporation rate of the o-xylene spray. As a consequence, the mass fraction of o-xylene in the gas mixture (after evaporation) is considerably larger in Case I compared to Case II, as it is is obviously seen from Fig. 7, bottom row. The differences are observed almost immediately after injection. Already for HAB/d j = 20 , the peak values are more than twice as high in Case I than in Case II. But these differences persist in the whole simulation domain. For HAB/d j = 40 , the peak value of Case I is still more than twice that of Case II. Also in the transverse direction, o-xylene is found in the gas phase over a much wider extent. Keeping in mind that this is the first step toward nanoparticle production, very significant differences are expected as well between Case I and Case II regarding final process outcome.

Conclusions
In this work, two different numerical settings regarding flow injection have been compared, considering in the DNS simulation the injector geometry including non-zero wall thickness-or not. As a result, the impact of taking into account the inner geometry of the injection nozzle with thick walls on flame structure and evaporation process has been examined using direct numerical simulations in conditions similar to those used in the real SpraySyn burner. It has been found that including the injector walls in the DNS leads to unexpectedly large differences, first in the flow field (velocity and vorticity). This is due to the generation of strong vortical structures immediately behind the relatively thick walls of the injector. These vortices impact the flame structure and-as an indirect consequence-lead to strong changes regarding the evaporation rate of the spray. Finally, evaporation is considerably increased (factor 2 on the peak mass fractions in the gas phase) when taking into account the injector geometry. At the same time, radial mixing is enhanced as well. It is therefore strongly recommended to take into account as far as possible the inner details of the injector geometry in numerical studies. Further comparisons with experimental measurements will be carried out as soon as those will be available. The additional influence of the nanoparticles on flow and thermochemical fields will be considered in future studies.