An anisotropic adaptive, Lagrange–Galerkin numerical method for spray combustion

Accepted Manuscript An anisotropic adaptive, Lagrange-Galerkin numerical method for spray combustion

Jaime Carpio, Juan Luis Prieto, Pedro Galán del Sastre

PII: DOI: Reference:

S0021-9991(19)30003-8 https://doi.org/10.1016/j.jcp.2018.12.022 YJCPH 8419

To appear in:

Journal of Computational Physics

Received date: Revised date: Accepted date:

7 July 2018 18 November 2018 23 December 2018

Please cite this article in press as: J. Carpio et al., An anisotropic adaptive, Lagrange-Galerkin numerical method for spray combustion, J. Comput. Phys. (2019), https://doi.org/10.1016/j.jcp.2018.12.022

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Highlights • • • • •

Eulerian-Lagrangian formulation for time-dependent, spray combustion problems. Lagrange-Galerkin numerical scheme in anisotropic adaptive finite element for gas phase. Lagrange-Galerkin scheme decouples Navier-Stokes from the energy/species mass fraction conservation equations. Ex-RKC tackles coupled, highly-stiff, non-linear systems with enhanced numerical stability. Simulations of canonical and real-world spray combustion configurations.

An anisotropic adaptive, Lagrange-Galerkin numerical method for spray combustion Jaime Carpio∗a , Juan Luis Prietoa , Pedro Gal´ an del Sastreb a) Departamento de Ingenier´ıa Energ´ etica. b) Departamento de Matem´ atica Aplicada a la Ingenier´ıa Industrial. E.T.S. Ingenieros Industriales. Universidad Polit´ ecnica de Madrid, 28006 Madrid, Spain.

Abstract We present in this paper a novel numerical method for the simulation of time-dependent, diluted spray combustion problems. Spray combustion is a complex numerical problem owing to the diversity of the phenomena (convection, diffusion, heat vaporization of droplets and chemical reaction) and the disparity of the scales involved. Our model considers a large number of tiny droplets (liquid phase) surrounded by a gas flow (gas phase), making up a special two-phase system where droplet interaction is neglected on account of the much larger inter-droplet distance compared to the droplet radius typically found in burners. The mathematical formulation of the problem follows a hybrid, Eulerian-Lagrangian approach: the continuous, gas phase is ruled by the conservation equations of mass, momentum, species mass fractions and energy, while the liquid phase is described by a set of ordinary differential equations (ODE) that govern the properties of each droplet along its fluid trajectory. For the resolution of the conservation equations of the gas phase, a Lagrange-Galerkin method in an anisotropic adaptive finite element framework is used. In combination with the aforementioned Lagrange-Galerkin methodology we apply a second order Explicit Runge-Kutta-Chebyshev scheme, also useful to solve the ODE systems of equations for the droplets movement. Explicit Runge-Kutta-Chebyshev preserves numerical stability while also facilitating a decoupled calculation of the unknown variables of the gas and liquid phases, a most beneficial feature from a computationally point of view. Finally, we show several examples to highlight the capabilities of our numerical algorithm in canonical and real-world configurations. Keywords: Eulerian-Lagrangian formulation, Lagrange-Galerkin scheme, finite element method, anisotropic mesh refinement, spray combustion, explicit Runge-Kutta-Chebychev scheme.

1. Introduction In our current world, liquid fuels play a prominent role in a wide array of industrial applications. From the optimal design of combustion chambers and injectors, to achieving an efficient and clean combustion, the numerical simulation and mathematical modeling of the physical processes involved in the combustion of liquid fuels proves instrumental in assessing, accurately and inexpensively, the performance of the ever more powerful industrial equipment. Since combustion is only possible in gaseous phase, when burning liquid fuels vaporization becomes imperative. To produce a fast, efficient vaporization, liquid fuels must be atomized in droplets of small size so as to favor a good heat transmission; thus, we can define a spray as a system of droplets immersed in a gaseous, continuous phase [1]. Although sprays may be produced in various ways, perhaps the most practical manner of achieving atomization is by inducing a high velocity between the liquid fuel stream and the surrounding gas, as it is the case with the pressure-swirl atomizers schematized in Figure 1. ∗ Corresponding

author Email address: [email protected] (Jaime Carpio∗a )

Preprint submitted to Journal Computational Physics

January 3, 2019

spray (vaporization and combustion) region α < 1 diffusion flame

liquid atomization region α  1 air coflow

fuel film

liquid fuel Figure 1: Schematic configuration of liquid fuel combustion with a pressure-swirl atomizer. Figure adapted from [2].

As it is explained in the work of S´ anchez et al. [2], in the system shown in Figure 1 we can define two different regions characterized by the local ratio of liquid to gas mass, α(x) = m ˙ l /m ˙ g , whose value varies along the spray jet, decreasing downstream from the injector rim. In the first zone, α(x)  1, heat and mass transport phenomena are negligible and the liquid fuel is clearly separated from the gas phase by a free surface; the injector (or atomizer) provides significantly different velocities to the liquid fuel stream and to the air stream, with the hydrodynamic instabilities initiating the breakup process. The atomization region can be mathematically described by the incompressible Navier-Stokes equations for both phases, liquid and gas, which interact through an interface where surface tension effects may be relevant, and whose accurate representation is critical to the correct understanding of this zone; see e.g. [3–5] for direct numerical simulations in the break-up zone. This region is not analyzed in the current paper; however, we plan to leverage our experience in anisotropic mesh adaptation for moving interfaces [6] to describe such atomization phenomena in future works. Rather, the focus of the present paper is on the other region, α(x) < 1, where droplet distribution becomes sufficiently diluted downstream from the atomization region, and the amount of gas entrained by the spray jet is enough to provide the heat of fuel vaporization. Then, the fuel vapor and the oxygen of the air react at high temperature, first in a partially premixed region, then in a diffusion flame that surrounds the spray jet. As it is explained in [2], two features relevant to liquid fuel combustion are: the large value of S, the mass of air needed to burn the unit mass of fuel (e.g., S = 6.5 for methane and S = 15 for dodecane); and the large ratio of liquid-to-gas densities found in combustion chambers (typically ρl /ρ0 ∼ 103 ). The first feature leads to a characteristic value α ∼ S −1 << 1 to produce vaporization and combustion in close to stoichiometric conditions. Considering the second issue together with the first one, we can conclude that −1/3 (n0 being the typical number of droplets per unit volume) is large the inter-droplet distance ld = n0 compared to the typical radius of the droplet a0 , i.e., ld  a0 . This means that each droplet vaporizes and moves with negligible direct interactions between neighboring droplets. See e.g. [7] as a reference textbook for a comprehensive, up-to-date presentation of the current level of understanding of the fluid dynamics and transport phenomena involving droplets and sprays. These previous considerations facilitate the mathematical description of this two-phase problem in the diluted droplet zone (α < 1), taking into account the disparity of spatial scales associated with this very particular category of two-phase flows. Therefore, the most natural formulation of spray problems deals with the gas phase in an Eulerian framework (i.e., v(x, t)), whereas the k = 1, ...., Nd liquid droplets, Nd the number of droplets in the domain, are tackled individually, computing their properties in a Lagrangian framework (i.e., vk (t)). In this work, we take advantage of the Eulerian-Lagrangian modeling strategy, also

2

used by other groups to tackle similar problems, such as those arising in pulverized coal combustion [9, 10], or aeronautical combustors [11–13]. In contrast to the previous Eulerian-Lagrangian formulation, the so-called “multi-continua EulerianEulerian” formulation considers the liquid phase as a continuum field, describing the droplet population in terms of the number of droplets per unit volume n(x, t), and using appropriate conservation equations for the continuous fields of velocity vd (x, t), temperature Td (x, t) and radius ad (x, t) of the droplets. See e.g. [2, 14, 15] as relevant works, from which we take not only the physical modeling, but also some canonical configurations for spray combustion which are later investigated in this paper, allowing us to perform a comparison with the literature before moving to more complex configurations. Using the same Eulerian-Eulerian approach, within a finite volume discretization referred to as “Multi-Fluid method”, we can find [16, 17] for studies on solid propellant combustion, or [18] for polydisperse reacting spray flows. Besides, we want to note that while the Eulerian-Eulerian formulation circumvents the need of following the evolution of a large number Nd of droplets, it is not adequate for addressing inertial sprays with multiple crossings of droplet trajectories, as it happens, for example, when the problem presents recirculating flow regions. To overcome this drawback, different droplet classes Nc may be introduced at every turning motion, thus defining nj (x, t), vdj (x, t), Tdj (x, t), ajd (x, t), for j = 1, ..., Nc , as shown in Figure 2. However, to our knowledge, this technique has just been applied successfully in straightforward 1D configurations so far [14].

vd1(x, t)

vd2 (x, t)

vd3 (x, t) vd3 (x, t)

vd4 (x, t)

vd2 (x, t)

Figure 2: Multi-continua Eulerian-Eulerian approach with inertial spray with multiple crossing of droplet trajectories. Dashed curves represent the flow streamlines; solid lines of different color represent different classes; thick, solid line, shows the stagnation plane.

Thus, the aim of this paper is to introduce a novel numerical method using a hybrid, Eulerian-Lagrangian formulation to solve spray combustion problems in the diluted droplet zone, α < 1. Accordingly, our numerical method follows, for the gas phase, the guidelines of our previous work [19] on time-dependent combustion problems at the low Mach limit, while we also include now the handling of liquid droplets in the combustion chamber. In [19], the gas phase is solved using a Lagrange-Galerkin formulation in which the convective terms are discretized explicitly by following, backwards in time, the characteristic curves of the convection operator (Lagrangian stage), while the diffusion-reaction terms define a highly-stiff, formally parabolic problem (Eulerian stage). In addition, now we have to solve the non-linear Nd systems of ordinary differential equations for the properties of each individual droplet, coupled to the gas phase. Nevertheless, the tracking of each droplet is a task that can be implemented straightforwardly taking advantage of the methodology already developed for the Lagrangian stage during the resolution of the gas phase, such as the optimized “search-locate” algorithm to spatial point location xk in unstructured meshes [8]. The main highlights of our numerical method are: a local anisotropic mesh refinement technique [20] for spatial discretization; an explicit Runge-Kutta-Chebyshev scheme for the resolution of the coupled, highly-stiff systems of equations [21]; and a Lagrange-Galerkin method for the treatment of convective terms [22]. All these techniques add up to provide a numerical code which is, in our opinion, both accurate and efficient, making possible running physically-relevant numerical simulations in personal computers, and opening the way for parallel computations in more complex, 3D configurations. The outline of the work is as follows: after this Introduction, section 2 is devoted to present the general formulation of the spray combustion problem and the non-dimensional form of the equations. In section 3, we outline the numerical algorithm to solve the previous equations. A series of numerical experiments is carried out in section 4 to show the capabilities of 3

the numerical algorithm, comparing the results with those found in the literature. Finally, section 5 includes the main conclusions and future lines of work. 2. General formulation for the spray problem In this section we describe the governing equations of the spray combustion problem, eventually writing them in non-dimensionalized form to better display its characteristic scales, while also producing the dimensionless numbers that rule the physics of the simulation. The mathematical model is taken from [2, 14, 15]; as a result, there is some overlap between the material presented here and that found in those works. Our rationale for including such material is based on producing a self-contained work facilitating the numerical reproducibility of the proposed examples. The gas phase is described by the conservation equations in an Eulerian framework, whereas the liquid phase consists of a set of Nd individual droplets whose properties evolve in time along the fluid trajectories xk (t), for k = 1, 2, ..., Nd , in a Lagrangian approach. 2.1. Gas phase equations For the gas phase we consider the well-known conservation equations at the low Mach number limit, which include the continuity and momentum equations, together with the conservation equations for the mass fraction of the chemical species and energy, which can be written in conservative form as: ⎧ Nd  ⎪ ∂ρ ⎪ ⎪ m ˙ k δ(x − xk ) + ∇ · (ρv) = ⎪ ⎪ ∂t ⎪ k=1 ⎪ ⎪ ⎪ ⎪ ⎪ Nd ⎪    ∂(ρv) ⎪ ⎪ (m ˙ k vk − fk ) δ(x − xk ), + ∇ · (ρvv) = −∇p + ∇ · (μ ∇v + ∇vT ) + ⎪ ⎪ ∂t ⎪ k=1 ⎪ ⎪ ⎪ ⎪ ⎨ Nd  ∂ (ρYF ) m ˙ k δ(x − xk ), + ∇ · (ρvYF ) = −∇ · jF + ω˙ F + ⎪ ∂t ⎪ k=1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ∂ (ρYi ) ⎪ ⎪ + ∇ · (ρvYi ) = −∇ · ji + ω˙ i , ∀i = F, N2 ⎪ ⎪ ∂t ⎪ ⎪ ⎪ ⎪ ⎪ I

