Transport and Mixing in Liquid Phase Using Large Eddy Simulation: A Review

Abstract

Many mixing processes in engineering applications are turbulent. At high‐Schmidt regime, the scalar scales are much lower than those of the velocity field, making difficult instantaneous measurements and direct numerical simulation for studying systems of practical interest. The use of large eddy simulation (LES) for analyzing transport and mixing of passive and reactive scalars at high‐Schmidt (Sc) regime is addressed in this article. We present two different approaches for studying scalar transport and mixing in LES: the conventional approach is based on the modeling of the unclosed subgrid‐scale scalar flux term in the filtered scalar equation by models commonly used for high‐Sc flows. The second approach presented in this review for dealing with high‐Sc flows is based on the use of a filtered mass density function (FDF) of the scalar field. Conclusions are presented about the relative merits of the two approaches.

1. Introduction

At molecular scale, the mixing of two or more fluids having different composition is driven by diffusion, a rather slow process. Turbulence increases the rate of mixing by the action of large‐scale motions, making higher the contact surface area between adjacent unmixed fluid “elements.” Turbulent mixing plays an important role in many engineering, biological, and environmental applications. In the case of reactive systems, chemical reactions can only take place after reactants are mixed at molecular level; in many of such systems, the availability of mixed reactants at the molecular scale limits the reaction rate; consequently, the mixing process can be a controlling rate.

In the chemical and pharmaceutical industries, an inadequate mixing can raise the production costs due to the reduction in selectivity, low yield, undesired product accumulation, and scale‐up and process development problems. Some costs related to mixing are reported by Paul et al. [1]: in 1989, the cost of poor mixing in the U.S. chemical industry was estimated at $1 to $10 billion; yield losses of 5% due to poor mixing are typical. In the pharmaceutical industry, the costs due to lower yield and due to problems in scale‐up and process development are on the order of $100 and $500 million, respectively. In the 1980s, the pulp and paper industry reported savings averaging 10–15% by introducing medium consistency mixer technology in their processes.

Numerous studies have dealt with transport and mixing of passive scalars in systems where the Schmidt number (Sc) is close to one, as it is the case for gas flows in which the smallest scalar scale (the Obukhov‐Corrsin scale) is of the same order of magnitude or higher than the smallest flow scale, the Kolmogorov scale (see Refs. [2–5]). However, in many industrial and biological applications that take place in the liquid phase, the ratio of the fluid momentum diffusivity to molecular diffusivity of the scalar is much greater than 1; so, the scalar field holds finer structures than the velocity field (see Refs. [6–9]).The smallest scalar length‐scale when the Schmidt number is greater than 1 is named the Batchelor scale (η_B) and is much smaller than the Obukhov‐Corrsin scale, η_oc, and the smallest velocity length‐scale,—the Kolmogorov scale (η). The small dimension of the Batchelor scale makes accurate measurement of the instantaneous concentration in turbulent liquid flows scarce and only possible with rather sophisticated systems [10].

To overcome many of the difficulties found in mixture‐based processes, numerical simulation appears as an important tool for design and process improvement, as it can predict important information about the flow and scalar fields, which is hard to measure. When numerical simulations of industrial process are carried out, key requirements are the computational efficiency, algorithm robustness, and accurate representation of the process.

Models that correctly represent the mixing processes in turbulent flows are important for the design of various engineering and environmental applications where turbulence plays an important role. While the current understanding of turbulence allows a reasonable description of the flow field, a universal description of turbulent mixing processes remains a challenge [11, 12].

Direct numerical simulation (DNS) gives the most detailed information about the flow, as it resolves all the scalar length‐ and timescales by setting up the spatial resolution according to the smallest scalar/flow scale. No empirical closure or turbulence model assumptions are required. DNS has been useful for the study of transition and turbulent flows [13]. Because physics and chemistry are properly represented, DNS has also been used for the analysis of complex phenomena in combustion systems [14]. The use of DNS has been increasing as a complementary means of study of turbulent mixing process [15]. Sophisticated discretization schemes are used by DNS, imparting low flexibility when complex geometries are used [13]. The full‐scale resolution requirement makes very expensive the computational cost, and sometimes prohibitive, even for future computational capabilities. The computational cost rises as Re³. Approximately 99% of calculations are used to solve the dissipation scales [16]. For these reasons, DNS is not a viable choice for systems of practical importance at high Reynolds and Schmidt numbers.

A more realistic alternative to DNS is the utilization of spatial filtering or temporal ensemble averages, such as large eddy simulation (LES) and unsteady Reynolds averaged Navier‐Stokes equations (U‐RANS) [17]. When LES and U‐RANS formulations are used, additional terms appear in the transport equations. These terms need to be modeled. The information given by DNS has been useful to validate the unclosed terms of LES and RANS approaches [18].

RANS is based on the application of the Reynolds decomposition to any quantity $Q\left(x,t\right)=\overline{Q\left(x,t\right)}+Q\text{'}\left(x,t\right)$, where $\overline{Q\left(x,t\right)}$ is the mean of Q, and $Q\text{'}\left(x,t\right)$ is the deviation from the mean. The temporal average must be done for time intervals greater than any flow timescales. Averaged equations contain additional unknown terms: Reynolds stresses, scalar fluxes, and averaged source terms. These terms are unclosed, and turbulence models must be provided in order to close the mathematical system. In RANS, the largest scales are solved, and the turbulent spectrum information is provided by turbulence models. By using RANS, a good balance between results and computational cost is obtained. RANS is useful in engineering, especially for industrial and environmental sectors [19]. It is possible to perform parametric studies, considering its low computational cost. The main limitation of RANS models is that they only give limited information about turbulence, because all turbulent scales are modeled. Although RANS has achieved reasonable good precision for simple flows, in the reproduction of large‐scale organized and nonstationary turbulence structures, when strong streamline curvature or nongradient transport is present, RANS has achieved limited success [20]. Some inaccuracies arise from the turbulent viscosity hypothesis, the equation for ε and the nonuniversality of the model's constants which have to be tuned in order to improve the simulation results [16]. Although RANS provides a reasonable computational cost‐precision ratio, it does not always predict the concentration fluctuation distributions [21] in turbulent reacting liquid flows; in the case of stirred tank reactors (STR), the turbulent kinetic energy is unpredicted in the impeller region and discharge stream [22–25].

On the other hand, LES can be viewed as an intermediate approach between DNS and RANS, because in the former all turbulent structures are calculated, while in the latter they are modeled. The idea of LES is to solve the large, nonuniversal, anisotropic turbulent scales and to model the small turbulent scales, which contain less kinetic energy and are nearly universal and isotropic. These small scales are easier to model than the whole turbulent spectrum. The LES approach was proposed by Leonard [26]. A higher fidelity on the representation of the flow structure than RANS is expected, since the geometry‐dependent, large turbulent scales are calculated.

