|∇fKSZ|=
rNf,KSZ
2D . (3.31)
The characteristic cross area of the wake structure is given by 2Dx1/Ud[22]. Follow-ing Eqs. (3.29) and (3.30), the PDF in the Kolmogorov scale zone can be rewritten as
Pf,KSZ∗ = 2cDx1/Ud qNf,KSZ
2D
. (3.32)
Note that fd is located at ∆x (mesh size) from the droplet in CP-DNS because fd should be interpreted as the value of the mixture fraction of the DNS cells where the droplets are located. These DNS cells can be understood as some virtual ”droplets”
or rather fuel sources with the mass fraction of evaporating species of fd, and the evaporated mass of the Lagrangian droplet is immediately distributed over the entire local Eulerian cell and then diffused to the inter-droplet space. In addition, assume that the average mixture fraction fg is approximately located at the inter-droplet distancerc. Then, the relationship betweenx1 andf can be approximately obtained by the interpolation
1 x1 = 1
∆x− fd−f fd−fg
1
∆x − 1 rc
. (3.33)
3.3 Summary
The scaling laws for the near droplet zone (NDZ, Sec. 2.4) and the Kolmogorov scale zone (KSZ, Sec. 2.4) are summarized in Tab. 3.1. There, the reader can find the location of the respective zones which will be quantified by the fully resolved DNS in Chapter 5, the determining parameters for mixing, the expressions for the characteristic mixing quantities and their corresponding equation numbers as listed in Chapters 2 and 3, and the chapters in this dissertation and the references where the scaling laws are introduced, derived and assessed.
Table 3.1: The scaling laws for scalar mixing in inter-droplet space
NDZ Expression Equation
#
Chapter Reference
Location x1<2Udτη Ch. 5 [145]
Parameters D,c,Jm,Ud,ρ f(x1, r) fN DZ=f2+J4πρDxm(fd−f2)
1φnexp
−4DxUd
1φnr2 , fmax,N DZ =f2+J4πρDxm(fd−f2)
1φn
Eqs. (2.37), (2.38)
Chs. 2, 5
[67,157]
Nf Nf,N DZ =(R2Ud(f−f2)2
c−∆x)φn
ln J
m(fd−f2) 4πρDφn(f−f2)
ln∆xRc +ln22∆x−ln22Rc
Eq. (2.41) Chs. 2, Apps.
A,B
[145]
Pf Pf,N DZ= cJρNm(fd−f2)
f,N DZ Eq. (2.43) Chs. 2,
Apps.
A,B
[67,145]
KSZ Expression Equation
#
Chapter Reference
Location x1>2Udτη Ch. 5 [145]
Parameters ε,c,Jm,Ud,ρ
f(x1, r) fKSZ=f1+ (φc−φ1)Jm(fρUd−f2)
dε
x1 Ud
−3
exp
−(φc−φεx1)πUd3
13 r2 ,
fmax,KSZ =f1+ (φc−φ1)Jm(fρUd−f2)
dε
x1
Ud
−3
Eqs. (2.39), (2.40)
Chs. 2, 3,5
[143, 145]
Nf Nf,KSZ =A1
1(f−f1)7/3J
m(fd−f2) ρUdε
−1/3
Eq. (2.42) Chs. 2, 3,5,6,7
[143, 144, 145]
Pf Pf,KSZ= cJρNm(fd−f2)
f,KSZ , Pf,G= √ 1
2πff02
exp
−(f−f˜)2
2ff02
(f−f1)a, β-PDF,
Pf,KSZ∗ = q2cDxNf,KSZ1/Ud
2D
Eqs. (2.44), (2.45), (3.32)
Chs. 2, 3,5,6,7
[67,143, 144, 145]
3.3. SUMMARY 30
Computational Approach
Chapter 4introduces the computational frameworks of fully resolved DNS and CP-DNS. Fully resolved DNS will form the basis for the detailed analysis of all essential local interactions within the gas phase including the boundary layer at the droplet surface and the smallest turbulence scales. The standard governing equations in both gas and liquid phases and the boundary conditions at the droplet surface and droplet center are presented in the first section. The random character of the droplet position is neglected in fully resolved DNS and CP-DNS can be used to overcome the limitation. The standard governing equations in fully resolved DNS are extended by the Lagrangian formulation where the droplets are treated as point sources of mass, momentum and energy in CP-DNS. The coupling between the two phases are described in the second section. Another two sections present the inflow turbulence generation method and the discretization schemes for solving the governing equations.
4.1 Fully Resolved DNS
Fully resolved DNS are used to simulate the combustion of droplet arrays in turbu-lent convective flows. The standard conservation equations for mass, momentum, energy, species and an additional transport equation for the conserved scalar mixture fraction are solved simultaneously for the gas and liquid phases. They are coupled by boundary conditions at the droplet surface which include the mass (flux) conser-vation, heat conduction and species diffusion to ensure the correct coupling. Heat transfer by radiation is neglected. The gas phase mixture is considered to be ideal with unity Lewis numbers and the Schmidt and Prandtl number are specified as Sc = Pr = 0.7. The thermodynamic properties of the liquid phase are assumed
con-31
4.1. FULLY RESOLVED DNS 32 stant. The effect of gravity and other body forces are small and can be neglected.
4.1.1 Governing Equations for Gas Phase
The mass continuity equation is given by [101, 126]
∂ρ
∂t + ∂(ρuj)
∂xj = 0, (4.1)
where ρ is the density, uj denotes the velocity.
The momentum transport equation is expressed as [101,126]
∂(ρui)
∂t +∂(ρujui)
∂xj =−∂p
∂xi +∂τij
∂xj, (4.2)
where p is the pressure. For Newtonian fluid, the viscous stress, τij, is given by τij =µ
∂ui
∂xj
+ ∂uj
∂xi
−2 3µ∂uk
∂xk
δij, (4.3)
where µ represents the dynamic viscosity, δij represents the Kronecker delta [101]
and its definition is given by Eq. (2.8).
With the assumption of unity Lewis number, the transport equation for each species is written as [22]
∂(ρYk)
∂t +∂(ρujYk)
∂xj = ∂
∂xj µ
Sc
∂Yk
∂xj
+ωk, (4.4)
where Yk is the mass fraction of the kth species and its chemical reaction rate is denoted by ωk. The chemical reaction source term, ωk, is the sum of reaction rates for all reactions involving the kth species, viz.
ωk =
Nr
X
i=1
(vki00 −vki0 )MkRi, (4.5) where Nr is the number of elementary reactions in the chemistry mechanism, vki00 and v0ki denote the reverse and forward stoichiometric coefficients of the kth species in the ith reaction, and Mk is the molar mass of thekth species. The reaction rate, Ri, is calculated by
Ri =kf,i
N
Y
k=1
ckv0ki −kr,i
N
Y
k=1
ckvki00, (4.6)
where N is the number of species and ck is the molar concentration of the kth species. The forward reaction rate constant of the ith reaction, kf,i, is generally modelled by the Arrhenius law [22, 127] which is expressed as
kf,i=AiTbiexp
− Ei RuT
, (4.7)
where Ai, bi and Ei are the pre-exponential factor, the temperature exponent and the activation energy, respectively, and they are specified for the ith reaction in the chemistry mechanism. Ru is the ideal gas constant. The reverse reaction rate constant of the ith reaction, kr,i, is given by
kr,i = kf,i Kc,i
. (4.8)
Here, Kc,i is the equilibrium constant of the ith reaction.
The energy transport equation is given by [22]
∂(ρhs)
∂t +∂(ρujhs)
∂xj = ∂
∂xj µ
Pr
∂hs
∂xj
+ωQ, (4.9)
where hs is the sensible enthalpy. The chemical reaction heat release rate, ωQ, is calculated by
ωQ=−
N
X
k=1
hkωk, (4.10)
where hk is the absolute enthalpy of the kth species which equals the sum of the sensible enthalpy and the standard enthalpy of formation [127].
An additional conserved transport equation for mixture fraction is solved as
∂(ρf)
∂t + ∂(ρujf)
∂xj = ∂
∂xj µ
Sc
∂f
∂xj
. (4.11)
4.1.2 Governing Equations for Liquid Phase
The liquid phase is represented by a single component fuel. Similar to gas phase, the governing equations for the mass continuity, momentum and energy transport are given by
∂ρl
∂t +∂(ρluj)
∂xj = 0, (4.12)
∂(ρlui)
∂t + ∂(ρlujui)
∂xj =−∂p
∂xi +∂τij
∂xj, (4.13)
4.1. FULLY RESOLVED DNS 34
∂(ρlhs)
∂t +∂(ρlujhs)
∂xj = ∂
∂xj µl
Prl
∂hs
∂xj
, (4.14)
where ρl, µl and Prl are the density, dynamic viscosity and Prandtl number of the liquid.
4.1.3 Boundary Conditions
The mass flux of evaporation, ˙m00, can be computed by a Robin type boundary condition [51] and is conserved on both sides of the droplet surface,
˙
m00YF l=YF g,sm˙00− µ Sc
∂YF g,s
∂n , (4.15)
˙
m00l,s = ˙m00, (4.16)
where YF l is a mass fraction of a specific fuel species in the liquid and set to 1 for a single component droplet, YF g,s represents the mass fraction of fuel vapour at the droplet surface on the gas side, ˙m00l,s is the mass flux on the liquid side and n is the normal distance from the droplet surface. The temperature at the droplet surface can be obtained by equating the energy fluxes between the liquid and gas phases,
µCp Sc
∂Tg,s
∂n =λl∂Tl,s
∂n + ˙m00hf g, (4.17) with
Tg,s=Tl,s, (4.18)
where Tg,s and Tl,s denote surface temperature on the gas side and surface tem-perature on the liquid phase side, respectively. The quantity hf g is the enthalpy of evaporation andλldenotes the heat conductivity of the liquid phase. The mass frac-tion of the evaporating species at the droplet surface, YF g,s, can be obtained using the Clausius-Clapeyron equation [127]. The mass fractions of the non-evaporating species (Yig,s) and the mixture fraction (fg,s) on the gas side are also given by a Robin type boundary condition,
˙
m00Yig,s− µ Sc
∂Yig,s
∂n = 0, (4.19)
˙
m00fd =fg,sm˙00− µ Sc
∂f
∂n, (4.20)
where the mixture fraction on the liquid sidefdis 1 as we consider a single component fuel. The boundary conditions at the droplet center can be expressed as
∂Tl
∂n = 0, (4.21)
∂ul
∂n = 0, (4.22)
whereTl andul represent the temperature and velocity inside the droplet. A quasi-steady state assumption is applied for the droplet size. The characteristic evapora-tion time scale (τv) is estimated to be up to one order of magnitude larger than the Kolmogorov time scale (τη) for the cases investigated here. Small turbulent scales affect the mixing fields faster than the evaporation process proceeds which facilitates the scaling law analysis. Temporal invariance of the droplet size is enforced and the mass evaporating at the droplet surface is replenished from the droplet center (cf.
Eqs. (4.21) and (4.22)) to ensure a realistic temperature profile within the droplet and thus energy balance at the droplet surface. This requires the additional solution of the mass, momentum and energy equation in the liquid and matching boundary conditions at the droplet surface (Eq. (4.16)). A truly transient evaporation pro-cess can then be investigated based on a series of quasi-steady states. Similarly, an analysis of the drag forces acting on the droplets shows that the relevant time scale governing relative motion between the droplets is large compared to the evaporation time and can therefore be neglected. This was also demonstrated in earlier compu-tations [158, 159] where both phases were fully coupled and a receding liquid/gas interface was considered in the computation. Therefore, a constant mean velocity of the inflow can be set during the simulations.