⎪ Nd I ⎪    ∂ (ρcp T ) ∂p0 ⎪ 0 ⎪ vT ) = ∇ · (ρD c ∇T ) − h ω ˙ − c j ∇T + [m ˙ k (Lv − cp Tk ) + qk ] δ(x − xk ) + ∇ · (ρc − p T p pi i ⎩ i i ∂t ∂t i=1

i=1

k=1

(1)

The unknown variables are: v(x, t) = (v1 , ..., vd ), gas velocity (where d = 2, 3 denotes the spatial dimension of the problem); YF , mass fraction of the fuel vapor F ; Yi (x, t), mass fractions of the remaining chemical species which compose the mixture for i = 2, ..., I − 1 (with I the total number of chemical species in the problem); T (x, t), temperature of the mixture; and p(x, t), hydrodynamic pressure accounting both for the spatial pressure variations and the isotropic component of the stress tensor. In addition, the state equations and the constitutive relations allow us to define the following magnitudes in terms of the previous unknown variables: ρ, gas density; μ, dynamic viscosity; cp , specific heat at constant pressure; and DT , thermal diffusivity. Regarding the I chemical species, h0i is the specific formation enthalpy, cpi the specific heat at constant pressure, and p0 (t) is the thermodynamic pressure which is spatially uniform at low Mach number flows. The term ji represents the diffusion mass flux of the chemical species i, while ω˙ i denotes the mass of species i generated by chemical reaction per unit volume and per unit time. The summations terms appearing on the right-hand sides of (1) are the source terms associated with the presence of droplets in the flow, which are modeled by the Dirac delta distribution δ(x − xk ) located at the position of each droplet xk (t). Thus, the mass and fuel conservation equations include the mass of fuel vapor produced by each individual droplet per unit time m ˙ k . The momentum exchange m ˙ k vk − fk between the liquid droplets and gas phase is also accounted for in the momentum equation, where vk is the droplet velocity and fk is the force of the gas on each droplet. Finally, in the energy conservation equation we need to include a term representing the heat transfer from the gas phase to the liquid phase, resulting in 4

droplet heating and subsequent vaporization; thus, qk denotes the heating rate experienced by each individual droplet, Lv the specific latent heat of vaporization, and Tk stands for the droplet temperature. 2.1.1. State and constitutive equations for the gas phase The description of the gas phase requires the integration of d + I + 1 conservation equations for the d + I + 1 unknowns {v1 , ..., vd , p, T , Y1 , ..., YI−1 }. The mass fraction of the Ith species (often N2 ) is I−1 computed by YI = 1 − i=1 Yi and the remaining variables of the gas phase in (1) can be written using the following state equations and constitutive relations: • We consider that the molecular mass of the air is close to that of the inert gas in the spray stream (N2 ), so that changes in the density due to the molecular weight of the mixture are only associated with the presence of fuel vapor. Assuming also that the mixture acts as a perfect gas in an open vessel, with a time-independent thermodynamic pressure p0 , the thermal equation of state can be written as: ρ=

p0 R0 T



(1 − YF ) YF + MN2 MF

−1 ,

(2)

where R0 = 8.314 J mol−1 K−1 is the universal gas constant, and MN2 and MF denote the molecular mass of nitrogen and fuel, respectively. • We define the diffusion mass flux ji considering Fick’s law with a binary diffusion coefficient Di . This coefficient is, for all chemical species i except fuel vapor, equal to the thermal diffusion coefficient, Di = DT ; for fuel vapor, we take a constant fuel Lewis number different from unity to account for the low diffusivity of most spray fuels, so that DF = DT /LF . Thus, the mass fluxes can be written as: jF =

1 ρDT ∇YF LF

and

ji = ρDT ∇Yi

∀i = F, N2 .

• The specific heat at constant pressure is assumed to be equal for all chemical species cpi = cp = cp0 , ∀i. • A simplified power law is used to describe the temperature dependence of the dynamic viscosity and thermal diffusivity: 0.7 T μ ρDT = = . μ0 ρ0 DT 0 T0 • The chemistry of the problem is modeled with a simple, irreversible one-step overall reaction between the oxygen of the air and the fuel vapor: F + sO2 → (1 + s)P + (Q), where s is the mass of oxygen consumed and Q the amount of heat released per unit mass of fuel burnt. The reaction rate (mass of fuel consumed per unit volume per unit time) follows Arrhenius’ law: ω˙ F = BρYF

YO2 exp(−Ta /T ), YO2 A

with Ta the activation temperature, YO2 A the mass fraction of oxygen in air, and B the frequency factor. Since the previous constitutive equations only involve the mass fractions of fuel and oxygen YF , YO2 , we just have to solve the species mass fraction conservation equations for these two species. Besides, ω˙ O2 = sω˙ F in the conservation equation for the mass fraction of O2 ; while in the energy conservation 

 I ∂p0 I = 0, and i=1 h0i ω˙ i = Qω˙ F . equation, i=1 cpi ji = 0, ∂t

5

2.1.2. Dimensionless numbers for the gas phase For a better understanding of a physical process, it is convenient to formulate the problem in a nondimensional form. This also shows in an explicit manner the dependence of all variables with the main controlling parameters (dimensionless numbers). Thus, for our time-dependent spray combustion problem, we consider as reference values the following fluid properties: temperature T0 , gas density ρ0 , viscosity μ0 , thermal diffusivity DT 0 , and specific heat at constant pressure cp0 ; we also select a characteristic length L0 , a typical velocity of the gas phase U0 , and a characteristic time t0 = L0 /U0 . With these characteristic scales we can define, for the continuous phase, the Peclet Pe, and Reynolds Re numbers as: Pe =

U0 L0 ; DT 0

We also define the Prandtl number as: Pr =

Re =

ρ0 U0 L0 . μ0

μ0 , ρ0 DT 0

(3)

(4)

so that Pe=RePr, and Pr = 0.7 for air. Further, the rate of chemical time to gas inertia time defines the Damk¨ ohler number Da as: Da =

L0 /U0 L0 t0 = B exp(−Ta /Tf ) , = −1 tc B exp(Ta /Tf ) U0

(5)

where Tf is named the adiabatic flame temperature, an estimation of the maximum value of the temperature at the flame, which can be computed by:

Q YO2 A . Tf = T0 1 + cp0 T0 s 2.2. Liquid phase equations The time evolution of the unknown variables for the droplets, namely velocity vk (t), temperature Tk (t) and radius ak (t) along their trajectories xk (t), is described by k = 1, ..., Nd ordinary differential equations [2]: ⎧ dxk ⎪ ⎪ = vk , ⎪ ⎪ dt ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 4 3 dvk ⎪ ⎪ ⎪ ⎨ 3 πak ρl dt = fk , (6) ⎪ ⎪ ⎪ 4 πa3 ρ c dTk = q˙ , ⎪ k ⎪ k l l ⎪ dt ⎪ 3 ⎪ ⎪ ⎪ ⎪ ⎪ 3 ⎪ ⎪ ⎩ 4 πρl dak = −m ˙ k, 3 dt where ρl and cl are the liquid density and specific heat, respectively. When writing (6), we suppose that the temperature inside the droplets Tk is uniform, a valid assumption when the thermal conductivity of the liquid fuel is much larger than that of the gas surrounding the droplet [7]. Moreover, system (6) cannot be solved without expressions for m ˙ k , fk and q˙k , which we obtain by considering the near-field solution of the velocity and temperature around the droplet; in particular, we assume that the droplet-based Reynolds numbers are small when compared to unity, and each droplet vaporizes and moves with no direct interactions with its neighbors, as in [2, 14, 15]. Although this approximation is accurate enough under most circumstances, more complete droplet models that explicitly include a dependence on droplet-based Reynolds numbers and other additional parameters (see e.g. [7]) can be added straightforwardly to our numerical scheme without introducing significant additional hurdles.

6

2.2.1. Constitutive equations for m ˙ k , fk and q˙k Using a quasi-steady analysis of the flow field near the individual droplet the drag acting on the droplet follows approximately Stokes’ formula: fk = 6πμak (v − vk ),

(7)

while the vaporization rate and the rate of heating of the individual droplets can be respectively computed by (8) m ˙ k = 4πρDT ak λk ,

cp (T − Tk ) q˙k = 4πρDT ak − Lv λk . (9) exp(λk ) − 1 In (8)-(9), λk is the dimensionless vaporization rate, determined by the balance between the convective and diffusive heat and mass transfer due to the temperature and concentration differences between the gas and the surface of the droplet [2]:

1 − YF 1 λk = ln , LF 1 − YFs with LF the Lewis number of the fuel vapor in the mixture, and YFs the mass fraction of the fuel vapor at the surface of the droplet. To compute YFs we consider the Clausius-Clapeyron equilibrium relation in terms of the droplet temperature, as follows:

Lv MF Lv , (10) exp − YFs = Ms RF TB RF Tk and

where TB is the boiling temperature of the fuel at chamber pressure p0 ; RF = R0 /MF is the fuel gas constant; and Ms represents the molecular mass of the gas at the surface of the droplet. The computation of the latter term Ms may be simplified by considering the large differences of molecular mass between fuel vapor and N2 , as we did when writing the gas perfect law (2); that is:

−1 YFs (1 − YFs ) Ms = + . MF MN2 2.2.2. Dimensionless numbers for the liquid phase Following our previous approach with the gas phase equations, we now explicitly show the scales of droplet movement. First, we define the local ratio of liquid to gas mass at injection α = m ˙ l /m ˙ g . At injection the spray is composed of a number of spherical droplets per unit volume n0 with radius a0 and density ρl , surrounding by a gas phase with the reference properties ρ0 , μ0 , and DT 0 ; so that a representative value of α can be defined as: ρl 4 α = πa30 n0 . (11) 3 ρ0 The second dimensionless parameter is defined as a ratio between a time scale associated with the droplet movement and the characteristic time of the gas phase. For the droplets we can define two time scales: a characteristic droplet accommodation time ta , and a vaporization time tv , obtained from the momentum and radius evolution of (6), respectively: ta =

2 a20 ρl , 9 μ0

tv =

1 a20 ρl , 3 DT 0 ρ0

with

ta =

2 tv , 3Pr

(12)

where both times are quite similar, being related by the Prandtl number. A further comparison of the accommodation time with the characteristic time L0 /U0 of the gas phase provides the Stokes number St: St =

2 U0 ρl a20 ta = . L0 /U0 9 L0 μ0

(13)

For small values of the Stokes number, St  1, the droplets experience a quick vaporization and their trajectories follow closely those of the gas particles; at the opposite end, St  1, the droplets move in nearly ballistic trajectories weakly influenced by aerodynamic forces, taking a longer time to vaporize. 7

2.3. Non-dimensional equations for spray combustion After all the previous considerations, we now present in dimensionless form the governing equations for the spray combustion problem, taking as characteristic values for the gas phase: temperature T0 , density ρ0 , viscosity μ0 , thermal diffusivity DT 0 , specific heat at constant pressure cp0 , characteristic length L0 , velocity U0 and characteristic time t0 = L0 /U0 . For the liquid phase, we consider the droplet radius a0 , and the number of droplets per unit of volume n0 at the injection area. For simplicity’s sake, in the following we dispense with any additional notation (prime ’) to represent dimensionless variables. Hence, the non-dimensional equations for the gas phase can be written as: ⎧ Nd  ∂ρ ⎪ ⎪ + v · ∇ρ = −ρ∇ · v + βk T σ ak λk δ(x − xk ) ⎪ ⎪ ⎪ ∂t k=1 ⎪ ⎪ ⎪ ⎪ ⎪   ⎪ ⎪ Nd    ⎪ ∂v Pr 3 ⎪ σ T σ ⎪ ρ + (v · ∇) v = −∇p + ∇ · (T ∇v + ∇v ) + βk T ak (vk − v) λk + Pr δ(x − xk ), ⎪ ⎪ ∂t Pe 2 ⎪ k=1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎨  ∂Y Nd  1 F ρ + v · ∇YF = ∇ (T σ ∇YF ) + Daω˙ F + βk T σ (1 − YF )ak λk δ(x − xk ), ⎪ ∂t LF Pe k=1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪  ∂YO Nd  ⎪ 1 2 ⎪ ⎪ ρ + v · ∇Y ∇ (T σ ∇YO2 ) + Dasω˙ F − βk T σ YO2 ak λk δ(x − xk ), O2 = ⎪ ⎪ ∂t Pe ⎪ k=1 ⎪ ⎪ ⎪ ⎪ ⎪   ⎪ Nd ⎪  ∂T 1 1 ⎪ σ σ ⎪ + v · ∇T = ∇ · (T ∇T ) − qDaω˙ F + δ(x − xk ) βk T ak λk (Tk − T ) 1 + ⎩ ρ ∂t Pe exp(λk ) − 1 k=1 (14) with appropriate boundary and initial conditions for the unknown variables. In system of equations (14), the convection terms on the left hand side of system (1) are written in their non-conservative form to show explicitly the total or material derivative, which will be dealt with by our Lagrange-Galerkin formulation. In addition, for convenience we define the new parameter βk as the following combination of dimensionless numbers: 1 2 α βk = , 3 StPr n0 L30 Given that n0 L30 is a dimensionless number, we can relate this parameter with the number of droplets injected in the domain per unit of dimensionless time, N˙ dI , a typical parameter in Lagrangian formulation, which can be computed as: N˙ dI = n0 (vkI U0 )(SI L20 )t0 = n0 U0 L20 L0 /U0 vkI SI ⇒ N˙ dI = n0 L30 vkI SI ,

(15)

where SI L20 is the injection surface and vkI U0 is the normal component of the velocity for droplets at injection (SI and vkI being dimensionless parameters). The value of parameter N˙ dI is a measure of the accuracy for the stochastic part of the numerical method, playing a role in a particle-based method (Lagrangian formulation) not unlike h, the size of the spatial discretization, does in an Eulerian framework. Hence, since the solution for the gas phase variables show statistical oscillations due to the random positions of the droplets in the domain, increasing the number of injected particles N˙ dI should reduce the statistical deviation of the numerical, gas phase solution. With the result of (15) in the above definition of βk , we obtain the final expression: βk =

2 α vkI SI . 3 StPr N˙ dI

(16)

Remark 1. For 2D planar and 3D configurations all droplets have the same weight βk , given by (16). However, in 2D axisymmetric configurations βk we have to change its definition to expression (17), since one

8

computational droplet injected at a point xkI = (zI , rkI ) represents the set of droplets injected at a differential surface 2πrkI dr, which is proportional to rkI . βk =

2 α vkI SI rkI , 3 StPr N˙ dI r¯I

(17)

with r¯I the mean value of the radius at the injection area SI , r¯I = (rI,max + rI,min )/2. For the liquid phase, the set of ordinary differential equations for each droplet k, for k = 1, ..., Nd , takes the non-dimensional form: ⎧ dxk ⎪ ⎪ = vk , ⎪ ⎪ dt ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ dvk 1 Tσ ⎪ ⎪ ⎪ = (v − vk ), ⎪ ⎪ St a2k ⎨ dt (18)  ⎪ (T − Tk ) dTk 2 Tσ ⎪ ⎪ ⎪ = − lv , λk ⎪ ⎪ dt 3StPr a2k cl /cp0 exp(λk ) − 1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ dak 1 Tσ ⎪ ⎪ ⎩ =− λk . dt 9StPr a3k with the initial values (xkI , vkI , TkI , akI ) taken at the time tkI of droplet injection. The dimensionless parameters Pe, Pr, Da, α, St must be defined for each numerical simulation; their physical meaning and relation with the dimensional fluid variables are given in (3), (4), (5), (11), and (13), respectively. Finally, the non-dimensional density ρ, reaction rate ω˙ F , and vaporization rate λk are given in terms of the unknown non-dimensional variables as: 1 , T [1 − YF (1 − MN2 /MF )] 

YO2 T0 Ta 1 − exp − , ω˙ F = ρYF YO2 A T0 T Tf 



−1 1 − YF (xk ) T0 1 MN2 1 cp0 λk = ln lv − , , with YFs = 1 − 1 − exp − LF 1 − YFs MF RF TB Tk ρ=

(19) (20)

(21)

where the dimensionless constants LF , TB /T0 , Ta /T0 , MF /MN2 , RF /cp0 , cl /cp0 , s, heat combustion q q = Q/(cp0 T0 ), adiabatic flame temperature Tf /T0 = (1 + qYO2 A /s), and latent heat of vaporization lv = Lv /(cp0 T0 ) depend on the fuel considered in the numerical simulation. 3. Numerical method In this section we present the numerical method devised to tackle the systems of equations describing the gas phase (14) and the liquid phase (18). While the liquid phase equations constitute a highly-coupled ODE system (18) for each of the droplets, the gas phase requires the solution of a system of partial differential equations (PDE) (14) in which the thermodynamical variables of the gas are coupled with the velocity of the fluid. Moreover, the systems of equations for both phases are further coupled by source terms involving gas and droplets variables. Accordingly, for the gas phase we have d + 4 spatial unknown variables {v1 , ..., vd , p, T , YF , YO2 }(x, t), while for the droplets 2d + 2 unknowns {xk1 , .., xkd , vk1 , ..., vkd , Tk , ak }(t) arise for each droplet k = 1, ..., Nd . To solve the continuous phase, we use a Lagrange-Galerkin formulation in an anisotropic adaptive, finite element framework, following an approach similar to our previous work [19]: we divide the whole time interval of integration in time subintervals [tn , tn+1 ] with a constant time step size Δt = tn+1 − tn , defining at each time subinterval, a time-marching algorithm to solve the system of equations. For simplicity’s sake, we note 9

that any conservation equation can be written as a generic convection-reaction-diffusion problem for the scalar variable c(x, t), subject to appropriate initial and boundary conditions: ⎧ Nd ∂c ⎪ ⎪ + v · ∇c = F (c) + k=1 f (c, ck )δ(x − xk ) in Ω × (tn−1 , tn ] ⎨ ∂t n−1 (22) (x) in Ω. ⎪ c(x, tn−1 ) = c ⎪ ⎩ on ∂Ω Bc(x, tn ) = 0 In (22), the right-hand side terms are generic non-linear operators: F (c) represents the diffusion and chemical reaction terms, which depend on the unknown, vector-valued variable of the continuous phase c(x, t) = {v, T, YF , YO2 }; f (c, ck ) denotes the droplet contribution to the gas phase, which depends both on c(x, t) and on the vector-valued variables of the discrete-phase ck (t) = {xk , vk , Tk , ak }. Further, B is a boundary operator, and Ω ⊂ Rd stands for the domain with a sufficiently smooth boundary ∂Ω. In our Lagrange-Galerkin scheme, we integrate the generic equation of the gas phase (22) backwards in time from tn to tn−1 along the characteristic curves X(x, tn ; t) of the material derivative operator D/Dt = ∂/∂t + v · ∇. With the notation X(x, tn ; t) we indicate the position of the fluid particle at time t which reaches the spatial point x at time tn , being the solution of the ordinary differential equation:  dX(x, tn ; t) = v (X(x, tn ; t), t) , tn−1 ≤ t < tn (23) dt X(x, tn ; tn ) = x. Thus, with the definition of X(x, tn ; t) by (23), the partial differential equation (22) can be rewritten as: ⎧ Nd ∂c(x, t) ⎪ ⎪ = F (c(x, t)) + k=1 f (c(x, t), ck (t))δ(X(x, tn ; t) − xk ) in Ω × (tn−1 , tn ) ⎨ ∂t (24) n−1 c(x, tn−1 ) = c (X(x, tn ; tn−1 )) in Ω. ⎪ ⎪ ⎩ Bc(x, tn ) = 0 on ∂Ω where c(x, t) = c (X(x, tn ; t), t) is the solution along the characteristic curves; in addition, the solution of the original problem (22) verifies c(x, tn ) = c (X(x, tn ; tn ), tn ) = c(x, tn ) at instant of time tn , since, by definition X(x, tn ; tn ) = x. Therefore, problem (22) is formally transformed into a parabolic problem, without convective terms. Such parabolic problem (24) is then solved together with the droplet equations (18) in an Eulerian stage, with the initial condition being computed in a Lagrangian stage. This highlights the main advantage of our Lagrange-Galerkin formulation, namely, its ability to deal with the convective terms, a well-known source of numerical problems, formally replacing the velocity of the fluid v(x, t) in (24) by the information along the characteristic curves X(x, tn ; t). We now build a regular, unstructured triangulation Th = {Kj ⊂ Rd , 1 ≤ j ≤ N E}, with Kj a generic element of the total of elements N E composing the mesh of the bounded spatial domain Ω ⊂ Rd . Since we perform adaptive mesh refinement, the triangulation changes with time; nevertheless, for clarity’s sake we present our algorithm in a fixed mesh. Thus, we define, associated with the triangulation Th , a conforming finite element space Vh composed of continuous, piecewise polynomials of order equal or less than m (m = 2 for all continuous variables c(x, t) = {v, T, YF , YO2 }, except the pressure p, for which m = 1). Hence, if Nh is the number of mesh points (or nodal points) of the triangulation Th , the numerical solution at time tn belonging to Vh can be written as: ch (x, tn ) = ch (x, tn ) = cnh (x) =

Nh 

Cin ϕhi (x),

i=1 h where Cin = cnh (xi ), xi denotes the i-th mesh point of Th , and {ϕhi }N i=1 is the set of basis functions of Vh satisfying ϕhi (xj ) = δij , with δij the Kronecker delta. Next, we describe the Lagrangian and Eulerian stages of our scheme to solve numerically (24) in a finite element framework.

10

3.1. Lagrangian stage The numerical solution at time tn−1 belongs to the finite element space of the current triangulation, cn−1 (x) ∈ Vh ; however, its evaluation at the departure points (also known as “feet of the characteristic h curves”), i.e., cn−1 (X(x, tn ; tn−1 )) does not, in general, belong to the current finite element space Vh . As a h consequence, to produce a sufficiently accurate initial condition we use the Galerkin projection of the solution in the finite element space Vh as follows:   ch (x, tn−1 ) = Πh cn−1 (X(x, tn ; tn−1 )) , h with Πh a projection operator onto the space Vh . Two main difficulties then arise in the calculation of the initial condition ch (x, tn−1 ). First, we need to compute the departure points of the characteristic curves X (xi , tn ; tn−1 ) (or Xn−1 (xi ) in abbreviated form) for each mesh point xi . To do that, we solve numerically the ordinary differential equation (23) considering that the velocity can be defined linearly (second order extrapolation formula vh (x, t) = v ˜h (x, t) + O(Δt2 )), using the previously computed velocities at times tn−1 and tn−2   t − tn−1 , v ˜h (x, t) = vhn−1 (x) + vhn−1 (x) − vhn−2 (x) Δt

tn−1 < t < tn ,

(25)

using an explicit, second-order Runge-Kutta scheme, with consistent accuracy O(Δt2 ). Second, once the departure points are known, we have to build cn−1 (x) ∈ Vh as a linear combination of the basis function of h the finite element space Nh  cn−1 (x) = C¯in−1 ϕhi (x). h i=1

According to the Lagrange-Galerkin formulation (described by Pironneau [23] in one of the pioneering works of the subject) this approximated solution is computed minimizing the difference between cn−1 (X(x, tn ; tn−1 )) h and cn−1 (x) in the L2 -norm: h    n−1 2 n−1 ch (X(x, tn ; tn−1 )) − ch (x) dΩ . min Ω

The solution to this minimization problem is the Galerkin L2 -projection of cn−1 (Xn−1 (x)) onto the finite h element space Vh given by: Nh  i=1

 C¯in−1

Ω

 ϕhi (x)ϕhj (x)dΩ =

Ω

cn−1 (Xn−1 (x))ϕhj (x)dΩ. h

(26)

The left hand side of (26) leads to the mass matrix of the finite element space, which can be computed exactly by quadrature rules of, at least, order 2m. In addition, the basis functions ϕi (x) can be easily xg ) between the evaluated at the quadrature points xg of each element K, via a linear mapping xg = FK ( ˆ and the mesh triangle K. In contrast, the right hand side of (26) entails the evaluation reference triangle K of the solution cn−1 at the feet of the characteristic curves for all quadrature points Xn−1 (xg ); this is a h potentially expensive procedure from a computational point of view since, to obtain sufficient accuracy in the 2D numerical simulations presented in this work, at least npg = 16 Gaussian quadrature points are used in each mesh element. As an alternative, we propose to use a cheaper “modified Lagrange-Galerkin” method which approximately computes Xn−1 (xg ) without loss of accuracy in the space discretization. This ∗ xg ) from procedure is based on the approximation of Xn−1 (xg ) via an isoparametric transformation xg∗ = F˜K ( ˆ onto the transported m-triangle K ˜ whose nodes are the feet of the characteristic the reference element K curves of the nodal points of the element K. Details of this procedure and a comparison with the classical, semi-Lagrangian interpolation method [24] applied to gas combustion problems can be found in [19]. As an additional advantage, the parallelization of the calculation of the characteristic curves as well as the projection of ch (x, tn−1 ), for each individual scalar variable of our problem, is straightforward. 11

3.2. Eulerian stage In the previous Lagrangian stage we have produced the initial condition cn−1 (x) that allows us to solve h (24) with appropriate boundary conditions. Before proceeding with the solution of (24), we now analyze the possibility to obtaining the thermodynamical variables {T, YF , YO2 } of the gas phase and then, the unknown mechanical variables {v, p} as proposed in [19] for a gas combustion problem. One of the built-in benefits of the Lagrange-Galerkin scheme is its ability to decouple the energy and species mass fraction conservation equations from the Navier-Stokes equations in gas combustion problems. However, for spray combustion this is not true: even though neither the velocity of the gas phase nor the velocity of droplets appear in the energy and species mass fraction conservation equations (24), they are needed to compute xk (t) by means of vk (t), which in turn requires the gas velocity v(xk ) from (18). To alleviate this problem while keeping the accuracy and stability of the numerical method, we use as gas velocity the explicit, linear extrapolation v ˜h (x, t) from (25). Once all gas phase variables {T, YF , YO2 }(x, tn ) and droplet variables {xk , vk , Tk , ak }(tn ) for k = 1, ..., Nd , are computed at instant of time tn , we then solve the Navier-Stokes equations to calculate the remaining unknowns {vn (x), pn (x)}. 3.2.1. Energy and species mass fraction conservation gas equations and droplet equations The aim now is to obtain, at time tn , the spatial, vector-valued variable of the continuous phase c(x, t) = {T, YF , YO2 }, and the set of vector-value variables of the discrete phase ck (t) = {xk , vk , Tk , ak }, ∀k = 1, ..., Nd . For that purpose, we have to solve a PDE system formed by the energy and mass species conservation equations (14) written in the form (24), along with Nd ODE systems (18) for the discrete phase, as follows: Continuous (gas) phase: ⎧ Nd ∂c(x, t) ⎪ ⎪ = F(c(x, t)) + f (c(x, t), ck (t))δk (X(x, tn ; t) − xk (t)) ⎨ k=1 ∂t n−1

c(x, tn−1 ) = c ⎪ ⎪ ⎩ B(cn (x)) = 0.

in Ω × (tn−1 , tn ],

(X(x, tn ; tn−1 ))

Discrete (liquid) phase: ⎧ ⎨ dck (t) = η (c(x (t), t), c (t)) k k k dt ⎩ . ck (tn−1 ) = cn−1 k

in (tn−1 , tn ],

in Ω × (tn−1 , tn ], (27)

(28)

In (27)-(28), F(c(x, t)), f (c(x, t), ck (t)) and ηk (c(xk (t), t), ck (t)) are generic vector-valued functions that may depend on the set of variables c = {T , Y F , Y O2 } and ck = {xk , vk , Tk , ak } in a non-linear way. To solve the continuous (gas) phase, we consider the weak formulation of (27): we multiply the equations by the discrete 3 3 vectorial test function ϕh ∈ (Vh0 ) and integrate in the domain Ω, looking for the solution ch (x, t) ∈ (Vh ) that verifies:

∂ch (x, t) , ϕh ∂t

Ω

= (F(ch (x, t)), ϕh )Ω +

Nd 

3

f (ch (xk (t), t), ck (t))ϕh (xk (t)) , in Ω×(tn−1 , tn ] ∀ϕh ∈ (Vh0 ) .

k=1

(29) In (29), we introduce the notation xk (t) used in the terms f (ch (xk (t), t), ck (t)) and ϕh (xk (t)) to denote the spatial point where the Dirac delta distribution provides non-zero contribution, X(xk (t), tn ; t) = xk (t). That is, xk (t) represents the gas fluid particle at time tn that was located at the position xk (t) at a previous instant of time t, with tn−1 ≤ t ≤ tn ; or, equivalently, xk (t) = X(xk (t), t; tn ). A sketch of this transformation is included in Figure 3, where we can observe that, at time t = tn , both spatial positions overlap, xk (tn ) = xk (tn ). Consequently, to evaluate ch (xk (t), t) and ϕh (xk (t)), the space point xk (t) is obtained integrating along a forward characteristic curve, using the same second order Runge-Kutta scheme employed in the calculation of the feet of the characteristic curves, with the explicit, (linear) extrapolated velocity v ˜h (x, t) from (25).

12

xk (tn−1) = X(xk (tn−1), tn−1; tn)

xk (tn−1) xk (t) = X(xk (t), t; tn) xk (tn) = X(xk (tn), tn; tn) xk (tn)

xk (t)

Figure 3: Sketch of the transformation xk (t) = X(xk (t), t; tn ).

Moreover, in the discrete (liquid) phase equations (28), we can replace ch (xk (t), t) by ch (X(xk (t), tn ; t), t) = ch (xk (t), t) to write the following ODE system: dck (t) = ηk (ch (xk (t), t), ck (t)) dt

in (tn−1 , tn ].

(30)

Equations (29)-(30) are highly coupled with non-linear relations. In particular, we remark the strong dependence on the temperature T of the chemical reaction rate of the gas conservation equations, and of the droplet source terms involving droplet radius ak and temperature Tk . The resulting, extremely stiff system of equations is hardly manageable using standard, explicit algorithms, due to the restriction in the time step size Δt as predicted by the Courant-Friedrichs-Lewy (CFL) condition. To address this issue while still keeping an explicit approach so as to decouple the discrete droplet equations from the continuous equations, we make use of the the Explicit, Runge-Kutta-Chebyshev (Ex-RKC) scheme presented in [19]. This scheme not only provides second-order temporal accuracy O(Δt2 ), but also displays a stability region that can be enlarged by adding, in a dynamic way, stages to the explicit scheme. Thus, the boundary of the real stability region has a quadratic dependence on the number of stages r. Though initially used for pure parabolic equations [25], we successfully tackled the conservation equations arising in combustion problems in [21] and [19] using a Lagrange-Galerkin formulation. Now, we propose its use in spray combustion problems, where the equations for the discrete and continuous phases must be solved simultaneously. 3 The Ex-RKC scheme with r stages used to calculate the numerical solutions cnh (x) ∈ (Vh ) of (29), and 3 ¯n−1 cnk ∈ R2d+2 , for k = 1, ..., Nd of (30), taking as initial conditions at instant of time tn−1 , c (x) ∈ (Vh ) , h

13

and cn−1 ∈ R2d+2 , respectively, can be written as: k ⎧ 0 ¯n−1 Wh = c (x), ⎪ h ⎪ ⎪ 0 ⎪ n−1 ⎪ wk = ck , ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ t0 = tn−1 , ⎪ ⎪ ⎪ ⎪ ⎪ compute x0k = X(xk0 , t0 ; tn ) ⎪ ⎪ ⎪ 1     0  ⎪ ⎪ ⎪ 1 Δt b0 , ϕh Ω , ⎪ Wh , ϕh Ω = Wh , ϕh Ω + μ ⎪   ⎪ ⎪ ⎪wk1 = wk0 + μ 1 Δtηk Wh0 (x0k ), wk0 , ⎪ ⎪ ⎪ ⎪ ⎪ for j = 2 to r ⎪ ⎪ ⎪ ⎨ tj−1 = tn−1 + αj−1 Δt, ⎪ ⎪ = X(xkj−1 , tj−1 ; tn ), compute xj−1 ⎪ k ⎪ ⎪ 





⎪   1 ⎪ j j−1 j−2 ⎪ ⎪ W W W , ϕ = (1 − μ − ν ) W , ϕ + μ , ϕ + ν , ϕ j j j j h h h h h ⎪ h h h Ω ⎪ Ω Ω Ω ⎪    j−1  0 ⎪ ⎪ ⎪ + μ  Δt b , ϕ + γ  Δt b , ϕ , ⎪ j j h h Ω Ω ⎪ ⎪

 ⎪     ⎪ j j−1 j−2 0 ⎪ ⎪ j Δtηk Whj−1 (xj−1 ), wkj−1 + γ wk = (1 − μj − νj )wk + μj wk + νj wk + μ j Δtηk Wh0 x0k , wk0 , ⎪ k ⎪ ⎪ ⎪ ⎪ ⎪ ⎪end ⎪ ⎪ ⎪ ⎪ cnh (x) = Whr , ⎪ ⎪ ⎩ ck (tn ) = wkr , (31)   3 where Wh0 , . . . , Whr ∈ (Vh ) and wh0 , . . . , whr ∈ R2d+2 are internal variables, and bj , ϕh Ω represents the right-hand side term of (29): Nd  

   j f Whj (xjk ), wkj ϕh (xjk ). b , ϕh Ω = F(Whj ), ϕh + Ω

k=1

 We note that, when the velocity of the gas phase shows up in the evaluation of the terms ηk Whj (xjk ), wkj , we use vjh (xjk ) = v ˜h (xkj , tj ). Besides, the coefficients appearing in (31) are given analytically by: ⎧ ε T  (w0 ) 2 ⎪ ⎪ , w0 = 1 + 2 , w1 = r , = ⎪ ⎪ 13 r Tr (w0 ) ⎪ ⎪ ⎪ ⎪ ⎪ Tj (w0 ) ⎪ ⎪ ⎪ b = 2 (2 ≤ j ≤ r), b0 = b1 = b2 , aj = 1 − bj Tj (w0 ),   ⎪ ⎨ j Tj (w0 ) ⎪ 2bj w0 bj 2bj w1 ⎪ ⎪ , νj = − , μ j = , γ j = −aj−1 μ j , μ 1 = b1 w1 , μj = ⎪ ⎪ bj−1 bj−2 bj−1 ⎪ ⎪ ⎪ ⎪ ⎪ Tj (w0 ) ⎪ α2 α2 ⎪ ⎪ α (2 ≤ j ≤ r), α1 =  = = w , α0 = 0, 1  ⎩ j Tj (w0 ) T2 (w0 ) 4w0

(32)

with Tj (x) the Chebyshev polynomial of the first kind and degree j, 2 ≤ j ≤ r. Finally, since the Ex-RKC algorithm (31) leads to an uncoupled, well-conditioned linear system of equations, one for each continuous variable and droplet property, a parallel implementation of this algorithm is straightforward. 3.2.2. Temporal discretization of the Navier-Stokes equations In addition to decoupling the energy and species mass fraction conservation equations from the NavierStokes equations, by taking a Lagrange-Galerkin approach we can transform the non-linear Navier-Stokes

14

equations into a linear Stokes problem, as follows:  Nd   ∂v Pr  σ  3 T σ ρ = −∇p+ ∇· T ∇v + ∇v + βk T ak (vk −v) λk + Pr δ(X(x, tn ; t)−xk ) ∂t Pe 2

in

Ω×(tn−1 , tn ].

k=1

(33) The continuity equation for spray combustion problems can be written, at time tn , as: d D log(ρ) 1  + βk T σ ak λk δ(x − xk ) Dt ρ

N

∇·v=−

in

Ω.

(34)

k=1

All the thermodynamical properties of the gas {Thn , YFn h , YOn2 h }, and those of the droplets {xkn , vkn , Tkn , ank }, along with the associated variables ρn and λnk , are known from the previous Ex-RKC algorithm. Therefore, the only unknowns in (33) and (34) are the velocity vh (x, t) and the (hydrodynamic) pressure ph (x, t). For the temporal discretization of the momentum equation (33), we use a BDF, second-order formula [26]; for the spatial discretization we employ Taylor-Hood finite elements, that is, quadratic m = 2 for the velocity, d vh ∈ (Vh ) , and linear m = 1 for the pressure, ph ∈ Qh , which satisfy the Ladyzenskaja-Babuska-Brezzi (LBB) condition. Thus, the weak formulation of the Stokes problem, with appropriate boundary and initial conditions, can be stated as: ⎧



n−2  n n Pr n σ  n ⎪ n 3vh n n T n 4vh − vh ⎪ ⎪ ρ , wh (T ) ∇vh + (∇vh ) , ∇wh , wh − (ph , ∇ · wh )Ω + = ρ ⎪ ⎪ 2Δt Pe h 2Δt ⎪ Ω Ω Ω ⎪ ⎪ ⎪  Nd ⎪  ⎨ 3 σ d + βk (T n (xk )) ank (vkn − v ˜hn (xk )) λnk + Pr wh (xk ), ∀wh ∈ (Vh0 ) , 2 ⎪ k=1 ⎪ ⎪ ⎪    ⎪ Nd ⎪  ⎪ D log(ρ) 1  ⎪ σ n  ⎪ , qh + n βk (Thn (xk )) ank λnk qh (xk ), ∀qh ∈ Qh , ⎪ ⎩ (∇ · vh , qh )Ω = − Dt tn ρ k=1 Ω (35) n−2 n−1 with vn−2 (x) calculated from v (X(x, t ; t )) using the same algorithm employed for v (x), with n n−2 h h h the caveat that in this case, the feet of the characteristic curves are computed two time steps backwards in time. In the saddle-point problem (35), the unknown variables are written on the left-hand side, whereas the right-hand side terms act as source terms computed from already available variables. In addition, we have used the extrapolated velocity v ˜hn (xk ) in the contribution of the droplets to the momentum conservation equations similarly as we did with the discrete, droplet momentum equations solved with the Ex-RKC scheme. Finally, the linear systems of equations arising from the spatial discretization are solved by an Uzawapreconditioned, conjugate gradient algorithm developed by Dean and Glowinski [27]. Since most of the computational time is consumed in the resolution of this Stokes problem, for future applications and 3D simulations we will take advantage of the PETSc (‘Portable, Extensible Toolkit for Scientific Computation’) tool [28], using the PCFieldsplit approach based on the Schur-complement, block-preconditioning of the Stokes matrix system as we proposed in [29] for multiphase flows. 3.3. Mesh refinement The numerical algorithm previously described can be embedded into a local, mesh adaptation procedure. Adaptive mesh refinement (AMR) is specially useful in flows showing a large disparity of scales, producing accurate results at a reduced computational cost. Moreover, if the problems display directional features such as the slender jets, thin flame fronts, or mixing layers typical of combustion problems, anisotropic mesh adaptation provides excellent accuracy with outstanding computational savings when compared to uniform or isotropic refinement, turning it into a tremendously appealing technique for combustion problems [30, 31]. Our aim with the local adaptation procedure is to adapt the computational domain to the particularities of the solutions as the integration time goes on. For that propose, we first compute the numerical solution of the problem at instant of time tn ; then, we compute a certain error indicator η n of the error incurred in the 15

numerical resolution of the equations from tn−1 to tn on the spatial triangulation Tnh ; finally, if the error is larger than a given tolerance η n  T ol, we build an optimal mesh Tn∗ h with the smallest number of elements satisfying the tolerance. For anisotropic mesh adaptation, we require the size, shape and orientation of the elements, as shown in Figure 4. This information is then collected into the so-called “metric tensor”, which v1 r2,K

λ2,K

r1,K

v2

v3 λ1,K

Figure 4: Parameters that define an anisotropic triangular element in 2D.

can be computed at each element K of the spatial triangulation Tnh as: −1 T MK = |K|−2/d RK SK RK ,

(36)

where the scalar |K| defines the size of the element; the matrix RK = {r1,K , ..., rd,K } with ri,K · rj,K = δij , identifies the principal directions of element K (orientation); and the stretching factors matrix SK = diag{s1,K , ...., sd,K } characterizes the shape of the element K, defined from the length of the semi-axis of

 −2/d d the associated d-ellipsoid, by si,K = λ2i,K λ for i = 1, ..., d. While there are several strategies j,K j=1 to define the metric tensor (see e.g. [32, 33]), we follow our previous works [19, 20]. Specifically, to compute the metric tensor we perform an analysis of the error incurred in the Lagrangian stage for the resolution of the gas phase, when we approximate the scalar variable cn−1 (X(x, tn ; tn−1 )) by cn−1 (x); this approach has h h proved useful in the solution of convection-dominated problems, as it is usually the case in fluid dynamic simulations. The error incurred in the Lagrangian stage can be evaluated in the L2 −norm by: 

1/2  n−1 2 n−1 n η = ch (X(x, tn ; tn−1 )) − ch (x) dΩ , Ω

being computed element-wise by: n

η =

 K

n 2 (ηK )

1/2

 , with

n ηK

= K



cn−1 (X(x, tn ; tn−1 )) h

−

2 cn−1 (x) h

1/2 dΩ

.

From this a posteriori error indicator, we can define the optimal size |K ∗ | of the element K so that the error indicator η n be below a given tolerance T ol in a triangulation Tn∗ h with as few elements as possible, according to: ⎞1/(2α) ⎛ 2 T ol n −2/(2α+1) ⎠ |K ∗ | = |K| ⎝  (ηK ) , (37) n )2/(2α+1) (η n K K∈T h

with |K| the current size of the element, and α = (m + 1)/d (m denotes the degree of the polynomials in (x); in our case, m = 2). We want to note, that expression (37) can be the finite element space of cn−1 h also written in terms of a certain degree of spatial resolution hmin (the minimum size of elements in the triangulation Tn∗ h ) instead of the tolerance T ol, as: |K ∗ | =

−2/(2α+1)

n ) |K| (ηK " # hdmin . −2/(2α+1) n minK∈Tnh |K| (ηK )

16

(38)

∗ ∗ From the analysis in [20], we can also define the optimal orientation RK and shape SK of the element K, n 2 using the a priori error between a function c ∈ H (K) and its linear interpolation. Here, the components r∗i,K and s∗i,K of the matrices defined above can be computed from the eigenvector-eigenvalue pair {li,K , gi,K } of the Hessian matrix HK (cnh ), according to:

s∗i,K =

$ d j=1

1/d    −1  gd+1−j,K  ,

|gj,K |

r∗i,K = ld+1−j,K ,

for i = 1, ..., d.

(39)

We then evaluate the Hessian matrix HK (cnh ) using the weak formulation of each derivative, as detailed in [34]. From (37) or (38) and (39), we define the optimal metric tensor (36), which is then passed on to an anisotropic mesh generator, such as BAMG [35] for 2D triangulations, or Mmg3d [36], PRAgMaTIc [37] for 3D meshes. Lastly, we want to remark the multivariate character of combustion problems. As a result, during the local mesh adaptation process we select the set of N unknown variables of our problem, φ = {β1 c1 (x, t), ..., βN cN (x, t)}, whose contribution to the global error is balanced using βi as the normalized weights associated with variable ci (x, t). Thus, for each scalar variable of the set φ, we define a metric tensor according to (36). The resulting metric tensor defining the optimal triangulation is formed as the intersection ∗ 1∗ N∗ = MK ∩ ... ∩ MK . Additional details on of the individual metric tensors of the scalar variables; that is, MK the metric intersection procedure can be found in [19]. 4. Numerical examples Though quite a few simulations on spray combustion at the low Mach limit can be found in the literature, it is difficult to establish a common, benchmark problem due to the variety of the constitutive equations of droplet source terms, turbulent motion and physical modeling of the combustion process. Nevertheless, in this section we try and test our numerical algorithm with a set of canonical problems taken from relevant works [2, 14, 15] in the laminar spray combustion field. The physical model in those works is quite similar to that presented in Section 2, with some differences that are pointed out when discussing each particular configuration; further, their numerical scheme is based on a Finite Difference implementation using an Eulerian-Eulerian (or multi-continuum) approach, as mentioned in the Introduction. Since real-world spray combustion is a very complex problem involving multiphase, turbulent reacting flows and largely different length and time scales, it is common, in the Fluid Mechanics community, to simplify the mathematical description of the flow by considering the characteristic scales that separate the dominant physical process in each region of the domain, thus examining several “canonical problems”. For instance, among the several canonical configurations arising in real-world, pressure-swirl atomizer we focus here in two of them: counterflow spray diffusion flames, and self-ignition of spray mixing layers, as depicted in the Figure 5. The first test we consider in this work is the counterflow configuration, which allows for an in-depth analysis of spray diffusion flames. This canonical problem produces local flow conditions with strain mixing layers where we can investigate the effects of strain rate and droplet inertia as measured by the Stokes number (see [14] and references therein). In a way, this problem is a local representation of a real turbulent flow in which the flames appear embedded in thin mixing layers that are locally distorted and strained by turbulent motion. We also note that, in applications involving spray combustion, flame-flow interactions depend on the presence of fuel droplets [7]. Next, we study the self-ignition of spray mixing layers: the self-ignition process leads to the onset of combustion downstream from the injector and upstream of the diffusion flame, combustion being facilitated by the high temperature of the surrounding gas. In many of the continuous-combustion systems, such as gas turbines or industrial furnaces, forced ignition (by electric sparks, torches,...) is not needed during normal steady operation. In these cases, injection velocities are much higher than the characteristic deflagration speeds, thereby precluding upstream triple-flame propagation. Thus, combustion stabilization must rely on the self-ignition of the fuel-air mixture caused by the high temperature of the surrounding gas. A canonical problem that helps us to understand this complex spray-ignition phenomenon consists of a mixing layer configuration of hot-air stream with a spray carried by an inert gas (see [15] and references therein). 17

spray( vaporization and combustion) region α < 1 diffusion flame

self-ignition of spray mixing layers

air

st

m rea

re y st spra

am

liquid atomization region α  1 fuel film y ra

re

am

ea

sp

st

m

air coflow

counterflow spray diffusion flames

ai

rs

tr

liquid fuel

Figure 5: Laminar canonical problems found in a real-world, pressure-swirl atomizer: counterflow spray diffusion flames and self-ignition in spray mixing layers. Figure adapted from [2].

Finally, to highlight the flexibility and full capabilities of our numerical algorithm we tackle a more complex problem with an arrangement similar to the pressure-swirl atomizer of Figure 5. We consider the axisymmetric, outer flow of swirling fuel injectors, in which a spray carried by a swirling inert gas is injected in a quiescent atmosphere of hot air. In this kind of devices, the swirling motion of the fuel stream could generate a hollow-cone spray jet. In all these examples we look for steady or stationary solutions, so that the initial conditions play no role; in any case, all numerical simulations start with “cold conditions”, with neither combustion nor droplets inside the domain. For further insight into the transient problem, see the provided supplementary material to this paper, with video clips showing the time evolution of the computed solutions along with the anisotropically-adapted triangulation. Regarding computational details, all simulations were run in a workstation with an Intel i7-3770k CPU and 16GB of DDR3 RAM@1600MHz, using the C-programming language and the GCC-4.7.3 GNU compiler with -O3 optimization and no explicit parallelization. As a rough estimation of the CPU cost of each simulation, we note that around 100 time steps per hour can be computed when the mesh consists of 10.000 elements and the number of droplets inside the computational domain is 250.000. The distribution of the total 100% CPU time among the different parts of the full algorithm is as follows: Lagrange-Galerkin scheme (4%), Ex-RKC scheme (35%), and linear Stokes problem with Uzawa-preconditioned, conjugate gradient algorithm (61%). As already mentioned in the previous section, CPU times for the Lagrange-Galerkin and Ex-RKC schemes can be reduced straightforwardly with parallel computing (e.g. OpenMP or MPI) achieving good scaling, while the most expensive part, the solution of the Stokes problem which is independent of the number of droplets, can be improved taking advantage of PETSc [28] as in [29]. 4.1. Test I: spray combustion in counterflow configuration In this example we analyze spray combustion in axisymmetric counterflow burners like that of Figure 6: the sketch represents a typical setup used in experimental studies, which involves two aligned opposed pipes of equal radius R placed at opposite ends and separated a distance 2H. From the nozzles, emerge two streams against each other: one stream of hot air flows from the bottom upwards at temperature TA , uniform axial velocity UA and oxygen mass fraction YO2 A = 0.233, while the other stream, formed by a monodisperse fuel spray carried by an inert stream of nitrogen, flows downwards from the top at temperature Ts and uniform axial velocity Us .

18

2R

Ts , Us Nitrogen + spray

z

RPe−1/2 ∼ δm

r

2H

Mixing layer

Diffusion flame

Air TA , UA Figure 6: Schematic configuration used in this test I, to characterized counterflow spray diffusion flames.

When the Peclet number Pe associated with the two gaseous, axisymmetric, coaxial counterflowing jets is moderately large, the resulting flow is inviscid in the first approximation. Viscous effects are confined to thin mixing layers of characteristic thickness δm ∼ RPe−1/2 . Of particular interest for combustion applications is the mixing layer developing at the inter-jet separating surface from the central stagnation point. As explained in [38], the solution at distances from the stagnation point that are much smaller than R is self-similar, both inside and outside the mixing layer. The self-similar solution involves a radial velocity increasing linearly with r according to vr = A(z)r/2, where A(z) is the strain rate, a function of the distance to the stagnation plane z. The accompanying axial velocity, temperature and species composition are a function of z only, as can be seen also in Figure 6. This one-dimensional counterflow solution has served as a basis for studies of flame-flow interactions in non-premixed and premixed combustion. Of relevance in connection with our study is the counterflow spray-flame analysis of Li˜ n´an et al. [14]. Their study analyzed the self-similar solution obtained by a one-dimensional Eulerian-Eulerian numerical code. Results [14] are compared with those obtained with our numerical code near the stagnation point in two-dimensional simulations for the opposed-pipe configuration of Figure 6. The non-dimensional variables are obtained using as reference fluid properties those of the spray stream of nitrogen: density ρ0 = ρs , viscosity μ0 = μs , thermal diffusion DT 0 = DT s and temperature T0 = Ts . As reference length we consider the nozzle radius L0 = R, while we take as characteristic velocity that of the spray stream U0 = Us , so that we can define a characteristic time t0 = R/Us . As fuel for this test, we use dodecane C12 H16 at atmospheric pressure and temperature Ts = 300K; all fuel properties are taken from [14], and are collected in Table 1.

Dodecane

cl /cp

MF /MN2

RF /cp

LF

TB /Ts

Ta /Ts

lv = Lv /(cp Ts )

q = Q/(cp Ts )

s

1.84

6.07

0.047

3.6

1.63

20

1.005

123.6

3.5

Table 1: Dodecane properties for Test I.

At the top nozzle we have the spray carried by a nitrogen stream, with temperature Ts < TB and axial velocity Us . Droplets are uniformly injected in a random fashion at zI = zmax along the nozzle section, in mechanical and thermal equilibrium with the nitrogen stream. At the bottom nozzle we consider hot air > TB > Ts , with an axial velocity UA so that the momentum of the spray stream is at temperature TA % balanced, UA = Us ρs /ρA . Additionally, the interior, lateral walls of the nozzle are considered slip-flow 19

boundaries to render all integration results independent of the axial dimension. Besides, the flow outside of the nozzles has no relevance for the analysis of the zone between the nozzles where the one-dimensional approximation is valid. Thus, we consider a coflow of cold air at temperature Ts with axial velocity 0.5Us , so as to convect the Kelvin-Helmholtz instabilities arising in the shear layer out of the domain. Finally, the no-slip flow condition is imposed at the lateral, exterior walls of the nozzle. Focusing now on the numerical approach used to solve this configuration, we consider a 2D axisymmetric, non-dimensional domain of radius rmax = 5 where we use an axial coordinate z ∈ [zmin , zmax ], with zmax = zI = −zmin = 3, so that the nozzle outlets are located at z = ±1. For time integration we consider a time interval t ∈ [0, tend ], with tend = 20 taken large enough to produce the self-ignition of the mixture along with a steady diffusion flame in the region between the nozzles. Numerical simulations are carried out with a constant, dimensionless time step size Δt = 0.01 that could temporarily be decreased when ignition appears; the number of injected droplets is N˙ dI = 10000 so that N˙ dI Δt = 100 new droplets enter the domain at each time step; for the spatial discretization, we use a dynamic, anisotropically refined domain with maximum √ mesh size hmax = 0.3 and minimum mesh size hmin = 0.5/ Pe. The error used in the mesh refinement procedure is weighted by a factor β(r) < 1 when r > 1, so as not to refine the shear layer produced between the inside and outside counterflow streams, which has no influence on the relevant, central zone between the nozzles. As a summary of the boundary conditions for the gas phase, we can write: z = zmax , r ≤ 1, vz + 1 = vr = YF = YO2 = T − 1 = 0, √ z = zmin , r ≤ 1, vz − 2 = vr = YF = YO2 − 0.233 = T − 2 = 0, ∂YO2 ∂T ∂YF = = = 0, 1 < |z| < zmax , r = 1− , vr = ∂r ∂r ∂r ∂YO2 ∂T ∂YF = = = 0, 1 < |z| < zmax , r = 1+ , vz = vr = ∂r ∂r ∂r |z| = zmax , r > 1, |vz | − 0.5 [1 − exp(10(1 − r))] = vr = YF = YO2 − 0.233 = T − 1 = 0,  Pr  μ ∇v + ∇vT · n = ∇YF · n = ∇YO2 · n = ∇T · n = 0, zmin < z < zmax , r = rmax . p·n+ Pe where n is the unit normal vector pointing outwards from the computational domain. In addition, the droplets are injected at zk = zI = zmax and rk < 1 uniformly, with initial conditions: vkz + 1 = ukr = Tk − 1 = ak − 1 = 0. Finally, the following dimensionless parameters define each numerical simulation: Peclet number Pe, Damk¨ohler number Da, liquid mass-loading ratio α, Stokes number St. Depending on the value of the Stokes number St, we analyze next two different regimes of spray combustion as in [14]. As shown by constant density models [2], for St < 1/4 (non-inertial spray) the droplets follow closely the gas particles, reaching the stagnation plane with almost zero axial velocity; on the other hand, when St > 1/4 (inertial spray) the droplets cross the stagnation plane and move into the opposing air stream. In contrast to [14], where an infinitely fast chemical reaction or “Bruke-Shumann limit” Da → ∞ is assumed for simplicity’s sake, we consider in this work a finite value of the Damk¨ ohler number. However, the main difference between our model and the one-dimensional counterflow model is the definition of the dimensionless time as the inverse of the axial strain spray-side A−1 s ; this quantity, as shown in Carpio et al. [38] for the two nozzle configuration of Figure 6, is close to our characteristic time, As ∼ Us /R, when Peclet is large enough and H/R = 1. Evolution of the ignition process and representation of the relevant fluid magnitudes may be found in the corresponding video clips included as supplementary material to this paper. 4.1.1. Non-inertial Sprays For small values of the Stokes number St  1 the droplets reach the gas stagnation plane with vanishing axial velocity. In this zone, they may vaporize due to the heat received from the hot air, producing fuel vapor that can burn with the oxygen in a diffusion flame embedded in the mixing layer. For this simulation we consider the set of parameters Pe = 1000, St = 0.2 and α = 0.1. 20

& ' and different values of the Damk¨ ohler number Da = 103 , 104 , 105 to assess the eventual convergence toward the infinitely fast chemical reaction model. Figure 7 shows the solution at the final instant of time tend , when the properties at the mixing layer have reached a steady-state solution: in the left panel, we plot the computational mesh and the temperature in color; further, pink dots represent droplets (only 10% of the total, for clarity’s sake) in a uniform distribution over the whole domain. The steady solution of this test displays N E ∼ 5750 mesh elements and Nd ∼ 43000 droplets scattered over the domain. The two right-most panels of Figure 7 represent two close-ups of the mixing zone: the top-most panel shows the spatially refined mesh with all 100% of the computational droplets, highlighting how the number of droplets found within a given computational element is large enough (typically, 20-50) to neglect statistical oscillation of the steady solution. The bottom-most panel shows a color map of the temperature, with dashed (white) lines representing the streamlines of the gas phase, and solid (pink) lines depicting the trajectories around the mixing layer and the diffusion flame. Due to the small inertia of the droplets St = 0.2 close to the stagnation plane, they follow the streamlines of the fluid without crossing the mixing layer, vaporizing then before reaching the stagnation plane, owing to the high amount of heat transferred from the diffusion flame. It may also be noticed from the bottom-right panel of Figure 7 the negligible dependence of the temperature on the radial coordinate for r < 1, a typical behaviour of self-similar structures used to validate the one-dimensional approximation of counterflow problems.

Figure 7: Numerical simulation of the counterflow problem, with Pe = 1000, St = 0.2, α = 0.1, and Da = 104 . The temperature is represented as a scalar color map. The insets show the mesh and droplet distribution (top-rightmost panel), whereas temperature (in color), streamlines (dashed lines) and droplet trajectories (solid line) are shown in the bottom-right panel.

Thanks to such negligible dependence of the temperature on the radial coordinate, we can represent profiles of the solution as a function only of the axial distance z for r < 1, taking the mean value of 11 cross sections from r = 0.25 to r = 0.75 to further decrease the statistical noise. Results are shown in Figure 8. The radius of the droplet ak remains almost constant until the droplets penetrate into the mixing layer. After that, vaporization occurs very rapidly, producing the fuel vapor that burns with the oxygen of the air in a diffusion flame located close to the stagnation plane marked by the vertical dashed line. The diffusion flame is found in a thin layer where fuel vapor and oxygen coexist. In the limit of high Damk¨ ohler number Da → ∞, oxygen and fuel cannot coexist and the temperature shows a peak with a discontinuous slope at the position of the flame, a behavior also found in [14]. In our case, the curves are rounded-off in the flame zone for finite values of the Damk¨ohler number; notice the'insets of Figure 8, which show profiles for & different values of the Damk¨ohler numbers Da = 103 , 104 , 105 . As the Damk¨ohler number increases, so does the stiffness of the system of equations due to the reaction terms, along with the computational cost; moreover, the zone where fuel vapor and oxygen coexist in the diffusion flame grows increasingly narrow and temperature shows a sharper peak as the Damk¨ ohler number increases. Even as we approach the infinitely fast chemical solution, our code remains numerically stable, refining the reaction zone and consequently 21

increasing the number of stages of the RKC scheme. We also note the good agreement between the profiles shown in Figure 8 and those presented in Figure 3 of [14]: the main characteristics of flame structure (shape, maximum values of temperature Tmax = 6.4 and fuel vapor YF max = 0.17) shown in [14] are quite similar to ours, although that work just considers a one-dimensional counterflow model with an infinitely fast chemical reaction limit. YO2 /YO2 A

1

ak

0.1

0.8

0.7

0.08

0.65

0.06

0.6

0.6

0.04

T/(10Ts )

0.02 0

0.4

-0.1

0.55 0.5 0.45

-0.05

0

0.05

0.4 -0.1

-0.05

0

0.05

0.2

YF 0 -1

-0.8

-0.6

-0.4

-0.2

0

0.2

z 0.4

0.6

0.8

1

Figure 8: Profiles of T , YF , YO2 , and ak for the counterflow test with non-inertial droplets St = 0.2, and Pe = 1000, α = 0.1 and Da = 104 . The insets compare the results obtained with different Damk¨ ohler numbers, Da = 103 (dashed curves), Da = 104 (thin solid curves) and Da = 105 (thick solid curves).

Before proceeding with the description of the inertial sprays, it could be useful to analyze the accuracy and sensitivity of our numerical code to the discretization parameters: the minimum size of the elements hmin , and the number of injected droplets in the domain per unit of dimensionless time N˙ dI . For the case of non-inertial sprays with a time step size of integration Δt = 0.01, and dimensionless flow parameters Pe = 1000, St = 0.2, α = 0.1 and Da = 104 , we evaluate the maximum value of temperature Tmax and fuel vapor mass fraction YF max in the central region at r = 0.5. These values as a function of time are ergodic stochastic processes, so that their statistical properties can be computed from a single, sufficiently long, random sample of both processes. We consider a simulation with a length of 5000 time steps during the stationary flame solution. Table 2 collects the mean value and the standard deviation (in parenthesis) of √ Tmax and YF max as a function of hmin Pe and N˙ dI . In Table 2 we analyze two different cases: in each case √ hmin Pe 0.5 0.5 0.5 0.5 8 4 2 1 0.5

N˙ dI 400 1600 6400 25600 25600 25600 25600 25600 25600

T max (σT ) 6.000 (0.261) 5.893 (0.149) 5.897 (0.074) 5.866 (0.042) 5.745 (0.095) 5.851 (0.065) 5.846 (0.051) 5.854 (0.043) 5.866 (0.042)

Y F max (σYF ) 0.2170 (0.0767) 0.1854 (0.0350) 0.1838 (0.0171) 0.1831 (0.0087) 0.1589 (0.0176) 0.1738 (0.0120) 0.1807 (0.0105) 0.1838 (0.0089) 0.1831 (0.0087)

Table 2: Mean value and the standard deviation (in parenthesis) of Tmax and YF max as a function of the discretization √ parameters hmin Pe, and N˙ dI for counterflow problem, with a time step size of integration Δt = 0.01, and flow parameters Pe = 1000, St = 0.2, α = 0.1 and Da = 104 .

we keep one discretization parameter constant and vary the another one to evaluate the convergence of the 22

statistical properties of the fluid magnitudes Tmax and YF max . For small value of hmin = 0.5Pe−1/2 , as we change the number of particles N˙ dI we observe how the standard deviation for both Tmax and YF max is −1/2 proportional to N˙ dI , as could be expected from a particle method. If we now focus on the sensitivity of the results to the spatial discretization parameter hmin , we observe a convergence of the mean value and a reduction of the standard deviation for both Tmax and YF max when hmin decreases, with N˙ dI = 25600 being constant. 4.1.2. Inertial Sprays For sufficiently large Stokes numbers, St  1, the droplets cross the stagnation plane and move into the opposing air stream, reaching large distances |z|  δm /R before they turn around. The vaporization of these crossing droplets, as well as the combustion of the so-generated fuel vapor, take place in the hot air stream, placing the flame far away from the stagnation plane and from the mixing layer. For the simulation of inertial sprays we consider the following set of dimensionless parameters: Pe = 1000,

St = 2.0 and α = 0.05, & ' along with different values of the Damk¨ ohler number Da = 102 , 103 , 104 . Figure 9 shows the steady solution reached in the case of inertial droplets with Damk¨ ohler number Da = 104 . As before, we present in the left panel of the figure the computational mesh and the temperature in color. Pink dots represent droplets (10% of the total). In this computation we have reached N E ∼ 6500 mesh elements and Nd ∼ 35500 droplets. In the right-most panels of Figure 9 we show two details of the flame zone between nozzles; the top-most panel shows the spatially refined mesh with 100% of the computational droplets, whereas the bottom-most panel plots the temperature (in color), dashed (white) lines for the streamlines of the gas phase, and solid (pink) lines representing the trajectories around mixing layer and the diffusion flame. Due to the large inertia, St = 2.0, droplet vaporization occurs outside of the mixing layer on the air side of the counterflow, in a zone of high temperature created by the flame. Thus, the droplets reach the flame and disappear soon after turning around.

Figure 9: Numerical simulation of the counterflow problem with Pe = 1000, St = 2.0, α = 0.05, and Da = 104 . The temperature is represented as a scalar color map. The insets show the mesh and droplet distribution (top-rightmost panel), whereas temperature (in color), streamlines (dashed lines) and droplet trajectories (solid line) are shown in the bottom-right panel.

In the Eulerian-Eulerian formulation [14], some modifications are required to accommodate for the presence of droplets with a reverse motion: at the same location (rk , zk ), we might find droplets with different values of velocity vk , temperature Tk , and radius ak . Hence, to account for advancing and returning droplets, the Eulerian-Eulerian formulation considers different classes of droplets, adding a new class to the spray description when the droplets reach null axial velocity and reverse their motion. As an advantage of our Eulerian-Lagrangian formulation, no special modification for turning droplets is needed. 23

YO2 /YO2 A

1

ak

0.1

0.8

0.9

0.08 0.06

0.6

T/(10Ts )

0.04

0.7

0.02 0.6

0

0.4

0.8

-0.2

-0.15

-0.1

-0.05

0.5 -0.2

-0.15

-0.1

-0.05

0.2

YF 0 -1

-0.8

-0.6

-0.4

-0.2

0

0.2

z 0.4

0.6

0.8

1

Figure 10: Profiles of T , YF , YO2 , and ak for the counterflow test with inertial droplets St = 2.0, and Pe = 1000, α = 0.05 and Da = 104 . The insets compare the results obtained with different Damk¨ ohler numbers, Da = 102 (dashed curves), Da = 103 (thin solid curves) and Da = 104 (thick solid curves).

Lastly, we present in Figure 10 the profiles of the solutions in the central zone (r < 1), as a function of the axial distance z. Here, we can see how the droplet radius ak starts to decrease once the droplet crosses the stagnation plane (marked by the vertical dashed line), first decreasing slowly due to vaporization, then with an abrupt drop in radius just after the droplet turns around. The diffusion flame (zone of maximum temperature, where fuel vapor and oxygen coexist) is located at the position where the droplet turns around, far away from the stagnation plane. Insets ohler & ' of Figure 10 show computations with different values of the Damk¨ numbers Da = 102 , 103 , 104 ; we can observe how the fuel vapor and oxygen coexist in an increasingly narrow area, and the flame temperature profile grows sharper as the Damk¨ohler number increases. Again, we can compare our results in Figure 10 with those of Figure 6 in [14] computed with the one-dimensional, Eulerian-Eulerian scheme. The shape and maximum values of temperature Tmax = 8.6 and fuel vapor mass fraction YF max = 0.33 found in [14] compare favorably with our own results. 4.2. Test II: self-ignition of fuel sprays in mixing layers Studies of laminar, mixing-layer configurations have been found to be instrumental in developing the understanding of the ignition of turbulent gaseous mixtures (see e.g. [15] and references therein). For the spray problem at hand, we consider the configuration shown in Figure 11, a coflow mixing layer formed by a stream of hot air (YO2 A = 0.233) moving at velocity UA with temperature TA , and a monodisperse spray moving at velocity Us < UA with temperature Ts < TA . A typical triple flame structure appears at the ignition point x = xign . As a reference for fluid properties we take those of the hot air stream, with density ρ0 = ρA , viscosity μ0 = μA , thermal diffusion DT 0 = DT A , temperature T0 = TA and a characteristic velocity at the injection U0 = UA . Since there is no characteristic length in this example, we follow [15] and consider the vaporization time to define an appropriate length scale L0 = U0 tv . With this assumption, and using expressions (12) and (13), we can define the Stokes number for this test as: St =

U0 2 2 ta tv = 0.9524, = = L0 /U0 3Pr L0 3Pr

where we have uses a Prandtl number Pr = 0.7. The liquid fuel used in this example is heptane C7 H16 , whose properties are taken from [15] and are evaluated at atmospheric pressure and temperature TA = 1000K, see Table 3. The remaining physical parameters in the governing equations are α = 1, and Da = 2.5 × 103 (this value differs from [15], because the Damk¨ohler definition (5) considers Tf instead of TA as in [15]). Moreover, the 24

UA , TA δ ∼ Pex −1/2 x

y

Lean premixed flame Diffusion flame

Hot air x

Rich premixed flame

Nitrogen + spray Us , Ts

xign

Figure 11: Schematic configuration used in test II: two-dimensional, laminar mixing layer. Black dots indicate fuel droplets.

Heptane

cl /cp

MF /MN2

RF /cp

LF

TB /TA

Ta /TA

lv = Lv /(cp TA )

q = Q/(cp TA )

s

2.2

3.57

0.083

2.6

0.37

10

0.34

39.5

3.5

Table 3: Heptane properties for Test II.

Peclet number Pe is different for each simulation. The inert stream of nitrogen which carries the fuel spray is injected with a velocity Us = 0.8UA at boiling temperature TB . The fuel spray droplets are injected at x = 0 with a uniform spatial distribution along ymin < y < 0, being in mechanical and thermal equilibrium with the inert gas. Consequently, droplet velocity at the injection is Us , and its temperature is to be computed imposing a null rate of droplet heating qk = 0 as per (9), with the expressions for the vaporization rate λ being given by (21); thus, for the heptane parameters of Table 3, we have Tk = 0.963TB at the injection of the droplets. Test II is solved in a 2D planar configuration with a non-dimensional spatial integration domain [0, xmax ]× [ymin , ymax ], with time interval t ∈ [0, tend ]. In all examples of this second test we consider tend √ = 40, xmax = 10,√while the y-coordinate is defined as a function of the Peclet number, ymin = −25/ Pe and ymax = 50/ Pe. As a summary of the boundary conditions for the gas phase, we can write: vx − 1 = vy = YF = YO2 − 0.233 = T − 1 = 0, vx − 0.8 = vy = YF = YO2 = T − 0.37 = 0, ∂YO2 ∂T ∂YF = = = 0, vx − 0.8 = vy = ∂y ∂y ∂y  Pr  μ ∇v + ∇vT · n = ∇YF · n = ∇YO2 · n = ∇T · n = 0, p·n+ Pe  Pr  p·n+ μ ∇v + ∇vT · n = ∇YF · n = ∇YO2 · n = ∇T · n = 0, Pe

x = 0,

y > 0,

x = 0,

y < 0,

x > 0,

y = ymin ,

x > 0,

y = ymax ,

x = xmax ,

with the following initial conditions for the droplets at the injection xk = 0, ymin < yk < 0: vkx − 0.8 = vky = Tk − 0.3563 = ak − 1 = 0. As mathematical parameters for the numerical integration, we consider a constant dimensionless time step size Δt = 0.02; for the spatial discretization, we use √ a dynamic, anisotropically refined domain with √ maximum and minimum element mesh sizes hmax = 10.0/ Pe, hmin = 0.2/ Pe, respectively; the number of injected droplets is N˙ dI = 5000, so that about N˙ dI Δt=100 new droplets enter the domain at each time 25

step. To obtain an accurate spatial resolution of the flame located at the mixing layer, we inject a different number of droplets N˙ dIj on our configuration depending on the injection surface SIj : √ for −1.5/√Pe < y < 0, N˙ dI1 = 0.50N˙ dI √ for −3.5/√ Pe < y < −1.5/√ Pe, N˙ dI2 = 0.15N˙ dI for −10/√Pe < y < −3.5/ Pe, √ N˙ dI3 = 0.15N˙ dI for −25/ Pe = ymin < y < −10/ Pe. N˙ dI4 = 0.2N˙ dI Since the spray should be uniformly distributed at x = 0 and y < 0, we then compute weights βkj , j = 1, ..., 4 by (16), which enter into the governing gas equations (14). Figure 12 shows, for Peclet number Pe= 50, the spatially adapted mesh for the steady solution, along with the droplet distribution. In this simulation we use N E ∼ 4000 mesh elements with roughly Nd = 50000 droplets in the entire computational domain; these values typically grow with the Peclet number. As can be observed in the insets of Figure 12, the number of droplets found within a given computational element is large enough (typically, 20-50) to neglect statistical oscillation of the steady solution, ensuring that the solution is independent of the specific value N˙ dI chosen for the simulation.

y 0.2

y

6

0 4 -0.2

-0.4

2

4

4.2

4.4

x

-1.5

0

y -2

-2

0

2

4

6

8

x

-2.5

10

6

6.5

x

7

Figure 12: Anisotropically adapted mesh obtained during the solution of test problem II for Pe=50. The insets show a snapshot of the droplet distribution at two specific locations.

We now finish our analysis of this test by comparing our results with those of [15]. In that work, the authors use a boundary layer approximation to study the structure of the ignition of laminar mixing layers in the limit of high laminar Reynolds number. Moreover, they consider a simplified model for the mass rate of vaporization λ given by the Spalding expression [15], and valid for fuels whose latent heat of vaporization is much larger than the fuel thermal energy Lv  RF TB (for heptane, Lv /(RF TB ) = 11.02). In Figure 13 we represent, for three different Peclet numbers Pe=50, 100, 200, the temperature of the flame (top-most panels), while the bottom-most panels show the reaction-rate isocontours ω˙ F given by (20), along with the droplet distribution (20% of the computational droplets). To highlight the convergence with the Peclet number, we also plot the solution at the bottom-most panels using mixing layer y-coordinates, yP e1/2 ; for values of the Peclet number above Pe> 250, the mixing layer becomes unstable as oscillations develop both in the mixing layer and in the flame proper. In particular, if we focus on Figure 3, panel a) of [15], and compare those results with our own in the last panel (Pe= 200) of Figure 13, we observe a striking 26