The cut‐off length must be proportional to the longitudinal integral length‐scale [16], L_EI. Typically, 92% of the kinetic energy is resolved for a one‐dimensional flow and 80% for a three‐dimensional flow [16]. The proportionality constant depends on the filter specification. The uncertainty about residual motions is lowered by reducing the filter width; therefore, constitutive LES equations are grid‐dependent, that is, LES is an incomplete model. References [27–29] review the LES literature.

Compared to DNS, LES saves the computational cost of solving the lower turbulent scales (99% of DNS calculations), and can be used for the simulation of systems of practical importance. However, it is more expensive than a k‐ε model for a given grid resolution.

In shear flows, LES prediction capabilities are diminished by the fact that in the viscous region, the energy‐containing structures are of the same order of magnitude as the viscous length‐scales (${\delta }_v=v\sqrt{\rho /{\tau }_w}$), and LES cannot solve them. To face this issue, it is possible to refine the computational mesh in near‐wall zones in order to consider the small coherent structures developed on those zones, but this refinement implies more computational effort. There are three additional possibilities: incorporating the fluid behavior in the wall region by wall models, simulating those zones with RANS (hybrid RANS‐LES), or incorporating wall‐adapted SGS closures.

Most research efforts in LES for reactive flows are concentrated on the SGS fluxes and Favre‐filtered source terms [30]. In diffusion flames, chemical reactions take place after mixing of reactive species is achieved at the smallest scales of turbulence (unresolved in LES); so, in LES the combustion process must be modeled.

According to Ref. [31], the uncertainty initially contained in the nonresolved scales is propagated to the resolved ones; hence, it is said that SGS modeling is not a well‐posed problem. LES results can be interpreted as a different realization of the flow. They could have the same statistical properties of the flow and may predict the same spatially organized structures but at a different location.

LES has been successfully applied to turbulent flows, in particular to complex geometries in liquid phase such as stirred tank reactors (see Refs. [10, 24, 32–34]). Because LES equations are unclosed, SGS models must be supplied in order to specify the SGS stresses, SGS scalar fluxes, and filtered reaction terms.

There are still unanswered questions about the behavior of high Schmidt turbulent flows. LES could be an alternative method for studying these flows, as long as the flow's physical behavior can be captured by the subgrid‐scale models [15]. The LES equations are presented in the following section.

2. Large eddy simulation (LES)

Spatial scale separation is done by applying a low‐pass filter to the governing balance equations. This set of filtered equations governs the dynamics of the large scales. Spatial fluctuations, lower than a defined filter cut‐off length ‐Δ‐, are smoothed or removed. The spatial filtering operation of a flow variable Q, being a function of time and space, is defined by the convolution integral:

${⟨Q⟩}_L\equiv \frac{1}{{⟨\rho ⟩}_l}{\displaystyle {\int }_D^{}\rho }G\left(x-x´\right)Q\left(x\right)dx\text{'}=\frac{{⟨\rho Q⟩}_l}{{⟨\rho ⟩}_l}$E1

where $Q={⟨Q⟩}_L+Q”$. Favre filtering has been applied to Eqs. (1) and (2). G is the filter kernel defined over the entire domain. The filter kernel depends on the filter width, Δ, and must satisfy a set of conditions (see Ref. [16]). Favre‐filtered mass, momentum, and scalars’ balance equations are as follows:

$\frac{\partial {⟨\rho ⟩}_l}{\partial t}+\frac{\partial }{\partial x_j}\left({⟨\rho ⟩}_l{⟨u_j⟩}_L\right)=0$E2

$\frac{\partial {⟨\rho ⟩}_l{⟨u_i⟩}_L}{\partial t}+\frac{\partial }{\partial x_j}\left({⟨\rho ⟩}_l{⟨u_i⟩}_L{⟨u_j⟩}_L\right)=-\frac{\partial }{\partial x_j}\left[{⟨\rho ⟩}_l\left({⟨{\tau }_{ij}⟩}_l+{\tau }_{ij}^{sgs}\right)\right]-\frac{\partial {⟨p⟩}_l}{\partial x_i}+{⟨\rho ⟩}_l{F}_i$E3

$\frac{\partial {⟨\rho ⟩}_l{⟨{\varphi }_{\alpha }⟩}_L}{\partial t}+\frac{\partial }{\partial x_j}\left({⟨\rho ⟩}_l\left({⟨{\varphi }_{\alpha }⟩}_L{⟨u_j⟩}_l+{⟨J_{\alpha ,j}^{}⟩}_L+J_{\alpha ,j}^{sgs}\right)\right)={⟨S_{\varphi }⟩}_L$E4

where the symbol ∂ denotes the partial differential operator, and the summation convention is used for repeated indices. Time (t) and spatial coordinates (x_j) are the independent variables; ρ is the fluid mass density; ${\varphi }_{\alpha }$ is the mass‐weighted value of the scalar field (e.g., internal energy, mass fraction); $u_j$ is the mass‐averaged velocity in the j‐direction; p denotes pressure; ${\tau }_{ij}$ is the shear stress tensor; and $J_{\alpha ,j}^{}$ is the diffusive flux of the scalar.

New terms appear in LES equations, as it occurs in RANS. These terms account for resolved and subgrid scales’ interactions. The SGS stresses $\left({\tau }_{ij}^{SGS}\right)$ and SGS scalar fluxes $\left(j_{\alpha ,j}^{SGS}\right)$ are defined by

${\tau }_{ij}^{SGS}={⟨u_i{u}_j⟩}_L-{⟨u_i⟩}_L{⟨u_j⟩}_L;j_{\alpha ,j}^{SGS}={⟨u_j{\varphi }_{\alpha }⟩}_L-{⟨u_j⟩}_L{⟨{\varphi }_{\alpha }⟩}_L$E5

In addition, a model is needed for filtered source terms. Subgrid‐scale (SGS) models must be supplied in order to specify the SGS stresses, SGS scalar fluxes, and filtered reaction terms.

The main task of SGS stress tensor model is to dissipate the energy transferred from turbulent large scales. The dynamics of the subgrid scales affect that of the resolved flow field through the subgrid‐scale stress tensor.

There are several SGS stress tensor models (for review, see Refs. [17, 35]). The simplest models are based on the gradient assumption. More complex models are of first and second order. As the model complexity increases, so does the required number of equations. Those models based on eddy viscosity show good balance between accuracy and numerical computation and present good prediction abilities in turbulent combustion and other highly interacting processes. It is also possible to close the LES approximation by PDF methods, such as the velocity‐filtered density function (VFDF) [36] and the velocity‐scalar filtered mass density function (VSFMDF) [37], but this approach increases considerably the number of equations to be solved. This procedure increases the computational cost 15–30 times longer than the Smagorinsky [38] and dynamic Smagorinsky [39, 40] models, respectively.

2.1. The Smagorinsky model

The Smagorinsky model [38] is one of the pioneer SGS models in the development of LES. Its simple formulation and good performance have made it very popular, and it is used in this work. The Smagorinsky model is an eddy‐viscosity based model, which assumes equilibrium between the turbulent kinetic energy dissipation and production rates to obtain a relation between the characteristic velocity and the resolved strain rate. A Boussinesq approximation is applied to the deviatoric part of the SGS stress.

${\tau }_{ij}^{SGS}-\frac{1}{3}{\tau }_{kk}^{SGS}{\delta }_{ij}^{}\approx -{\nu }_T\left[2{⟨S_{ij}⟩}_L-\frac{2}{3}⟨S_{kk}⟩{\delta }_{ij}^{}\right]$E6

where ${\nu }_T$ is the subgrid eddy viscosity of residual motions and $\delta$ is the Dirac delta function. The filtered strain rate is given by

${⟨S_{ij}⟩}_L=\frac{1}{2}\left(\frac{\partial {⟨u_{ij}⟩}_L}{\partial x_j}+\frac{\partial {⟨u_j⟩}_L}{\partial x_i}\right)$E7

The subgrid eddy viscosity is modeled in a similar way as the mixing length:

${\nu }_T={\left(C_s\Delta \right)}^2{\left|S\right|}_L,{\left|S\right|}_L\equiv {\left(2{⟨S_{ij}⟩}_L{⟨S_{ij}⟩}_L\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$E8

where C_s is the Smagorinsky coefficient. This Model parameter is not universal and depends on the flow configurations. While for isotropic turbulence, the Smagorinsky coefficient is about 0.17 [41], for other configurations this value may not be correct; for a channel flow, the Smagorinsky coefficient is about 0.1 [42]. In addition, the spatial variation of this coefficient makes it difficult to find a proper value. This model parameter can be dynamically calculated by using the dynamic procedure [39, 40], in which, by assuming scale invariance, a test filter greater than the filter width, Δ, can be used for calculating C_s from the resolved flow field.

2.2. The dynamic Smagorinsky model

The major drawback of the Smagorinsky model is that a single universal constant cannot correctly represent different turbulent flows. Germano et al. [39] proposed a dynamic procedure for computing the Smagorinsky model coefficient, based on the instantaneous information given by the filtered velocity field. The coefficient is locally recalculated during the simulations; so, it is no longer necessary to specify its value as an input parameter.

The dynamic procedure is based on the idea that the information given by the smallest resolved scales can be used to model the largest unresolved scales, as they have a similar behavior. The latter can be done by employing a test filter, with the test filter width larger than the filter width (Δ), usually taken as $\widehat{\Delta }=2\Delta$.

A model for the SGS stress is

${\tau }_{ij}^{SGS}-\frac{1}{3}{\tau }_{kk}^{SGS}{\delta }_{ij}^{}=-2C_{dyn}{\Delta }^2\left|{⟨S_{ij}⟩}_L\right|\left({⟨S_{ij}⟩}_L-\frac{1}{3}{⟨S_{kk}⟩}_L{\delta }_{ij}^{}\right)=-2C_{dyn}{\alpha }_{ij}^{}$E9

where C_dyn is a Model parameter, and the term should be modeled. If the test filter is applied to the filtered momentum equation, the filter stress has the form:

$T_{ij}^{}={\widehat{⟨u_i{u}_j⟩}}_L-{\widehat{⟨u_{ij}⟩}}_L{\widehat{⟨u_j⟩}}_L$E10

As it was the case for the SGS stress, a model for the filter stress is

$T_{ij}^{}-\frac{1}{3}{T}_{kk}^{}{\delta }_{ij}^{}=-2C_{dyn}{\widehat{\Delta }}^2\left|{\widehat{⟨S_{ij}⟩}}_L\right|\left({\widehat{⟨S_{ij}⟩}}_L-\frac{1}{3}{\widehat{⟨S_{kk}⟩}}_L{\delta }_{ij}^{}\right)=-2C_{dyn}{\beta }_{ij}^{}$E11

The coefficient has been assumed to be independent of the filtering process. By applying the test filter on the unresolved SGS scalar flux (${\tau }_{ij}^{SGS}$) and subtracting it from the test filter stress (T_ij), the resolved part of the SGS stress tensor is obtained. This is the Germano identity [40]:

$L_{ij}^{}=T_{ij}^{}-{\widehat{\tau }}_{ij}^{}=\widehat{{⟨u_{ij}⟩}_L{⟨u_j⟩}_L}-{\widehat{⟨u_{ij}⟩}}_L{\widehat{⟨u_j⟩}}_L$E12

The Germano identity can be partially satisfied by replacing model equations for the SGS and filter stresses [Eqs. (9) and (11), respectively] into Eq. (12), giving an approximation for the model coefficient (C_dyn):

$L_{ij}^{}-\frac{1}{3}{L}_{kk}^{}{\delta }_{ij}^{}=-2C_{dyn}\left({\beta }_{ij}^{}-\widehat{{\alpha }_{ij}}\right)=-2C_{dyn}{M}_{ij}^{}$E13

Equation (13) is overspecified, because there are five independent equations and one unknown value. There are different ways to calculate the coefficient. By applying a least‐squares methodology, Lilly [40] minimized the error incurred in the calculation of the coefficient. The error is

$E_{ij}^{}=L_{ij}^{}-T_{ij}^{}+\widehat{{\tau }_{ij}}=L_{ij}^{}-\frac{1}{3}{L}_{kk}{\delta }_{ij}^{}+2C_{dyn}{M}_{ij}^{}$E14

The minimum corresponds to

$\frac{\partial \left(E_{ij}^{}{E}_{ij}^{}\right)}{\partial C_{dyn}}=0$E15

leading to

$C_{dyn}=-\frac{1}{2}\frac{L_{ij}^{}{E}_{ij}^{}}{M_{ij}^{}{M}_{ij}^{}}=0$E16

The calculation of the parameter generates spatial and temporal fluctuations, as well as negative values [43]. Regions with a negative coefficient may be interpreted as regions where backscatter takes place. According to Carati et al. [44], the dynamic Smagorinsky model does not have information about the amount of energy that is available in the subgrid scales.

3. Subgrid‐scale scalar flux models

In the context of mixing of scalars, the interaction between resolved and nonresolved scalar structures is accounted for by the unknown SGS scalar flux term, which must be provided via SGS scalar flux models. In the context of mixing, the aim of LES closures is to express the SGS scalar flux in terms of the known filtered values in order to close the numerical set or partial differential equations. Accounting for the nature of turbulent mixing, where various regimes are present in the scalar spectrum, modeling of the SGS scalar flux models is a complex task that is far from trivial [45].

Macromixing, mesomixing, and micromixing happen simultaneously as the mixing process is taking place. The spatial and temporal transport of large‐scale structures describes the macromixing, and they can be solved in LES. Mesomixing is driven by turbulent fluctuations in the energy‐containing range, and viscous convective deformations of fluid elements and molecular diffusion are responsible mechanisms for micromixing. Because the inertial range of the velocity spectrum is modeled in LES, mesomixing can be expressed by gradient correlations. On the other hand, characterization of micromixing can be based on the second statistical moment (variance) of the local concentration distribution [46–48]. A physical meaning of the concentration variance is that it provides a measure of the scalar distribution from small‐scale homogeneity. The variance production, then, is dependent on the scalar flux.

There are many published works proposing SGS scalar flux models; most of them were developed for gas flow. For review of SGS scalar flux models in gas flows, see Refs. [27, 49]. Even though gas‐based SGS scalar flux models have some limitations on the prediction of mixing where high Schmidt number effects are important [48], physically based models for the SGS scalar flux have been rarely proposed. One exception is the study of Jaberi and Colucci [50].

According to experimental data and DNS results, the assumption of isotropy at inertial and dissipation scales of the scalar field is no longer valid for structured functions and derivative skewness, when a mean scalar gradient is taking place [4, 8, 35]. Different models have been proposed in order to take into account the anisotropic behavior of micromixing, mainly focused on gas‐phase (see Refs. [50–52]), but their application for fluids with high Sc numbers is scarce.

Few models have been developed or used for high Schmidt flows. The eddy diffusivity model [53] with constant and dynamically calculated turbulent Schmidt number [39, 54] has been widely used on the simulation of mixing at high Schmidt numbers [40, 55].

Other models used in the simulation of high Schmidt flows have been adapted from SGS stress tensors or extended from their RANS counterpart. These models include the following:

1. The dynamic mixed model, originally developed by Zang et al. [55] to represent the SGS stress tensor extended to the LES of scalar field by Na [56] in a turbulent channel (Sc/Pr =1, 3, 10); The model was also used by Tkatchenko et al. [57] in the near‐field of a coaxial jet mixer (Sc = 1000).

2. The dynamic structure model proposed by Chumakov et al. [58] and a priori tested using DNS data of a nonreacting mixing layer and decaying isotropic turbulence.

3. The multifractal SGS stress tensor model introduced by [59, 60] and extended to model the SGS scalar flux in a low Schmidt flow [61, 62] in a homogeneous isotropic turbulence at Sc = 100 [15] and a jet mixer [63].

4. The implicit method (non‐SGS model), in which discretization methods are developed, so that the truncation error itself acts as an implicit SGS model. This approach was applied to a stirred tank reactor by Revstedt et al. [32] and isotropic turbulence and channel flow (Sc = 1, 3, 10, and 25) by Hickel et al. [45].

On the other hand, the anisotropy model [48] was proposed for the simulation of mixing of passive and active scalars. Simulation results of mixing at low and high Schmidt regime have shown its superior performance among other SGS scalar flux models. The anisotropy and the eddy diffusivity models are presented in the following subsections.

3.1. The eddy diffusivity model

In analogy to the Smagorinsky model, most of the SGS scalar flux models are based on the eddy diffusivity assumption. The eddy diffusivity model [53] can be viewed as a counterpart of the Smagorinsky model for the scalar field. The eddy diffusivity model relates the SGS flux to the local filtered scalar strain rate, and the transport is assumed to be aligned to the filtered scalar gradient [15]. The proportionality constant is called the turbulent Schmidt number, which appears as an adjustable parameter. However, the turbulent Schmidt number can be dynamically calculated by using the methodology proposed by Germano et al. [39] and further refined by Moin et al. [54].

$J_{\alpha ,j}^{SGS}=-C_{\varphi }{\Delta }^2\left|{⟨S⟩}_L\right|\frac{\partial {⟨\varphi ⟩}_L}{\partial x_j}=-\frac{v_T}{S{c}_T}\frac{\partial {⟨\varphi ⟩}_L}{\partial x_j}$E17

where Sc_T is the turbulent Schmidt number, which is given by the ratio of SGS viscosity to the SGS diffusivity ($S{c}_T=v_T/D_T$).

The eddy diffusivity model has been often used due to its simplicity, and it shows good results for gas flows. The model has been widely used in simple and complex geometries, some involving reactive scalars.

The turbulent Schmidt number appears as an adjustable parameter that can be tuned in order to minimize the error with reference data. According to Durbin and Patterson [63], the value is dependent on the type of flow, and is in the range 0.1–1 [27]. Various authors (see Table 1) have proposed different values for the nondimensional parameter Sc_T in high Schmidt flows.

Configuration	ScT	Reference
Stirred tank; reactors	0.8	[23]
	0.7	[24]; [113], [114]
Channel	0.25	[110]
Inclined channel	0.7	[109]
Confined; impinging jets	0.4	[108]
Coaxial jet	0.7	[107]
Jet in channel	1.0	[48]
Mixing layer	1.0	[48]
Co‐flowing jet	0.7	[99]

Table 1.

Turbulent Schmidt number used in LES of high Schmidt flows.

In turbulent shear flows, different orientations of the mean scalar gradient yield different values for the turbulent Schmidt number [49]. A constant value of Sc_T is, therefore, not adequate. The turbulent Schmidt number can be thought of as a parameter that characterizes the dissipative/diffusive cut‐off scales of the velocity and scalar fields. The behavior of the scalar mixing spectrum, presented in Section 2.2, indicates that a universal distribution of an effective turbulent Schmidt number cannot exist [27].

3.2. The eddy diffusivity model with dynamic procedure

Similar to the dynamic Smagorinsky model, in the dynamic model, the Sc_T in Eq. (17) can be dynamically calculated by using the methodology proposed by Germano et al. [39], which was further extended to scalar transport and compressible flow by Moin et al. [54]. By applying a test filter over the filtered velocity and scalar fields, the turbulent Schmidt number can be computed as follows [64].

The filter scalar flux is

$\frac{1}{S{c}_T}=\frac{1}{2C_{dyn}^{}}\frac{F_i{H}_i}{H_i{H}_i}$E18

where

$F_i=Q_i-\widehat{J_i^{SGS}}=\widehat{{⟨u_i⟩}_L{⟨\varphi ⟩}_L}-{\widehat{⟨u_i⟩}}_L{\widehat{⟨\varphi ⟩}}_L,H_i=\widehat{{\Delta }^2\left|{⟨S⟩}_L\right|\frac{\partial {⟨\varphi ⟩}_L}{\partial x_i}}-{\widehat{\Delta }}^2\left|{\widehat{⟨S⟩}}_L\right|\frac{\partial {\widehat{⟨\varphi ⟩}}_L}{\partial x_i}$E19

Configuration	Sc; (Pr)	Reference
Channel	0.1–200	[105], [106]
	1–100	[104]
	0.1–100	[103], [104]
Mixing layer	600	[10]
	600	[34]
	3300	[48]
Coaxial jet	1000	[57]
Jet in channel	3300	[48]
Grid‐generated turbulent flow	600	[102]
Co‐flowing jet	2000	[98], [99]; [102]

Table 2.

Previous works of turbulent mixing of Sc >> 1 using the dynamic procedure.

On implementation of the dynamic procedure, negative values are clipped to zero, and a relaxation process is imposed to instantaneous values.

The dynamic procedure is one of the most popular SGS scalar flux models in LES of high Schmidt flows. Table 2 summarizes the geometry and Schmidt numbers of previous works, using the dynamic eddy diffusivity model.

Theoretical considerations suggest that the dynamic procedure should be adequate for LES, where the scalar fluctuations are resolved but the velocity fluctuations are not, so that SGS kinetic energy transfer takes place [27]. This is the case of low Schmidt flows.

Since the eddy diffusivity model assumes the alignment of the scalar flux with the scalar gradient, this approach does not predict realistic values of the scalar flux components [45, 48]. In addition, errors are introduced, as the model does not account for the different dynamics of the velocity and scalar fields.

3.3. The dynamic anisotropy model

Another recent SGS scalar flux model is the anisotropy model [48]. Contrary to the eddy diffusivity and dynamic models, the anisotropic model [48] accounts for the nonlinear contributions of the SGS to the turbulent scalar flux. If it is assumed that the turbulent flow field (mean velocities and turbulence characteristics) is available from LES using the dynamic‐based Smargorinsky SGS model [39], the anisotropy model is the simplest explicit anisotropic‐resolving algebraic form, and consists of a cubic formulation in terms of the scalar gradients. This model is thermodynamically consistent as it agrees with the irreversibility requirement of the second law of thermodynamics [47, 65]. It combines the linear eddy diffusivity model with an additional term coupling the (deviatoric) SGS stress tensor and the gradient of the filtered scalar field. Pantangi et al. [65] present a detailed analysis:

$J_{\alpha ,i}^{SGS}=-D_T\frac{\partial {⟨\varphi ⟩}_L}{\partial x_j}+D_{dev}{T}_{sgs}{\tau }_{ij}^{SGS(dev)}\frac{\partial {⟨\varphi ⟩}_L}{\partial x_j}$E21

where D_dev is the anisotropy model coefficient. Both D_T and D_dev depend on the invariants of ${\tau }_{ij}^{SGS}$ and $\frac{\partial {⟨\varphi ⟩}_L}{\partial x_j}$. The tensor diffusivity is defined by

$D_{i,j}^{SGS}=-D_T{\delta }_{ij}+D_{dev}{T}_{SGS}{\tau }_{ij}^{SGS(dev)}$E22

The anisotropy model may lead to other models proposed in the literature, according to the modeling level used for the deviatoric part of the SGS stress tensor. Note that the additional term is a measure of how local mixing, dependent on the molecular Schmidt number through the SGS timescale, is influenced by the nonresolved flow structures (see Eq. (9) in Ref. [65]).

If the turbulent flow field (mean velocities and turbulence characteristics) is available from the LES solver using the dynamic‐based Smargorinsky SGS model, the anisotropy model is the simplest explicit algebraic model that solved accounts for the anisotropy found in turbulent flows. In the context defined by the Smagorinsky model, the anisotropy model can be developed by expressing the SGS timescale using the filter size and the subgrid viscosity. The anisotropy model reads:

$J_{\alpha ,i}^{SGS}=-\frac{v_T}{S{c}_T}\frac{\partial {⟨\varphi ⟩}_L}{\partial x_j}+D_{an}{\Delta }^2\overline{S_{ij}^{SGS}}\frac{\partial {⟨\varphi ⟩}_L}{\partial x_j}$E23

The anisotropy coefficient, D_an, can be either specified or calculated by the dynamic procedure as well as the Sc_T. A preliminary evaluation of this procedure was reported in Ref. [47] for chemical liquid flows and for gases in Ref. [46]. In order to compute the anisotropy coefficient, the turbulent Schmidt number must be calculated using Eqs. (21)–(23). Then, the calculation of the anisotropy model coefficient is done as follows:

$D_{an}=(F_i-\frac{C_s}{S{c}_T}{H}_i)H_{2i};H_{2i}=\widehat{{\Delta }^2\left|{⟨S⟩}_L\right|\frac{\partial {⟨\varphi ⟩}_L}{\partial x_i}}-{\widehat{\Delta }}^2{\widehat{⟨S⟩}}_L\frac{\partial {\widehat{⟨\varphi ⟩}}_L}{\partial x_i}$E24

4. The filtered density function (FDF) approach

A second approach in the LES mixing context is to solve the joint probability function of the SGS scalars [66], named filtered density function (FDF), in conjunction with the filtered momentum equations from LES. In this approach a filtered density function (FDF) is used to quantify the probability to find a filtered variable of the flow (Y) in the range (Y*–ΔY/2, Y*+ΔY/2). The fundamental property of the FDF method is to account, in a probabilistic way, for the effects of SGS fluctuations. The advantage of the FDF approach over other methods is that the chemical source term appears in a closed form, so that the turbulent/chemistry interaction can be correctly included. The counterpart of the FDF equation in the RANS context is named probability density function (PDF) equation. The main difference between PDF and FDF is that temporal fluctuations over different flow realizations are characterized by the use of PDF, while the instantaneous subgrid‐scale fluctuations are characterized by the use of FDF [67]. The PDF is the expected value of the FDF in the limit of vanishing filter width [68].

The statistical information of the reactive scalar fields can be obtained explicitly from a transport equation. Pope [66] introduced the mathematical definition of the FDF transport equation. Pope [66] highlighted one of the major advantages of the FDF approach: the chemical reaction term in the FDF transport equation is closed; so, the modeling of this term is no longer required. Drozda et al. [69] present an overview of the state of progress in FDF.

The mass‐filtered mass density function (FDF) is defined as

$F_L\left(\psi ;x,t\right)={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}G\left(x-x\text{'},t\right)\rho \left(x\text{'},t\right)f\text{'}\left(\psi ;x\text{'},t\right)}dx\text{'}$E25

where Ψ is the sample space variable for each composition and $f\text{'}$ is the fine‐grained joint FDF of compositions [16], defined as

$f\text{'}\left(\psi ;x,t\right)\equiv {\displaystyle \prod _{\alpha =1}^{n_{\varphi }}\delta \left({\varphi }_{\varphi }\left[x,t\right]-{\psi }_{\varphi }\right)}=\delta \left(\varphi \left[x,t\right]-\psi \right)$E26

Based on the properties of Dirac delta functions and, after some manipulation, the FDF equation is [70]

$\frac{\partial F_L}{\partial t}+\frac{\partial }{\partial x_j}\left[F_L{⟨u_j|\psi ⟩}_L\right]=-\frac{\partial }{\partial {\psi }_{\alpha }}\left[F_L{⟨\frac{D{\varphi }_{\alpha }}{Dt}|\psi ⟩}_L\right]$E27

The introduction of the scalar balance equation into Eq. (19) gives

$\frac{\partial F_L}{\partial t}+\frac{\partial }{\partial x_j}\left[F_L{⟨u_j|\psi ⟩}_L\right]+\frac{\partial \left[F_L{S}_{\alpha }\right]}{\partial {\psi }_{\alpha }}=-\frac{\partial }{\partial {\psi }_{\alpha }}\left[F_L{⟨-\frac{1}{\rho }\frac{\partial J_{\alpha }}{\partial x_j}|\psi ⟩}_L\right]$E28

The zeroth, first, and second order moments of FDF are

$\begin{array}{l}{\displaystyle \int F_L\left(\psi ;x,t\right)d\psi }={⟨\rho ⟩}_l\\ {\displaystyle \int {\psi }_{\alpha }{F}_L\left(\psi ;x,t\right)d\psi }={⟨\rho {\varphi }_{\alpha }⟩}_l\equiv {⟨\rho ⟩}_l{⟨{\varphi }_{\alpha }⟩}_L\\ {\displaystyle \int {\psi }_{\alpha }{\psi }_{\beta }{F}_L\left(\psi ;x,t\right)d\psi }={⟨\rho ⟩}_l{⟨{\varphi }_{\alpha }{\varphi }_{\beta }⟩}_L\end{array}$E29

.

Conditionally filtered terms of equations are unclosed. The conditional advection term can be modeled by the conventional gradient diffusion [70]:

$\left({⟨u_j|\psi ⟩}_l-{⟨u_j⟩}_L\right)F_L=-{⟨\rho ⟩}_l{D}_T\frac{\partial }{\partial x_j}\left(\frac{F_L}{{⟨\rho ⟩}_l}\right)$E30

where D_T is the SGS diffusion coefficient. The conditional diffusion term in Eq. (30) accounts for transport in the physical space and mixing in the composition space. It has to be modeled by mixing models. The interaction by exchange with the mean model (IEM) or linear mean square estimation (LMSE) [71, 72] is described by Eq. (32) [70]:

$\frac{\partial }{\partial {\psi }_{\alpha }}\left[{⟨-\frac{1}{\rho }\frac{\partial J_{\alpha }}{\partial x_j}|\psi ⟩}_L{F}_L\right]=\frac{\partial }{\partial x_j}\left[{⟨\rho ⟩}_l{⟨D⟩}_L\frac{\partial }{\partial x_j}\left(\frac{F_L}{{⟨\rho ⟩}_l}\right)\right]+\frac{\partial }{\partial {\psi }_{\alpha }}\left[{\Omega }_m\left({\psi }_{\alpha }-{⟨{\varphi }_{\alpha }⟩}_L\right)F_L\right]$E31

where $C_{\Omega }$ is the frequency of mixing within the subgrid, and it can be modeled as

${\Omega }_m=\frac{1}{{\tau }_m}=C_{\Omega }\left({⟨D⟩}_L+D_T\right)/{\Delta }_L^2$E32

where $C_{\Omega }$ is the SGS mixing time constant. The solution of the FDF equation by conventional difference schemes is intractable due to the high‐dimensionality of the equation: the computational cost rises exponentially as the number of scalars increases [73]. To overcome this problem, stochastic methods have been used for the solution of the equation, because the computational cost rises linearly with the number of scalars [74].

The FDF approach also has a modeling requirement, because the conditional convective flux and conditional diffusion terms in the FDF transport equation are unknown. By deducing an equation for the conditional convective flux based on the velocity‐scalar filtered density function formulation, Ref. [75] showed that a gradient transport hypothesis model performs well under many conditions. This model is usually implemented in FDF methods.

On the other hand, the conditional diffusion term accounts for transport of the FDF in the physical space and mixing in the composition space by the action of molecular diffusivity, and it has to be modeled by mixing models. From a theoretical point of view, McDermott and Pope [30] developed a set of desirable properties that an ideal FDF mixing model should fulfill. These properties differ from the corresponding PDF mixing models (see Ref. [16]). Most of the FDF mixing models have been adapted from their RANS counterpart.

A widely used mixing model is the interaction by exchange with the mean (IEM) model [71, 72]. The IEM model has a deterministic origin, initially formulated for PDF methods in RANS. The IEM states that the stochastic particle only interacts with itself, and the rate of change is proportional to the distance to the mean in the scalar space. Equivalently, the model relaxes the solution in the composition space toward the mean value; the relaxation is done over a subgrid‐scale mixing time. The major drawback of the IEM model is that it preserves the shape of the PDF, avoiding the relaxation toward a Gaussian distribution [76], being an unphysical situation. However, as pointed by [77], large‐scale mixing in inhomogeneous turbulent mixing problems also affects the shape of the distribution, and then limitation of the preservation of the distribution becomes less critical.

Numerical solution of the RANS/PDF equations has shown that the IEM model does not have good predictive capabilities (see Refs. [78, 79] for comparison of the accuracy of different PDF mixing models). On the contrary, the IEM model performs reasonably well in FDF/LES simulations. It has been used in most LES/FDF studies (see Refs. [70, 74, 78, 80–84] and others). There are few LES/FDF works in which other mixing models were used.

Mitarai et al. [78] and Olbricht et al. [74] implemented the modified Curl model—MCurl [85], which is a particle interaction model based on Curl's model [86]. In the Curl model, the PDF of the composition field is described by N stochastic particles. Mixing results from the movement of these particles in the composition space. According to the “simple” model, the mixing processes take place randomly between two stochastic particles, and the resulting composition of the new particle in the following time step is taken as the mean of the preceding particles. The Curl model produces a correct decay rate, but it relaxes the PDF to a bell‐shaped curve instead of the Gaussian shape. This problem can be reduced by using improved models. In the modified Curl model of [85], the particle gradually mixes with the other one at a random mixing rate (for each pair of particles), instead of the originally instantaneous mixing assumption.

The Euclidean minimum spanning tree model (EMST) from Ref. [87] was evaluated by Mitarai et al. [78] and Shetty et al. [77]. The EMST is a sophisticated particle interaction model, where local mixing in composition space occurs between the selected group of particles, in contrast with the previous models (IEM, MCurl, etc.). Briefly, a Euclidean minimum spanning tree is constructed on the ensemble of particles which are selected based on a defined state variable. The Euclidean minimum spanning tree is formed by minimizing the length of the spanning tree. Mixing occurs between pairs of particles connected by a common branch of the tree. The EMST model has a superior performance over other PDF mixing models; however, it demands more computational time. In comparison to the IEM and MCurl, the implementation in the numerical solvers is more complicated.

Cleary et al. [88] implemented the multiple mapping conditioning (MMC) mixing model [89]. In the MMC model, a reference space is mapped to the physical space in order to attain local mixing in the composition space. This is done by choosing one or more reference variables (e.g., mixture fraction, sensible enthalpy, scalar dissipation, and others). Mixing occurs between particles that are close in both physical and reference spaces. The main difference of the MMC model with the EMST model is that local mixing in the former model is indirectly enforced in the composition space by assuring localness in the reference space, while it is directly enforced in composition space in the latter model [90]. Therefore, the principle of independence of scalars is assured in MMC. An advantage of the MMC model in FDF is that the averaged joint scalar distribution can be well predicted by using fewer particles in the computational domain than other models. However, the model has two unknown constants, and the reference space has to be solved in LES.

A recent PDF mixing model is the parameterized scalar profile (PSP), developed by Meyer and Jenny [91] and used by Shetty et al. [77] in LES/FDF. In this approach, it is assumed that one‐dimensional effects dominate the dynamics of molecular diffusion at molecular interfaces. The PSP model is based on a parameterization of one‐dimensional scalar profile in high Reynolds number flows. These scalar profiles are composed of scaled sinusoidal shaped sections, which move with the fluid within a reference frame. There are three parameters used by the scalar profiles; so, they become as additional properties of the particles. The latter makes the model very expensive, because three additional equations must be solved.

A mixing model developed in the LES framework is the fractal IEM (FIEM) [77]. The FIEM is a modified version of the IEM in which it is considered that the behavior of the system at the small scales is repeated at the large scales. It is done by splitting the control volume into many smaller control volumes, each one having the same size, allowing strong mixing in the smaller control volumes and further weak mixing in the larger control volume. The model uses two additional empirical parameters that account for the relative importance of the two mixing steps into the global mixing process. These constants are tuned by trial and error.

The performance of different mixing models was evaluated by Mitarai et al. [78]. These authors compared the numerical results of IEM, EMST, and MCurl in RANS/PDF and LES/FDF with DNS data of a gas diffusion flame with one‐step reaction. The predicted filtered temperature was in an acceptable agreement with DNS results; the results from MCurl and IEM models were similar, while a better agreement was achieved with the EMST model. The mean and variance of the mixture fraction were accurately predicted by all mixing models. Their results also showed that LES/FDF approach provides much better results than the RANS/PDF in which marked differences between PDF mixing models were observed.

Shetty et al. [77] compared the predictions with IEM, EMST, PSP, and the FIEM mixing models in LES/FDF when modeling a low Schmidt three‐stream turbulent jet. Comparisons with experimental data were done in the near‐field of the configuration [max. downstream location = 7.2 jet diameters]. The simulation results of the mean and RMS of the radial scalar distribution obtained with all models followed the experimental observations. However, the FIEM results were in better agreement with the mean experimental observations at the largest downstream location. Results showed fast mixing of the scalar of the inner jet and slow mixing of the scalar of the annular jet, but, conversely, this did not happen when the FIEM was used. Large scale motions play an important role in mixing at the unstable shear layers formed between the inner, annular, and co‐flowing streams, mainly due to entrainment processes. In the inner part of this jet, the scalar field may behave similar to the velocity field. Because the filtered velocity field is supposed to be the same for all tests [not showed in the chapter], little variations are expected in the performance of different mixing models, as was the case for the IEM, EMST, and PSP models. This suggests that the performance of the FIEM could be the result of artificial diffusion transport in the calculations.

Another aspect in the development of FDF is the fact that in LES there may be computational zones where the velocity field is locally fully resolved (the DNS limit). The mixing model has to take into account this aspect. McDermott and Pope [30] showed that the latter can be accomplished by using a mean drift term in the particle composition equation, rather than a random walk in physical space. The model developed by McDermott and Pope [30] correctly reduces to DNS in the limit of vanishing filter width and is able to deal with differential diffusion.

On a discrete representation of the FDF at a given time, the particles in the sample space correspond to particles in physical space, and distance between particles in the subgrid volume is small [67]. In a LES/FDF computational domain, the numerical solution of a FDF equation using a particle‐based method may involve 106–108 stochastic particles. Therefore, the scalar resolution increases and may approximate to that of DNS, depending on the Schmidt number.

According to Klimenko [92], in the DNS limit, when the physical distance between particles becomes infinitesimally small (particle‐based solutions of the FDF equation), many mixing models behave essentially the same. Furthermore, recalling the differences between PDF and FDF methods, the modeling requirements are different: in the RANS context, the PDF aims to solve mean scalar values, and fluctuations about the local mean value must be represented by the mixing models; contrarily, spatially filtered scalar values are sought in the LES approach, and the role of mixing models is to deal with SGS fluctuations about the local spatially filtered scalar. For these reasons, it is expected that simple mixing models used by the FDF method can provide a good approximation [68, 93]. In contrast, RANS has higher levels of sensitivity on the mixing model selection; so, the closure of the conditional diffusive term is the main difficulty in the application of PDF methods [88].

For review of the FDF method, see Haworth [93] who presents a comprehensive paper about theoretical aspects, physical models, numerical algorithms, and applications of PDF‐FDF methods. Drozda et al. [69] present an overview of the state of progress of LES/FDF methods applied to the particular case of turbulent combustion.

The application of the FDF method to high Schmidt flows is scarce. Schwertfirm et al. [94–96] used the FDF approach to study high Schmidt number flows using the DNS of the flow field (a priori analysis). Schwertfirm et al. [94–96] showed for three different Schmidt numbers (Sc = 3, 25), that if a correct definition of the SGS mixing frequency is done then the statistical and instantaneous behavior of the scalar field can be well predicted by using the IEM mixing model. One important finding is that the SGS scalar variance transport can be neglected in the SGS mixing frequency definition.

Van Vliet et al. [97] applied the LES/FDF approach to a tubular reactor with a moderated Reynolds number (Re = 4000) and high Schmidt number (Sc = 2000). They used IEM mixing model. They compared the mean and variance concentration distribution in the axial direction with a few experimental data, and without reaction (conserved scalar). From a qualitative point of view, simulation results are in agreement with experimental observations. However, the quantitative comparison of simulation results with experimental information shows that the simulated concentration (mean) decays faster. A similar behavior was exhibited by the scalar variance, showing a poor agreement with experiments. The discrepancy between experiments and simulations indicates that the mixing rate is overpredicted due to an inaccurate estimation of the subgrid‐scale mixing time, rather than a feature of the IEM model. The performance of the reactor was studied by the LES/FDF approach by changing the Damköhler number over eight (8) orders of magnitude. A LES/FDF numerical study of a low‐density polyethylene tubular reactor was done by Van Vliet et al. [81], at the vicinity of the initiator injection point. Multiple reactive scalars were solved (concentration and temperature), and a mechanism of six chemical reactions was used.

Experimental studies indicate that in gas flows, the subgrid mixing time constant is not universal, having values between 0.6 and 3.1 in the RANS context [76]. Different values of the SGS mixing time constant have been reported, mostly in the range 2–10 for gas flows. For liquid flows, Van Vliet et al. [97] used a value of 3. Mejía et al. [98] evaluated four values (between 4 and 10) in a turbulent round jet. They found that as the subgrid mixing time constant increases, the jet decay rate as well as the mixing rate decreases. The nonuniversality of this constant can affect the predictive capability of LES.

The works of Van Vliet et al. [81, 97] and Schwertfirm et al. [94–96] demonstrate the feasibility of performing (advanced) numerical simulations of systems of practical interest (e.g., industrial reactors), with minimum assumptions and complex chemistry. The numerical simulation is a valuable tool for designing, optimizing, and evaluating such processes. Their simulations showed that the calculation of the SGS mixing time, a common parameter in many mixing models, is an Achilles heel in their simulation of the high Schmidt flow, and it may be also the case for other high Schmidt systems.

5. Comparison between different SGS scalar flux models and the FDF method

In previous works [98, 99], we simulated a turbulent round jet (Re = 10,000) discharging diluted rhodamine B in a co‐flowing stream of water (Sc = 2000), using large eddy simulation. The flow configuration is detailed by Antoine et al. [100]. Three different models for the unknown subgrid‐scale (SGS) scalar flux term were used for closing the filtered scalar balance equation: eddy diffusivity model with both constant turbulent Schmidt number (Sc_T = 0.7) and dynamically calculated turbulent Schmidt number, and the dynamic anisotropy model. In addition, the filtered density function (FDF) method was implemented. The interaction by exchange with the mean (IEM) mixing model was used for closing the conditional diffusion term in the transported FDF equation. The FDF transport equation was solved using a Monte Carlo method.

Figure 1.

Radial mean concentration distribution across the jet. Experimental data (symbols) from Ref. [100]: + = 70; ● = 80; × = 90. Simulation results (lines) from Refs. [98, 99]: Solid line: eddy; dashed line: dynamic; dotted line: anisotropy; dashed‐dotted line: FDF.

Figure 2.

Radial mean (a) and fluctuation (b) concentration distribution across the jet. Experimental data (symbols) from Ref. [100]: + = 70; ● = 80; × = 90. Simulation results (lines) from Refs. [98, 99]: Solid line: eddy; dashed line: dynamic; dotted line: anisotropy; dashed‐dotted line: FDF.

Figures 1 and 2 compare the simulated radial mean and fluctuation distribution of the scalar concentration [98, 99], normalized by the centerline value (η = r/x, where r is the radial coordinate and x is the downstream distance from the nozzle exit), with the experimental data [101].

Having in mind the limitations of the eddy diffusivity model, it is noted in Figure 1 that the mean concentration distribution is predicted reasonably well. On the other hand, Figure 2 reveals that the three SGS scalar flux models fail to reproduce the concentration fluctuations at higher values of the nondimensional radial coordinate (η). Although the LES simulation results of the concentration distribution in radial direction are similar for the tested SGS scalar flux models in Figures 1 and 2, the dynamic anisotropy model provides a better performance. In LES, a mesh grid refinement would improve the simulation results. Mejía et al. [98] also showed that if the mesh grid is further refined, the dynamic eddy diffusivity improved its predictive capabilities.

Figure 3.

Streamwise velocity‐concentration correlation across the jet. Experimental data (symbols) from Antoine et al. [100]: + = 70; ● = 80; × = 90. Simulation results (lines) from Refs. [98, 99]: Solid line: eddy; dashed line: dynamic; dotted line: anisotropy; dashed‐dotted line: FDF.

Figure 4.

Radial velocity‐concentration correlation across the jet. Experimental data (symbols) from Antoine et al. [100]: + = 70; ● = 80; × = 90. Simulation results (lines) from Refs. [98, 99]: Solid line: eddy; dashed line: dynamic; dotted line: anisotropy; dashed‐dotted line: FDF.

The simulation results of the FDF approach are in good agreement with the mean and fluctuations concentration distribution for all regions of the jet. The radial distributions of the axial and radial scalar fluxes are presented, respectively, in Figures 3 and 4.

Figures 3 and 4 show that the dynamic eddy diffusivity and anisotropy models reproduce the behavior of the scalar fluxes for low and intermediate values of the radial coordinates.

The agreement of the simulations of the FDF method (Figures 3 and 4) with the experimental data is good, and much better than that obtained with LES simulations using advanced SGS scalar flux models, such as the anisotropy model (see, e.g., Figures 1 and 2 ). The LES/FDF can capture strong intermittency effects that take place on the surroundings of the superviscous layer formed between the jet and the co‐flowing streams, leading toward anisotropy of the scalar field. Contrary to the SGS scalar flux models evaluated in [98, 99], the FDF equation does not spread the scalar away from the superviscous layer, as seen in the predicted rms and scalar fluxes distributions.

The computational costs of the dynamic and anisotropy models were 1.2 and 1.7 times that of the eddy diffusivity model, respectively. When a higher level model has to be considered, it is usual that more computational resource is required. From an engineering point of view, the dynamic model captures the most important mixing characteristics at a rather low computational cost. The computational time of the FDF, compared to the one required by using the eddy diffusivity model is approximately 1.9 times higher though [101].

6. Conclusions

Large eddy simulation of transport and mixing in high Schmidt flows remains a challenge, in particular for engineering applications. The appraisal of conventional and advanced SGS scalar flux models as well as the FDF method has to done based on the object of study. However, some of the merits of the approximations can be highlighted. The SGS scalar flux models have the advantage over FDF methods in that the simplest ones (e.g., eddy diffusivity model) are available in many CFD solvers. This closure strategy uses the same mesh grid of the flow solver as well as the numerical methods used for the numerical solution of the N‐S equations. Therefore, one important consequence is that the programming of advanced SGS scalar flux models is straightforward. Having in mind the computational cost, simple SGS scalar flux models can reproduce important characteristics of the scalar field such as mean quantities on a moderate grid resolution. However, if the application requires a deeper understanding of mesomixing and micromixing (such as combustion and flow instability studies), the use of advanced SGS scalar flux models becomes mandatory. On the other hand, the main advantage of FDF method is that it allows for a detailed description and simulation of the scalar field. The FDF transport equation can be solved using particle‐based methods, which are very simple to implement and to couple to LES solvers. The computational domain can be simpler than the one used by LES, since it only has to account for the spatial domain where the scalar field evolves.

Acknowledgments

We thank Deutscher Akademischer Austausch Dienst (DAAD) and Departamento Administrativo de Ciencia, Tecnología e Innovación de Colombia (COLCIENCIAS) for logistic and financial support. The first author acknowledges the Institute of Energy and Power Plant Technology (EKT) at the Technische Universität Darmstadt (Germany) for its support during a 1‐year academic stay. We thank the reviewers for their helpful comments on an earlier version of the manuscript.