Figure 13: Evolution of the solution increasing Pe number from left to right Pe=[50, 100, 200]. The color isolines show dimensionless reaction-rate contours ω˙ F = Da−1 [0.01, 0.02, 0.05, 0.1, 0.2, 0.25]. Dashed lines show streamlines, solid black line the stoichiometric line.

resemblance. That is all the more remarkable since the boundary layer simplification used in [15] does not allow for upstream propagation of the information; such effect is quite relevant to this configuration, as thermal expansion decelerates the flow upstream, and the premixed flame has a propagation velocity that moves the flame upstream towards the unburnt region. Nevertheless, the resulting ignition distance, identified by the local maximum of the reaction rate, is just a little bit smaller in our simulations, xign = 4.67, than that obtained in [15] (xign = 4.95). 4.3. Test III: spray combustion in swirling fuel injectors In this test, we consider a more complex problem, the axisymmetric, outer flow of swirling fuel injectors, whose sketch can be seen in Figure 5. At the initial instant of time, a spray jet carried by an inert (nitrogen) gas, both at temperature Ts , starts to discharge impulsively (impinging jet) into a stagnant region of hot air (YO2 A = 0.233 and TA = 2Ts ) through a circular orifice of radius R. The axial discharge velocity is assumed to be uniform with a value Us small enough to render all compressibility effects negligible (low Mach number limit). Moreover, we consider cases with a Peclet number Pe moderately large, but small enough that the resulting jet remains laminar and axisymmetric. In an effort to accurately replicate the conditions found in actual engineering applications, both the inert gas and the droplets are discharged in the camera with azimuthal velocities vθ and vkθ , respectively, like a swirl injector. Using both linear momentum J and angular momentum L of the injection jet we define the swirl ratio of the jet at z = 0 as S = L/(RJ). If the axial discharge velocity is Us uniform, and the azimuthal velocities of both the droplet and the inert gas are equal to a rigid-body rotation with angular velocity ωs , the swirl ratio results S = ωs Us /(2R). The azimuthal momentum equation for the gas phase and for the droplets determine both azimuthal velocities vθ and vkθ , respectively. However, these equations depend in a non-linear manner on the radial velocities of the gas and the droplets, through the terms vθ vr and vkθ vkr . To overcome this drawback, we define as unknown variables the non-dimensional circulations Γ and Γk , instead of using vθ and vkθ , respectively, according to: Γ(x, t) =

1 r vθ (x, t) , S R Us

and

Γk (t) =

1 rk vkθ (t) . S R Us

To define the non-dimensional variables relevant to this example, we make use of the properties of the nitrogen gas stream at temperature Ts ; further, we take as characteristic length and velocity, L0 = R and U0 = Us , respectively. With the previous definitions, and assuming an axisymmetric configuration, the

27

azimuthal momentum equations to compute the non-dimensional circulations can be written as: ⎧

N  c Γ Pr 3 DΓ ⎪ 2 σ σ ⎪ = ∇ · r Pr δ(x − xk ), T ∇ β T (x )a (Γ − Γ(x )) λ + + ρ ⎪ k k k k k k ⎪ ⎨ Dt Pe r2 2 k=1 ⎪ ⎪ dΓk 1 Tσ ⎪ ⎪ = (Γ − Γk ). ⎩ dt St a2k

(40)

System of equations (40) must be solved together with the energy and species mass fraction conservation equations for the gas (14), and the droplet equations (18), as presented in Section 2.3. As the non-dimensional circulations of (40) are independent of the radial and axial velocities, we can compute them later by writing the standard equations in axisymmetric form v = (vz , vr ), and introducing a force term in the radial component induced by the swirl motion: ⎧  2 Nc    Pr 3 Dv ⎪ σ T σ 2Γ ⎪ = −∇p + ∇ · (T Pr δ(x − x ∇v + ∇v ) + β T a (v − v) λ + ) + ρS er , ρ ⎪ k k k k k ⎪ ⎨ Dt Pe 2 r3 k=1 ⎪ ⎪ Γ2 dvk 1 Tσ ⎪ ⎪ ⎩ = (v − vk ) + S 2 3k er . 2 dt St ak rk

(41) The swirl motion induces a radial velocity in the flow as well as a significant flow deceleration of the axial component due to mass conservation. If the swirl ratio S is large enough, S ∼ 0.6 (as reported by [39] for submerged, laminar swirling jets), the flow structure changes and a vortex breakdown phenomenon occurs [40], producing a hollow-cone spray jet. Thus, we integrate the main equations (14)-(18), together with the corrections (40) and (41) for swirl motion at the injection, in a 2D axisymmetric, non-dimensional domain of radius rmax = 50 and axial coordinates z ∈ [0, zmax ], with zmax = 100; the injection of the spray takes place at z = 0 with a radial coordinate r < 1. The numerical simulations run during a time interval large enough, t ∈ [0, tend ] with tend = 500, that spray jet self-ignition take place and a stationary solution of the fully developed jet can be obtained; a constant, dimensionless time step size Δt = 0.02 is used throughout. The number of droplets at the injection is N˙ dI = 2000 (40 new droplets each time step). For the spatial discretization, we use a dynamic, anisotropically refined mesh with maximum and minimum element sizes hmax = 2.5 and hmin = 0.02, respectively. As fuel, we consider methanol CH3 OH at the atmospheric pressure and temperature Ts = 300K; all fuel properties are taken from [15], and are collected in Table 4.

Methanol

cl /cp 2.5

MF /MN2 1.14

RF /cp 0.25

LF 1.2

TB /Ts 1.123

Ta /Ts 20

lv = Lv /(cp Ts ) 3.67

q = Q/(cp Ts ) 66.43

s 1.5

Table 4: Methanol properties for Test III.

A summary of the appropriate boundary conditions for the gas phase used for this configuration is as follows: vz + 1 = vr = Γ − 2r2 = YF = YO2 = T − 1 = 0, ∂YO2 ∂T ∂YF = = = 0, vz = vr = Γ = ∂z ∂z ∂z  Pr  p·n+ μ ∇v + ∇vT · n = ∇Γ · n =∇YF · n = ∇YO2 · n = ∇T · n = 0, Pe  Pr  p·n+ μ ∇v + ∇vT · n = ∇Γ · n =∇YF · n = ∇YO2 · n = ∇T · n = 0, Pe

z=0

and r ≤ 1,

z=0

and r > 1,

z = zmax , r = rmax .

For the droplets at the injection zk = 0 and 0 < rk < 1, the following initial conditions apply: vkz − 1 = vkr = Γk − 2rk2 = Tk − 1 = ak − 1 = 0. 28

Finally, to completely define the simulations, we must provide the Peclet number Pe = 500, Damk¨ ohler number Da = 104 , and mass loading-ratio α = 0.05. In the next paragraphs, we study the influence of the Stokes number St and the swirl ratio S, showing results for the stationary solution, while the temporal evolution from t = 0 to t = tend may be found in the supplementary material to this paper. 4.3.1. Spray jets without swirl First, we perform spray numerical simulations without swirl motion S = 0, considering both non-inertial (St = 0.2) as well as inertial (St = 2.0) droplets. The impinging jet starts developing a head vortex which travels downstream and eventually leaves the domain, the jet reaching a stationary or quasi-steady solution. In Figure 14 we can observe the stationary solution at the final instant of time tend , with the fully developed jet. Although in both simulations the number of droplets at the injection is the same, N˙ dI = 2000, the final stationary solution has a larger number of droplets when inertial droplets are considered: Nd = 74000 for St = 2.0, and Nd = 47500 for St = 0.2. This is so because inertial droplets have a longer vaporization time (since vaporization and accommodation times are of the same order tv ∼ ta ), which means longer penetration distances. However, the length of the jet, marked by the stoichiometric (solid black) line, is slightly shorter for inertial droplets St = 2.0; as to jet length, it depends linearly on the Peclet number.

Figure 14: Fuel injected without swirl S = 0 and Stokes numbers St = 0.2 (above) and St = 2.0 (below). Temperature is shown in colormap, with pink dots denoting 10% of the computational droplets in the domain. The solid black line represents the stoichiometric line where the diffusion spray flame is situated. The anisotropically-adapted mesh is also shown in the lower half of each panel.

In the stationary state, the spray jet flame has the typical structure of a triple flame located at the self-ignition coordinate zign , with zign ∼ 7 for St = 0.2, and zign ∼ 15 for St = 2.0. At that point, a premixed flame (lean in the side of air and rich in the side of fuel vapor) appears, producing downstream a diffusion flame placed along the stoichiometric jet line. A detail of this region is shown in Figure 15: for each panel, the top half shows isolines of the reaction rate along with the droplet distribution, while the bottom half depicts the adapted mesh in the mixing layer and the triple flame region. We also notice that the inertial spray (St = 2.0) presents more instabilities in the mixing layer and leads to a larger number of elements N E = 6500 than the non-inertial spray (St = 0.2), for which N E = 3500. Such difference may be ascribed to combustion, since the high temperature zone starts later for inertial droplets (St = 2.0); thus, the instabilities at the shear layer (formed by the jet stream impinging at unity velocity in the quiescent atmosphere) have more space to develop, since at high temperature, molecular diffusion increase (T σ ) and instabilities are reduced. 4.3.2. Spray jets with swirl Next, we provide both the spray and the carried inert stream of nitrogen with an azimuthal velocity whose intensity is given in terms of the swirl rate S. For S = 0.5, the impinging jet does not develop a 29

3

r

3

St = 0.2

2

2

1

1

0

0

1

1

2 3

4

6

8

10

St = 2.0

2

z 2

r

3 12

12

z 14

16

18

20

22

Figure 15: Detail of the self-ignition zone of the fuel injector test without swirl (S = 0) and Stokes numbers St = 0.2 (left) and St = 2.0 (right). The top half of each panel show the isolines of reaction rate in color with the stoichiometric line in solid black; pink points represent 100% of the droplets. The bottom-half panels show the anisotropically-adapted mesh with the elements clustered around the mixing layer and the triple flame location.

head vortex as it happened before (S = 0) at the first stages of the simulation; instead, the jet presents a small, hollow-cone structure that owing to instabilities, evolves once again into a slender jet structure. The stationary solution for S = 0.5, considering both non-inertial and inertial droplets is shown in Figure 16.

Figure 16: Fuel injected with swirl S = 0.5, and Stokes numbers St = 0.2 (above) and St = 2.0 (below). Temperature is shown in colormap, with pink dots denoting 10% of the computational droplets in the domain. The solid black line represents the stoichiometric line where the diffusion spray flame is located. The anisotropically-adapted mesh is also shown in the lower half of each panel.

If we compare these results for S = 0.5 (Figure 16) with the previous case where no swirl motion was considered S = 0 (Figure 14), we observe how the ignition distance is now located closer to the ignition plane zign ∼ 0, and the developed jet is a little bit shorter and wider than for S = 0. This can be explained as swirl motion speeds up the vaporization process, so the number of droplets in the domain is sharply reduced, while reaching higher temperatures for S = 0.5 than for S = 0. We also observe, for this test with motion, a similar degree of oscillations when non-inertial St = 0.2 and inertial droplets St = 2.0 are considered; the oscillations in the shear layer are induced by the discontinuities in both the axial and azimuthal velocities between the jet and air coflow streams. As a consequence, the computational requirements in both cases, non-inertial and inertial droplets, are quite similar: for St = 0.2 we end up with N E = 8500 mesh elements and Nd = 14100 droplets, whereas for St = 2.0 we have N E = 7800 elements and Nd = 11300 droplets. Finally, and to highlight the capabilities of the numerical method, we consider a situation with a more intense swirl motion, S = 0.75; under such circumstances, the swirl motion induces an attached vortex downstream of the jet, which reduces the axial velocity of the fluid, appearing a steady, hollow-cone configuration for the spray jet flame. This is the so-called “vortex breakdown”, a fluid mechanics phenomenon 30

well known in the literature [40]. Thus, we show in Figure 17 the stationary solution at the final instant of time tend . The colormap for the temperature shows much larger values for St = 2.0 than for non-inertial droplets St = 0.2. In addition, when the hollow-cone configuration appears droplet penetration in the domain is severely dampened (droplets are heated and vaporized very close to the swirl injector), and numerical simulations may be carried out scattering as low as Nd ∼ 7000 droplets in the domain, in either case St = 0.2 and St = 2.0. Another important feature of this configuration is that the size and the shape of the hollow-cone is weakly dependent on the Peclet number, which also diminishes the requirements for spatial resolution, with some N E ∼ 4500 mesh elements being used. Each inset of Figure 17 represents one region of interest in the flow developed outside of the swirling fuel injector: the counterflow configuration previously analyzed in test I, and the self-ignition zone with triple flame structures studied in detail in test II. Thus, the lower insets show a detail of the central region of the hollow-cone, with the counterflow structure of the solution developed due to the stagnation point generated by the vortex. In solid pink line we plot the trajectories of the droplets, while dashed white lines represent the streamlines of the gas phase. We can see that for the non-inertial case St = 0.2, the droplets follow the streamlines and do not cross the stagnation surface located at the diffusion flame (marked by the stoichiometric, solid black, line) inside the mixing layer. In contrast, inertial droplets St = 2.0 do cross the stagnation surface, being fully vaporized just after turning around, moving the stoichiometric line (and the diffusion flame) outside the mixing layer. A similar behavior was observed in test I. The upper insets of Figure 17 plot a detail of the self-ignition zone, following our previous analysis for test II; we also show the spatially-adapted mesh, the isolines of the chemical reaction rate and the droplet distribution. We can observe an increase in the ignition distance caused by the axial velocities of the spray and coflow air streams being of similar magnitude, an interesting effect that could be a relevant design parameter of industrial burners.

Figure 17: Fuel injected with swirl S = 0.75 and Stokes numbers St = 0.2 (left) and St = 2.0 (right). Temperature is shown in colormap, with pink dots denoting 10% of the computational droplets in the domain. The solid black line represents the stoichiometric line where the diffusion spray flame is located. Insets show details of the counterflow configuration analyzed in test I, and the self-ignition zone with triple flame structures studied in test II.

5. Conclusions and future work In this work, we have presented a novel, Eulerian-Lagrangian method to solve time-dependent problems of diluted spray combustion at low Mach and moderately high Peclet numbers (laminar flows). The mathematical model of the problem considers a special two-phase system, in which a large number of tiny droplets with negligible neighboring interaction (liquid phase) are surrounded by a gas flow (gas phase). The conservation

31

equations for the gas phase are solved via a Lagrange-Galerkin formulation involving a Lagrangian stage for the convective terms (using a “modified Lagrange-Galerkin” scheme) along with an Eulerian stage to deal with the diffusion-reaction part of the equations; this formulation allows for a decoupling of the Navier-Stokes equations from the thermochemical problem (energy and species mass fraction conservation equations). A set of ordinary differential equations further describe the evolution of the properties of each droplet in a Lagrangian formulation; together with the thermochemical conservation equations for the gas phase, these droplet equations are solved by an explicit Runge-Kutta-Chebyshev method suitable for highly-stiff, non-linear systems of equations featuring a strong coupling between the variables. All the previous numerical techniques are implemented within an anisotropic, locally adaptive mesh procedure which reduces the computational cost of the simulations to such an extent that they may be run in commodity personal computers. To gain insight into the complex dynamics of spray flames and for the purposes of code validation we have investigated a series of canonical problems: in test I, we analyze the structure of diffusion flames in counterflow axisymmetric configurations, studying the different behavior of non-inertial (St < 1) and inertial (St > 1) droplets; in test II, we carefully examine the self-ignition region of the spray produced by the high temperature of the surrounding gas in planar mixing layers; in both test cases, our results compare very satisfactorily with those of the literature. Both test problems serve as an introduction to more complex, real-world configurations such as the outer flow of swirling fuel spray injectors presented in our last example test III, whose flow field encompasses the two phenomena previously analyzed. Nevertheless, our current numerical scheme will benefit from future work; thus, we plan to devise a more realistic physical mode, accounting for droplet interaction and including a more complex chemical kinetic mechanism. Such additional complexity should present no challenge from a numerical point of view, as its calculation would be handled straightforwardly by the explicit Runge-Kutta-Chebyshev method. In addition, a high-performance computing implementation of the current algorithms is in the works: the parallelization of both the Lagrangian droplet stage and the Ex-RKC method should be effortless; and the solution of the linear systems of equations, including the saddle-point, Stokes problem, can be further improved using optimized libraries such as PETSc which works seamlessly with MPI. We also aim at replacing the Dirac delta distribution by other, smoother functions to reduce random oscillations in the solution, trying to arrive at a minimal statistical dispersion even when using a small number of droplets at the injection N˙ dI . Finally, we expect interface capturing methods in combination with anisotropic mesh adaptation to be extremely helpful in accurately defining the moving interface when studying the zone closer to the swirl atomizer α(x)  1 where the liquid fuel film is separated from the gas phase by a free surface. Acknowledgments This research has been funded by project MTM2015-67030-P of the Spanish “Ministerio de Econom´ıa y Competitividad”. The authors would like to thank Professors A. S´ anchez and A. Li˜ n´an their priceless dedication and fruitful discussions, which have tremendously helped in our understanding of spray combustion phenomena. Supplementary data Supplementary data can be found in the online version of this article. References [1] A.H. Lefebvre. Atomization and Sprays, Hemisphere Publishing Corporation (1989). [2] A.L. S´ anchez, J. Urzay, A. Li˜ n´ an(2015). The role of separation of scales in the description of spray combustion. Proc. Combust. Inst. 35 (2015) 1549-1577. [3] O. Desjardins, V. Moureau, H. Pitsch. An accurate conservative level set/ghost fluid method for simulating turbulent atomization. J. Comp. Phys. 227. (2008) 8395-8416.

32

[4] V. Le Chenadec, H. Pitsch. A monotonicity preserving conservative sharp interface flow solver for high density ratio two-phase flows. J. Comp. Phys. 249. (2013) 185-203. [5] J. Shinjo, A. Umemura. Detailed simulation of primary atomization mechanisms in Diesel jet sprays (isolated identification of liquid jet tip effects). Proc Comb. Inst. 33 (2011) 2089-2097. [6] J.L. Prieto, J. Carpio. A-SLEIPNNIR: a multiscale, anisotropic adaptive, particle level set framework for moving interfaces. Application to advection-dominated flows. J. Comp. Phys. 377. (2019) 89-116. [7] W.A. Sirignano. Fluid Dynamics and Transport of Droplets and Sprays. Cambridge University Press. 2nd Ed.(2010). [8] A. Allievi, R. Bermejo. A Generalized Particle Search–Locate Algorithm for Arbitrary Grids. J. Comput. Phys. 132 (1997) 157-166. udez, J.L. Ferr´ın, A. Li˜ na´n, L. Saavedra, Numerical simulation of group combustion of pulverized [9] A. Berm´ coal, Combust. Flame 158(9) (2011) 1852–1865. [10] X.Y. Zhao, D.C. Haworth, Transported PDF modeling of pulverized coal jet flames, Combust. Flame 161 (7) (2014), 1866–1882. [11] N. Patel, M. Kirtas, V. Sankaran, S. Menon, Simulation of spray combustion in a lean-direct injection combustor. Proc. Combust. Inst. 31, (2017) 2327–2334. [12] G. Hannebique, P. Sierra, E. Riber, B. Cuenot, Large eddy simulation of reactive two- phase flow in an aeronautical multipoint burner, Flow Turbul. Combust 90 (2) (2013), 449–469. [13] K. Luo, H. Pitsch, M.G. Pai, O. Desjardins. Direct numerical simulations and analysis of three-dimensional n-heptane spray flames in a model swirl combustor. Proc. Combust. Inst. 33 (2011) 2143-2152. na´n, D. Mart´ınez-Ruiz, A. L. S´ anchez, J. Urzay. Regimes of spray vaporization and combustion in [14] A. Li˜ counterflow configurations. Combust. Sci. Tech., 187 103–131 (2015). n´ an, F. A. Williams. Dynamics of thermal ignition of [15] D. Mart´ınez-Ruiz, J. Urzay, A. L. S´anchez, A. Li˜ fuel sprays in mixing layers. J. Fluid Mec., 734 387–423 (2013). [16] F. Doisneau, F. Laurent, A. Murrone, J. Dupays, M. Massot. Eulerian Multi-Fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion. J. Comput. Phys. 234 (2013) 230–262. [17] A. Sibra, J. Dupays, A. Murrone, F. Laurent, M. Massot. Simulation of reactive polydisperse sprays strongly coupled to unsteady flows in solid rocket motors: Efficient strategy using Eulerian Multi-Fluid methods. J. Comput. Phys. 339 (2017) 210–246. [18] F. Doisneau, M. Arienti, J.C. Oefelein. A semi-Lagrangian transport method for kinetic problems with application to dense-to-dilute polydisperse reacting spray flows. J. Comput. Phys. 329 (2017) 48–72. [19] J. Carpio, J.L. Prieto, M. Vera. Local anisotropic adaptive algorithm for the solution of low-Mach transient combustion problems. Journal of Computational Physics, Vol. 306 (2016), pp. 19 - 42. [20] J. Carpio, J.L. Prieto. An anisotropic, fully adaptive algorithm for the solution of convection dominated equations with semi-Lagrangian schemes. Computer Methods in Applied Mechanics and Engineering, Vol. 273 (2014), pp. 77 - 99. [21] R. Bermejo, J. Carpio. An adaptive finite element semi-Lagrangian Runge-Kutta-Chebychev method for convection dominated reaction-diffusion problems. Applied Numerical Mathematics, 58 (2008) 16 - 39. [22] R. Bermejo, J. Carpio. A semi-Lagrangian-Galerkin projection scheme for convection equations. IMA Journal of Numerical Analysis, Vol. 30(3) (2010), pp. 799 - 831. 33

[23] O. Pironneau. On the transport-diffusion algorithm and its applications to the Navier-Stokes equations, Numer. Math. 38(3) (1982) 309–332. [24] A. Allievi, R. Bermejo, Finite element modified method of characteristics for Navier-Stokes equations, Int. J. Numer. Meth. Fl. 32 (2000) 439–464. [25] B.P. Sommeijer, L.F. Shampine, J.G. Verwer, RKC: An explicit solver for parabolic PDEs, J. Comput. Appl. Math. 88(2) (1998) 315–326. [26] R. Bermejo, P. Gal´an del Sastre, L. Saavedra, A second order in time modified Lagrange-Galerkin finite element method for the incompressible Navier-Stokes equations, SIAM J. Numer. Anal. 50(6) (2012) 3084–3109. [27] E.J. Dean, R. Glowinsky, On some finite elements methods for the numerical simulation of incompressible viscous flow, in: M.D. Gunzburger, R.A. Nicolaides (Eds.), Incompressible Computational Fluid Dynamics: Trends and Advances, Cambridge University Press, Cambridge, 1993, pp. 17–66. [28] S. Balay, S. Abhyankar, M.F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout, W.D. Gropp, D. Kaushik, M.G. Knepley, L.C. McInnes, K. Rupp, B.F. Smith, S. Zampini, H. Zhang, and H. Zhang, PETSc Users Manual, Argonne National Laboratory, 2016, ANL-95/11 - Revision 3.7, http://www.mcs.anl.gov/petsc. [29] J.L. Prieto. SLEIPNNIR: A multiscale, particle level set method for Newtonian and non-Newtonian interface flows. Computer Methods in Applied Mechanics and Engineering, Vol. 307 (2016), pp. 164 192. [30] P. Frey, F. Alauzet. Anisotropic mesh adaptation for CFD computations. Comput. Methods Appl. Mech. Eng. 194(48) (2005) 5068-5082. [31] M.J. Castro-D´ıaz, F. Hecht, B. Mohammadi and O. Pironneau. Anisotropic unstructured mesh adaptation for flow simulations. Int. J. Numer. Methods Fluids 25(4) (1997) 475-491. [32] T. Coupez, Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing, J. Comput. Phys. 230(7) (2011) 2391–2405. [33] W. Huang, Metric tensors for anisotropic mesh generation, J. Comput. Phys. 204 (2005) 633–665. [34] J. Carpio, J.L. Prieto, R. Bermejo, Anisotropic “goal-oriented” mesh adaptivity for elliptic problems, SIAM J. Sci. Comput. 35(2) (2013) A861–A865. [35] F. Hecht, BAMG: Bidimensional Anisotropic Mesh Generator, INRIA-Rocquencourt, France, 2006, https://www.ljll.math.upmc.fr/hecht/ftp/bamg. [36] P. Frey. Mmg3d: 3d Delaunay-based anisotropic mesh adaptation and mesh moving, 2013, https://www.mmgtools.org. [37] G. Rokos, G. Gorman. PRAgMaTIc-parallel anisotropic adaptive mesh toolkit, in: Facing the MulticoreChallenge III, vol. 7686 of Lecture Notes in Computer Science, Springer Berlin Heidelberg, 143-144, 2013. https://github.com/meshadaptation/pragmatic [38] J. Carpio, A. Li˜ n´an, A.L. S´ anchez, F.A. Williams. Aerodynamics of axisymmetric counterflowing jet atmospheres. Combut. Flame 177 (2017) 137 - 143. anchez, and A. Li˜ n´an. The quasi-cylindrical description of submerged laminar [39] A. Revuelta, A.L. S´ swirling jets, Phys. Fluids 16(3) (2004) 848–851. [40] S. Leibovich. The structure of vortex breakdown, Annu. Rev. Fluid Mech. 10 (1978) 221–246

34