Mixing mechanisms in the view of mixing indicators: from passive-scalar mixing to variable-density mixing

Nomenclature

1.   Introduction

1.1.   Engineering background

Efficient mixing is essential for maintaining stable and efficient combustion in air-breathing power systems. The characteristics of mixing are critical determinants of key combustion factors, such as combustion efficiency, flame extinction, and pollution emissions. Despite significant advances in ramjet and gas turbine engines, achieving mixing enhancement remains a significant challenge in these two advanced representative engines.

In the case of a ramjet engine, the environment is characterized by a high-speed incoming flow ($Ma = 2 \sim 6+$), and achieving stable and efficient combustion in this setting is a critical challenge ^[3]. While the feasibility of supersonic combustion was theoretically analyzed in 1964 ^[4], Ferri further noted in 1973 that efficient supersonic combustion is mixing-controlled ^[5]. However, a unified design methodology for organizing the flow in the ramjet combustion chamber to achieve efficient mixing has not been established ^[3].

The primary focus with turbo-based engines is to reduce fuel consumption and carbon emissions. Over nearly half a century of exploration, the design of turbo-based combustion chambers has been based on swirling flow ^[6]. The approach is characterized by using the recirculation region to stabilize the flame and enhance mixing ^[7]. To achieve significant thrust, a high temperature rise is required, which necessitates a high fuel-air ratio in the combustion chamber. Consequently, additional intake air should be introduced into the primary zone, as discussed by Bahr ^[8]. In that case, the combustion chamber utilizes a strong swirl flow to mix these additional fuel and intake air, that may result in flame extinction caused by excessive mixing in some lean fuel conditions ^[9].

The design of these two advanced power systems demonstrates that while ramjet combustion is mixing-controlled, there is currently no established methodology for organizing the flow inside the combustion chamber to enhance mixing. Although the turbo-based engine primarily relies on the use of swirling flow, problems such as flame extinction resulting from excessive mixing still persist. The practical engineering design problems in these two systems can be largely attributed to the need for a better understanding of the mixing mechanisms and enhancement methods in the combustion chamber.

1.2.   Scientific problem

From a more fundamental perspective, mixing is a general and fascinating flow behavior that appears in numerous fields, including ocean flow ^[10], biological reproduction ^[11], virus diffusion ^[12], and natural/forced heat convection ^[13]. Accordingly, various schools of mixing theory have emerged. An important branch of mixing mechanism can be traced back to 1958 when based on the Kolmogorov turbulent cascade theory ^[14], Batchelor ^[15] proposed the existence of a diffusion scale smaller than the Kolmogorov scale at low diffusion coefficients. This type of mixing is called passive-scalar (PS) mixing as it is passively affected by the flow ^[16]. This theory paved the way for large eddy simulation (LES) and was later verified by subsequent direct numerical simulation (DNS) results ^[17]. It is commonly accepted that turbulent flow promotes mixing. Therefore, a criterion to predict the transition of turbulent mixing is crucial. It was not until 2000 that Dimotakis proposed a scale separation physical concept ^[18] and put forward the criteria for the onset of turbulent mixing. In recent years, Villermaux has provided theoretical solutions to the advection-diffusion equation based on simple flows such as single vortex flow ^[19] and shear flow ^[20]. Moreover, Villermaux also pointed out that stretching on moving lamellas enhances diffusion, confirming that turbulent flow may exponentially increase the mixing indicator. However, some researchers recently suggested that turbulence may not be the optimal mixing organizer at low diffusion coefficients, i.e., large Péclet numbers ($Pe$) ^[21]. In some specific cases, laminar mixing may actually be more efficient than turbulent mixing ^[22]. Overall, the complexity of mixing mechanisms results from the diversity and intricacy of flows, the analytical difficulty of solving equations, and the disagreement among views. Mixing can only be analytically described in simple flows and is still far from the controllable stage.

Regarding the internal flow of engine combustor, a significant density change occurs. In 2005, Dimotakis categorized mixing into three levels ^[23]. The first level is passive-scalar (PS) mixing, where the mixing process does not affect the flow. The second level is variable-density (VD) mixing, characterized by mixing in fluid with varying densities. The misalignment of the density gradient and the pressure gradient generates baroclinic vorticity, altering the local flow field. The third level is combustion mixing, characterized by the flame changing quantities such as temperature, density, pressure, viscous coefficient, etc. Although the second and third levels of mixing are the main types that occur in actual combustion chambers, Dimotakis ^[23] concluded in the Annu. Rev. Fluid Mech. that, "Unfortunately, however, there is no (even empirical) theory to describe the variable density mixing problem." The lack of connection between basic theory and practical problems underlies the urgent need for more theoretical guidance in the design of mixing enhancement devices, which is still mainly based on experience.

This review is centered around exposing the mixing mechanisms in PS mixing and VD mixing from the perspective of mixing indicators. The fundamental step in mixing research is defining objective mixing indicators. The two commonly utilized indicators can be represented by $c^i$, where $c$ represents concentration, and $i = 1$ denotes first-order statistics, while $i = 2$ stands for second-order statistics ^[24]. The following sections are organized as follows. Section 2 focuses on PS mixing by introducing the first-order and second-order statistics used to describe the mixing phenomena and the corresponding mixing mechanisms discovered through related mixing indicators. Section 3 is devoted to VD mixing, discussing new mixing indicators and phenomena and reviewing the unique characteristics of VD mixing compared with PS mixing. Finally, Section 4concludes this review and points out future research directions.

2.   Passive-scalar mixing mechanisms

PS mixing, according to the classification proposed by Dimotakis ^[23], is the first-level mixing characterized by the mixing process not affecting the flow. Research on PS mixing mainly focuses on low-velocity flows, thus, PS mixing can be described by the advection-diffusion equation:

where $\mathscr{D}$ is the diffusion coefficient. There is also a restriction condition that the velocity divergence is 0, given by

Therefore, how does PS mixing occur, and what is the mechanism behind it? To answer this, we can illustrate the mechanism by examining the mixing process in particular flows, as shown in Figure 1.

The mixing processes depicted in Figure 1 are single-vortex mixing ^[19], random flow mixing ^[25], jet mixing ^[26], and coaxial jet mixing ^[27], respectively. These mixing processes demonstrate that the initially clustered high-concentration region is stretched into various structures, including lamellas and bubbles. During the stretching and elongation of the lamellas, the scalar concentration on these structures decreases until it is conserved with the surrounding fluid. On the bubble structures, the concentration remains virtually unchanged. These examples underscore that the mixing process can be influenced by stretching and that the mixing dynamics are distinct on different kinds of flow structures. Describing these phenomena has always been the primary focus of mixing research. As mentioned earlier, mixing indicators can be categorized into two groups: first-order statistics and second-order statistics. From the following review, one can find that the first-order statistics describe the mixing mechanisms on the lamella structures, and the second-order statistics describe the mixing mechanisms on the flow structures, which is described by the invariants of velocity gradient tensors and the rotation of principal strain axis.

2.1.   First-order statistics: mixing mechanisms on the lamella structures

The first type of mixing indicators is the first-order statistics, which typically refers to the scalar concentration $c$. This theory focuses on the concentration distribution $p(c)$ on the lamella structures and aims to establish the relationship between the evolution of the concentration and lamella stretching. Concentration is a vital physical quantity in various mixing scenarios. For instance, factories need to monitor whether the concentration of pollutants in wastewater exceeds the standard level. In 1989, Vidick observed that the residue of additives in glass or cement could lead to mechanical defects ^[28]. Poulain discovered that the lifespan of the liquid film, such as bubbles in seawater, is determined by the concentration of impurities ^[29]. These examples underline the importance of describing the distribution of concentration $c$ over time and space and its probability density distribution $p(c)$.

The theory of first-order statistics was firstly introduced by Marble in 1985 to explore the relationship between product generation and circulation in the single vortex flame ^[30]. Since Marble's proposal, using this theory to analyze PS mixing has drawn significant attention, and Villermaux further developed this approach, which has now become a popular research paradigm. Villermaux explains this method in detail in his publication ^[31], and only a brief outline of the framework will be presented here.

In this theoretical framework, a Lagrangian coordinate system is established on a two-dimensional moving lamella, as shown in Figure 2. The $x$ axis points in the direction of the largest scalar gradient, while the $y$ axis is perpendicular to it. The width of the lamella is denoted by $s(t)$, and the length of the lamella is denoted by $l(t)$. The velocity in the $x$ direction is denoted by $u_1$, and the velocity in the $y$ direction is denoted by $u_2$. Given that $\frac{\partial u_j}{\partial x_j} = 0$, the velocity $(u_1, u_2)$ can be associated with the stretching of this lamella using

where $\dot{s}$ represents the time derivative of the width $s(t)$.

The advection-diffusion equation on the lamella is given by $\frac{\partial c}{\partial t} + (\dot{s}/s)x\frac{\partial c}{\partial x} = \mathscr{D} \frac{\partial^2 c}{\partial x^2}$, where a Ranz transform ^[32,33] is applied to convert it into the Fourier equation for analytical solutions, given as $\xi = \frac{x}{s(t)}, \tau = \mathscr{D}\int_{0}^{t}\frac{dt'}{s(t')^2}$. Thus, we have:

Eq 2.4 describes the corresponding relationship between the concentration $c$, spatial position $\xi$, and stretch time $\tau$ to obtain the probability density function $p(c)$. The solution of Eq 2.4 is

When considering the interaction between two lamellas, the linearity of the Fourier equation allows for direct addition of the concentration, giving:

As the result the direct addition of the concentration in Eq 2.6, the probability density function can be calculated using convolution as

This method is described in-depth in the literature ^[34,35].

In this theoretical framework, two important mixing indicators can be defined: mixing time and Batchelor scale. From the solution of Fourier equation Eq 2.4 as described in Eq 2.5, $x = 0$ is the position with the highest concentration, which can be expressed as

When $\tau \ll 1$, $c_{max}(\tau)$ approaches 1, whereas when $\tau \gg 1$, $c_{max}(\tau) = \frac{1}{\tau}$ decreases with time. The distribution of concentration on the lamella is Gaussian ^[31]:

where the spatial scale is given by:

If the lamella is stretched exponentially, $s(t) = s_0 e^{-\gamma t}$, then $\eta(t)$ is constant, indicating that the concentration distribution is no longer affected by the flow. Therefore, the time $t_s$ satisfying the condition $\tau(t_s) \sim O(1)$ is defined as the mixing time. For instance, the mixing time of PS mixing in a single vortex can be expressed as:

where $\varGamma$ is the circulation of the vortex, $\mathscr{D}$ is the diffusion coefficient, and $R_0$ and $s_0$ are the shape factors ^[19]. The mixing time signifies that after time $t_s$, mixing transitions from the stretch-dominant stage to the diffusion-dominant stage. The scale $\eta(t)$ corresponds to the Batchelor scale, which has a similar physical meaning as proposed in Batchelor's turbulent mixing theory ^[15].

Villermaux et al. have applied this theoretical framework to study three types of problems:

The first type of research involves establishing stretching models in different flows and examining the corresponding concentration $c$. For example, in 2003, Meunier et al. studied the concentration distribution and mixing time in a single vortex flow ^[19]. Villermaux et al. investigated mixing in jet flow and obtained its $p(c)$ function ^[34]. Duplat et al. examined the mixing behavior in random and turbulent flows. Duplat concluded that the concentration follows a gamma distribution in random flows ^[25,35], and highlighted the self-convolution behavior resulting from the interaction between lamella ^[36].

The second type of research aims to investigate more accurate ways to describe the stretching indicator. In 2019, Gotzfried et al. examined whether the Finite-time Lyapunov exponent (FTLE) proposed by Haller ^[37] could be used as the stretching indicator for the concentration distribution of turbulence. FTLE has been used in other areas such as vortex ^[38], mixing layer ^[39], and vortex interaction ^[40]. The concentration distribution predicted by using FTLE as the stretching indicator was consistent with the results of DNS ^[41]. In 2022, Meunier attempted to describe the stretching behavior of "diffuselet" accurately in three-dimensional flows by using the velocity gradient tensors or the FTLE number ^[42].

The third type of research is focused on predicting scalar distribution using stretching information of the flow without using computational fluid dynamics (CFD) methods. In 2010, Meunier and Villermaux proposed the Diffusive Strip Method (DSM), which accurately predicted the concentration distribution of a single vortex system or sine flow ^[43]. In 2018, Martinez-Ruiz et al. extended this method from two to three dimensions, resulting in the three-dimensional Diffusive Sheet Method, which accurately predicted scalar distribution in Taylor-Couette flow ^[44]. Recently, Sen et al. used DSM to analyze the dynamics of solute mixing and reaction in a mixing-limited reactive flow in both shear flow and Rankine vortex cases ^[45].

To summarize, the first-order statistics theory provides an analytical solution to the Fourier equation using the Ranz transform, which accurately describes the relationship between scalar distribution and lamella stretching on high Péclet number lamellas. This theory effectively describes the PS mixing mechanisms on lamella structures. Recent research on first-order statistics has focused on two frontiers for studying mixing: 1) This theoretical framework is extended to reactive flows, such as mixing-limited reactive flow in linear shear and Rankline vortex cases ^[45], porous media flow ^[46], and reactive flow with different Damköhler numbers ^[47,48]. 2) The stretching indicator was given by modeling in the research before. Recently, discussions have focused on developing a more fundamental stretching indicator using basic physical quantities in the flow field such as velocity and pressure ^[42].

2.2.   Second-order statistics: mixing mechanisms on flow structures

Second-order statistics in mixing research typically include mixedness $f$ and scalar dissipation rate $\chi$. Mixedness is a measure of the degree of mixing and is expressed as a second moment of concentration. Depending on the specific research problem, different forms of mixedness may be used. For example, Cetegen proposed the mixedness given by Eq 2.12 in 1993 to study laminar flow single-vortex mixing ^[1].

This definition assumes $f = 1$ at $c = 0.5$. Another commonly used form of mixedness is $f_s$ in Eq 2.13,

which was introduced by Verzicco in 1995 for combustion problems ^[2]. $f_s$ assumes $f_s = 1$ when the concentration $c$ is equal to the stoichiometry concentration $c_s$.

Scalar dissipation rate $\chi$ is defined as the dot product of the scalar gradient $\nabla c$ ^[49]

If mixedness is defined as Eq 2.12, the relationship between the mixedness and scalar dissipation is

Eq 2.15 shows that the scalar dissipation rate is related to the rate of change of mixedness with respect to time, so this quantity is an indicator of the mixing rate.

Mixedness and scalar dissipation rate have significant importance in the field of mixing research.

1) These two quantities objectively evaluate the degree of mixing and mixing rate. By using mixedness and scalar dissipation rate, researchers have discovered new mixing phenomena. For instance, Cetegen found a positive linear correlation between mixedness and time and vortex intensity for single vortex and vortex pairs ^[1,50]. Basu further developed an empirical relationship between mixedness and temperature, circulation, and interaction time ^[51].

2) The scalar dissipation rate provides a crucial link between flow and mixing. Researchers have utilized the scalar dissipation rate to gain insights into the relationship between flow and mixing. One representative study is that Flohr ^[52] investigated mixing cascade affect by a given single vortex by analyzing the evolution of scalar dissipation rate in the spectral space, and this work demonstrated the usefulness of scalar dissipation rate in understanding the relationship between flow and mixing. In this section, the studies on this relationship will be reviewed in detail.

3) Scalar dissipation rate is a fundamental physical quantity in combustion theory that is highly relevant to various combustion factors. In non-premixed combustion, the scalar dissipation rate positively correlates with combustion efficiency ^[53,54]. In premixed combustion, the scalar dissipation rate affects flame speed and thickness ^[55,56]. Consequently, the scalar dissipation rate has been the subject of much research in the field of combustion, with numerous studies examining this quantity in the context of practical combustion problems ^{[57,58,59,60,61,62]}.

In this section, we will review studies that have explored the dynamic properties of scalar dissipation rate and how this physical quantity can be used to establish the relationship between the flow structure and mixing.

As shown in Figure 3, the starting point of relating flow and mixing is the scalar dissipation rate equation. Buch et al. ^[49] derived the transport equation for scalar dissipation rate for the first time

where ${\boldsymbol{S}}$ is the strain rate tensor:

In Eq 2.16, the first term on the right side is called the "stretching term". This term describes how the fluid stretches the scalar gradient, thus increasing or decreasing it. In turbulent mixing and combustion research, the stretching term is often called "turbulence-scalar interaction" ^[68,69], which is the core of studying the kinematic properties of the scalar dissipation rate equation. The second term on the right is always negative, indicating that the scalar dissipation rate decreases due to the effects of diffusion.

The stretching term in the scalar dissipation rate equation, $(\nabla c \cdot {\boldsymbol{S}} \cdot \nabla c)$, bears a resemblance to the stretching term of the element line stretching equation in incompressible, inviscid flow ^[65], given by

where ${\boldsymbol{l}}(t)$ is the element line vector. When this element line is stretched, its width decreases, leading to an increase in the scalar gradient on this element line. The similarity between the stretching term in the scalar dissipation rate equation and the element line stretching equation is thus evident. The stretching term also occurs in the enstrophy equation ^[70], given by

where ${\boldsymbol{\omega}}$ represents the vorticity vector of the flow.

By comparing the equations for scalar dissipation rate (Eq 2.16), element line (Eq 2.18), and enstrophy (Eq 2.19), it can be observed that the stretching terms of these equations have a similar mathematical form, except for the opposite sign in the scalar dissipation rate equation. This suggests that stretching is an essential fluid characteristic which occurs not only in the evolution scalar dissipation rate but also in the element line stretching and vortex stretching phenomena. Thus, studying the mechanisms of element line and vortex stretch can contribute to a deeper understanding of the scalar dissipation rate.

The stretching term $\nabla c \cdot {\boldsymbol{S}} \cdot \nabla c$ can be expressed in tensor form as $\nabla c_i S_{ij} \nabla c_j$. This term can also be written in vector form as $\nabla c^T {\boldsymbol{S}} \nabla c$, where $\nabla c = [\nabla c_1, \nabla c_2, \nabla c_3]^T$ and ${\boldsymbol{S}} = [S_{ij}]$. Since ${\boldsymbol{S}}$ is a real symmetric matrix, one can perform an orthogonal decomposition such that ${\boldsymbol{S}} = X^T {\boldsymbol{D}} X$, where ${\boldsymbol{D}} = diag\{ s_i\}$ is a diagonal matrix, $X = [{\boldsymbol{\xi}}_1, {\boldsymbol{\xi}}_2, {\boldsymbol{\xi}}_3]^T$ is an orthonormal basis, and ${\boldsymbol{\xi}}_i$ indicates the direction of the principal strain axis. Here, $s_i$ represents the principal strain. Without loss of generality, one can define $s_1 > s_2 > s_3$, where $s_1, s_3$ are principal strains in the stretching and compression directions, respectively. If $\nabla \cdot {\boldsymbol{U}} = 0$, then $s_1+s_2+s_3 = 0$. Expressing $\nabla c^T {\boldsymbol{S}} \nabla c$ using the orthogonal decomposition yields,

Here $\lambda_i = \frac{\nabla c}{|\nabla c|} \cdot {\boldsymbol{\xi}}_i$ represents the cosine of the angle between the scalar gradient $\nabla c$ and the principal strain axis ${\boldsymbol{\xi}}_i$. The principal strains $s_i$ and the angle $\lambda_i$ determine the stretching term.

Early research on the stretching mechanism focused mainly on the stretching of the element line. However, due to limitations in experimental and numerical methods, the angle $\lambda_i$ could not be directly measured, and thus, the first-stage research on the stretching mechanism was based primarily on hypotheses, as shown in Figure 3. One such hypothesis was proposed by Batchelor and Townsend. They suggested that the element line vector ${\boldsymbol{l}}(t)$ was aligned with the principal strain in the stretching direction, denoted by ${\boldsymbol{\xi_1}}$, with $\lambda_1 = 1$ and $\lambda_2 = \lambda_3 = 0$.

Townsend first proposed the hypothesis that the element line vector was aligned with the principal strain in the stretching direction in his study of the cooling of heat spots in grid turbulence in 1951 ^[71]. In 1952, Batchelor studied the element line stretching in isotropic turbulence and found that the element line was exponentially stretched. Assuming this hypothesis ^[63], Eq 2.18 can be simplified to:

Here, $\overline{\delta L(t)}$ represents the modulus of the element line vector ${\boldsymbol{l}}(t)$, then the exponential stretch rate $\zeta _{line}$ can be described by

Cocke provided a strict description of the exponential stretching of the element line over a limited period of time in 1969 ^[72]. However, after the 1980s, with the widespread use of CFD technology in fluid mechanics research, numerical results showed that the hypothesis proposed by Batchelor and Townsend was not valid, as represented in the second stage of research shown in Figure 3. For instance, as in Figure 4, using DNS simulations, Ashurst observed that the vorticity tends to align more closely with the intermediate principal strain direction, whereas the scalar gradient direction has a high probability of aligning with other directions ^[73]. Vincent obtained similar conclusions when studying the angle between vorticity and principal strain in isotropic turbulence using DNS in 1991 ^[74]. Girimaji and Pope et al. simulated isotropic turbulence with $Re_T = 90$ to study the element line stretching, and found that the average exponential stretching rate was only one-third of the stretching rate proposed by Batchelor, as given by Eq 2.22 ^[64]. The contradiction between the hypothesis and these phenomena suggests that the element line vector is not consistent with the principal strain direction:

This discrepancy has prompted researchers to investigate which flow properties affect the alignment dynamics.

In 1991, Dresselhaus derived the alignment equation for three-dimensional element line stretching problems for the first time ^[65]. The alignment equation is given by:

$\xi_{line}$ represents the stretching rate of the element line:

${\boldsymbol{\Omega}}$ is the angular velocity:

where $\varepsilon _{ijk}$ is the alternating tensor, and ${\boldsymbol{\Omega'}}$ is the rotation of principal strain axis:

where ${\boldsymbol{\dot S}}$ represents the time derivative of strain rate tensor.

The first two terms of the alignment equation given by Eq 2.24 attract the element line vector to the principal strain axis at an exponential rate, while the third term rotates the element line around the "effective rotation vector" ${\boldsymbol{\mathfrak{R}} } = {\boldsymbol{\Omega - \Omega'}}$.

Dresselhaus emphasized that the "effective rotation" is crucial for the alignment dynamics. However, determining its direction in three-dimensional flows is challenging, adding considerable complexity to the alignment problems. Therefore, simplifying the alignment equation given by Eq 2.24 to obtain the general properties of a particular flow has attracted significant attention, and the two-dimensional stretching theory is an excellent example of simplification and represents the fourth-stage research in Figure 3. The alignment equation Eq 2.24 can be simplified because the direction of ${\boldsymbol{\mathfrak{R}}}$ in a 2D flow is always perpendicular to the plane.

Early studies on the two-dimensional stretching problem were represented by Okubo's study on particle aggregation ^[75], and preliminary conclusions were obtained in Weiss's research on the enstrophy transport represented by Eq 2.19 in two-dimensional turbulence ^[76]. In incompressible, two-dimensional flows, the principal strains have the relationship $s_1 = -s_3, s_2 = 0$. Under the assumption that the material derivative of the strain is negligible, i.e.,

they introduced the parameter $\theta = s_1^2 - \omega^2$, similar to the widely used vortex identification method Q criterion ^[77]. Based on this parameter, they proposed the Okubo-Weiss criterion: in the region of $\theta > 0$, the fluid has hyperbolic characteristics, and the vorticity gradient increases exponentially, while in the region of $\theta < 0$, the fluid has elliptic properties, and the vorticity gradient oscillates periodically.

Following Weiss's proposal, many researchers verified the validity of the hypothesis in Eq 2.28. Representative studies include Basdevant ^[78] and Hua ^[79], who found that Weiss's hypothesis was only valid near the core of the vortex and not at the edge of the vortex. They also discovered that the strain acceleration tensor and pressure Hessian matrix play a crucial role in stretching problems. Similar phenomena were also observed in three-dimensional turbulence ^[80,81,82].

Subsequently, Klein et al. studied the problem of passive scalar gradient in two-dimensional turbulence and defined dimensionless numbers $r$ ^[66] and $n$ ^[67] that determine the two-dimensional stretching properties. They accomplished this by using the geometric model shown in Figure 5 to simplify the alignment equation given by Eq 2.24.

The stretching equation of the scalar gradient without diffusion effect, $(\frac{\partial}{\partial t} + u_j\frac{\partial}{\partial x_j}) \frac{1}{2} (\nabla c \cdot \nabla c) = -(\nabla c \cdot {\boldsymbol{S}} \cdot \nabla c)$, can be transformed using the geometric model shown in Figure 5. The strain and vorticity are defined as:

The scalar gradient and strain are represented by their modulus and orientation:

where $\rho$ is the modulus of the scalar gradient, $\theta$ represents the angle between the scalar gradient and the $x$ axis, and the angle between the principal strain and $x$ axis is $-\phi - \frac{\pi }{4}$.

Using this geometric model, Eq 2.24 can be simplified,

where $\frac{1}{2}\zeta + \frac{1}{4}\pi$ is the angle between the scalar gradient and the principal strain axis in compression direction.

Klein et al. defined a dimensionless number $r$ to discuss the two-dimensional alignment equation given by Eq 2.31,

where the numerator part $\omega+ 2({\rm d}\phi / {\rm d}t)$ represents the sum of vorticity and rotation of the principal strain axis, which is a simplification of the "effective rotation" ${\boldsymbol{\mathfrak{R}}} = {\boldsymbol{\Omega - \Omega'}}$ in the two-dimensional case. The denominator part $\sigma$ is the modulus of the principal strain. Therefore, this dimensionless number measures the ratio of "effective rotation" to the principal strain. If $r^2 < 1$, Eq. 2.31 has a stable point, and the flow is hyperbolic. The scalar gradient vector will converge to the stable direction. If $r^2 > 1$, Eq 2.31 is unstable, and the flow is elliptical. The scalar gradient will not converge and rotate around the $z$ axis.

The second-order time-derivative equation ${\rm d^2}(\nabla c \cdot \nabla c)/{\rm d}t^2 = -\nabla c \cdot {\boldsymbol{N}} \cdot \nabla c$ can also be transformed using the same geometric model. Here, $N$ involves ${\rm d}{\boldsymbol{S}}/{\rm d}t$ and is directly relate to pressure Hessian, the specific expression of which can be found in Klein's publication ^[67]. Another dimensionless number $n$ can be defined, given by

which measures the time change rate of the principal strain. This dimensionless number determines the alignment dynamics in the elliptic region. The derivative of the scalar gradient angle is the smallest when it aligns with the eigenvector ${\boldsymbol{N_-}}$ of matrix ${\boldsymbol{N}}$, which ${\boldsymbol{N_-}}$ depends on the dimensionless number $r$ and $n$. The above conclusions were confirmed by Klein's observations in two-dimensional turbulence ^[66,67].

Based on the discussion of the alignment dynamics of the element line (Eq. 2.31) and the scalar gradient (Eq. 2.24), one can conclude that the alignment dynamics are influenced by the flow properties, which are determined by the velocity gradient tensors that contain two parts: the strain rate tensor in Eq 2.17 and the vorticity in Eq 2.26 ^[84], and they are closely related to the "flow structure" ^[77,85]. One way to study the "flow structure" is through the application of velocity gradient tensor invariants, which define flow topologies ^[84,86,87]. Different flow topologies can be defined based on the eigenvalue properties of the velocity gradient tensor matrix:

Flow topologies with focus, UF/C (unstable focus/compressing), SF/ST (stable focus/stretching), SF/C (stable focus/compressing), and UF/ST (unstable focus/stretching), are dominated by rotation, while flow topologies with saddle points, UN/S/S (unstable node/saddle/saddle), SN/S/S (stable node/saddle/saddle), SN/SN/SN (stable node/stable node/stable node), and UN/UN/UN (unstable node/unstable node/unstable node), are dominated by stretching/compression, as shown in Figure 6. While scalar dissipation rate has been investigated in different flow topologies in compressible turbulence ^[88], shock turbulence interaction ^[89], premixed flame ^[83,90] and non-premixed flame ^[91], the alginment dynamics cannot be determined solely by the invariants of velocity gradient tensors from Eq. 2.24. The rotation of the principal strain axis also plays a vital role in influencing the alignment dynamics of the scalar gradient, which is often neglected in flow structure research. To capture the full dynamics of the alignment process, it is necessary to consider both the velocity gradient tensor invariants and the rotation of the principal strain axis simultaneously. One of the frontiers in mixing research is to establish the complete stretching dynamic equation to explore the relationship between mixing and "flow structure" by scalar dissipation rate.

Regarding the second-order statistics, significant progress has been made in studying mixedness and scalar dissipation rate.

This section mainly reviews the studies that aim to establish a deep understanding of flow and mixing through scalar dissipation rate, by examining the dynamic properties of the stretching term in Eq 2.16. These studies mainly focus on the alignment dynamics of the scalar gradient. The research history can be divided into four stages: hypothesis ^[15], experimental verification ^[73], theoretical derivation ^[65], and equation simplification ^[66]. It has been shown that scalar dissipation rate exhibits different dynamical properties in flows with different flow characteristics. To gain a more thorough understanding of the relationship between "flow structure" and mixing by scalar dissipation rate is a crucial frontier in mixing research.

2.3.   A brief summary on passive-scalar mixing

The PS mixing phenomena are described by two types of mixing indicators: first-order statistics and second-order statistics. These indicators reveal PS mixing mechanisms from different levels.

First-order statistics describe the PS mixing mechanism on the lamella structures. This theory solves the Fourier equation to obtain concentration distribution $p(c)$, and then defines two significant mixing characteristic parameters: the Batchelor scale and mixing time ^{[25,30,31,34,35,36]}. At present, the first-order statistics theory has two priorities: 1) applying this theory to combustion flows ^{[45,47,48,92]}; 2) defining an objective stretch indicator by basic flow indicators such as velocity and pressure ^[41,42].

Second-order statistics theory describes the PS mixing mechanism on the general flow structures. This theory focuses on the stretching term of the scalar dissipation rate equation ^[49]. The angle between the scalar gradient and the principal strain axis determines this term. By discussing the alignment dynamics ^[65,66,67], it is observed that the stretching mechanism, i.e. the mixing mechanism is different on the flow with different flow characteristics, as defined by velocity gradient tensors ^{[83,88,89,90,91]}. At present, one of the frontiers of second-order statistics theory research is to explain the relationship between "flow structure" and mixing thoroughly by scalar dissipation rate.

3.   Variable-density mixing mechanisms

Theories of single vortex PS mixing have guided the design of mixing enhancement devices such as lobe mixers ^[93] or strut mixers ^[94,95] that generate interacting streamwise vortex to enhance mixing ^[96]. However, as mentioned before, the flow density in the engine combustion chamber is always non-uniform. Therefore, it is necessary to understand the role of density gradient ^[97] and shock compression ^[98] in compressible mixing enhancement. The main focus of research on VD flows is on Richtmyer-Meshkov instability (RMI), Rayleigh-Taylor instability (RTI), VD jet/shear layer, and VD isotropic turbulence. As a classical VD flow, RMI/RTI has been reviewed before ^[99,100,101]. This section focuses on the characteristics of VD mixing compared to PS mixing. According to the mixing classification proposed by Dimotakis ^[23], VD mixing is the second-level mixing, which should be described using the component-transport equation:

where $Y_i$ is the mass fraction of $ith$ component.

Compared with PS mixing, VD mixing has two characteristics that significantly affect the mixing mechanism:

1) The misalignment between density gradient and pressure gradient can generate baroclinic vorticity $\frac{1}{\rho^2}(\nabla \rho \times \nabla p)_i = \frac{1}{\rho^2}(\varepsilon _{ijk}\frac{\partial \rho}{\partial x_j}\frac{\partial p}{\partial x_k})$. This vorticity can change the flow field ^[102]. As Dimotakis proved in Annu. Rev. Fluid Mech. ^[23]: "An important characteristic that distinguishes such flows from Level-1 mixing is the generation of baroclinic vorticity that derives from misalignments between pressure and density gradients, or, equivalently, temperature and entropy gradients in the flow."

2) The mixing governing function changes from the advection-diffusion equation to the component-transport equation in VD mixing. This change leads to two significant features that distinguish VD mixing from PS mixing.

The first feature is the velocity divergence $\frac{\partial u_j}{\partial x_j}$ is not 0 due to density. Livescu ^[103] provided a detailed review of the expressions of the velocity divergence in different VD flows. One of the expressions is given by Eq 3.3:

The second feature is that when the component-transport equation is written in the form of the advection-diffusion equation, a non-zero term appears on the right side of Eq 3.4, leading to the loss of hyperbolic persistence of the advection-diffusion equation ^[104].

where $\eta' = {\rho_1}/{\rho_2}$ is the density ratio. These two features brought by the change of mixing governing equation can lead to the difference of mixing indicators in VD mixing compared with PS mixing.

The VD mixing mechanism, similar to PS mixing, is influenced by fluid stretching. However, the misalignment between the density gradient and pressure gradient induces baroclinic vorticity, which causes alterations in the flow field, leading to changes in the evolution of mixing indicators. Consequently, research on VD mixing can be categorized into two parts shown in Figure 7: 1) The first part is dedicated to understanding the flow characteristics introduced by the density effect in VD mixing through linear stability theory and vortex dynamics. 2) The second part examines the evolution of VD mixing indicators, which includes first-order statistics focused on concentration, and second-order statistics focused on mixedness and scalar dissipation rate.

3.1.   Variable-density flow characteristics

VD mixing, for example RMI/RTI etc., classified as a second-level mixing process, is distinguished by its ability to alter the flow field. Richtmyer initially discovered the instability of shock interface interaction through experiments conducted in 1960 ^[109], with Meshkov later corroborating this phenomenon through numerical simulation ^[110]. Subsequently, this instability was named after these two researchers. The RMI phenomenon involves the amplification of a perturbed density discontinuity under shock waves causing the interface to develop unsteadily. Two key controlling parameters in RMI are the shock Mach number $Ma$, and the Atwood number $At = (\rho_2-\rho_1)/(\rho_2+\rho_1)$, which represents the density gradient at the interface.

The generation of baroclinic vorticity at the disrupted interface triggers the formation of complex structures, such as bubbles, spikes, and large-scale vortex structures. These structures then develop into more intricate turbulence. The RMI phenomenon occurs in different scenarios, including supernova explosion ^[111,112], inertial confinement fusion ^[113], detached-to-detonation transition ^[114], supersonic mixing ^[115,116] and supersonic combustion ^[117,118].

To study the flow characteristics of RMI, two primary methods are employed. The first method is to understand the disturbance growth law through linear stability theory ^[105]. The second approach is to describe the VD flow utilizing vortex dynamics and attempt to model and analyze the vortex ^{[106,119,120]}.

The earliest investigations of RMI were conducted utilizing linear stability theory. Richtmyer's assumption of a discontinuous shock wave revealed that the disturbance growth rate is influenced by the shock wave and initial disturbance ^[110]. Mikaelian went on to propose a nonlinear perturbation growth rate ^[105]. More recently, Luo et al. developed a nonlinear perturbation growth theory on column shock wave and interface interaction problems ^[121]. Similarly, the growth of disturbance in RTI garnered significant attention as with RMI. The Cook et al. study pointed out that the growth rate of RTI disturbance is related to $t^2$ ^[122,123]. Clark, in the same year, gave a comparable belief of $t^2$ for the RTI disturbance growth rate based on self-similarity theory ^[124].

In many RMI, VD jet, and shear layer problems, the VD flow characteristics are also described from the perspective of vortex dynamics. RMI, especially single-mode RMI, has a clear vortex structure, with several circulation prediction models proposed for this reason ^[125]. In comparison to constant density flow, Zabusky et al. discovered for the first time that secondary baroclinic vortex structure is a critical flow characteristic influencing the evolution of VD flows ^[102]. This structure significantly affects the middle and late-period flow evolution, and the phenomenon is also present in variable-density jet and shear layers. Marugan et al. modeled the impact of secondary baroclinic vortices on circulation and vortex structure by introducing the $Fr$ number ^[126]. Soteriou and Ghoniem further suggested that this structure would reduce the entrainment from heavy to light fluids ^[127]. This phenomenon is called asymmetric entrainment by Livescu, which means in RMI mixing layer, the light-density fluids tend to occupy a larger proportion than heavy-density fluids ^[103,128]. Much work has been done on the circulation and vorticity distribution in different types of RMI by Luo et al. ^{[129,130,131,132]}.

In PS mixing theories, what affects the mixing indicators are the local flow characteristics. In first-order statistics theory, scalar concentration is the primary focus, and the stretching of the lamellas influences mixing. Second-order statistics theory, which emphasizes scalar dissipation rate, is affected by the velocity gradient tensors represented by vorticity and strain. Whilst current linear stability theory and vortex dynamics theory can describe global VD flow characteristics, it remains challenging to link these flow characteristics and mixing indicators. As such, it is crucial to describe VD flow characteristics from the perspective of local characteristics.

3.2.   First-order statistics: dimensionless mixing time

In the framework of first-order statistics mixing theory, predicting mixing time using the stretching model is the primary focus of VD flow studies. This type of research primarily concentrates on shock bubble interaction (SBI). Though SBI's configuration is straightforward, it can be seen as a single vortical mixing once the shock interaction stage concludes. Therefore, SBI mixing is a distinct and particular branch of VD mixing problems.

When not considering the diffusion process, the bubble shape determines the mixing time of SBI. In general, the bubble shape transition can be observed by the stretching rate of the bubble boundary contour, as provided by Eq 2.22 ^[115,133], to obtain the mixing time. The stretchomg rate is a global stretching indicator comparable to the first-order statistics theory in PS mixing. However, this method is only effective when the boundary of the mass fraction contour can be identified, ensuring the corresponding stretching rate obtained is suitable to model the mixing time ^[107]. It is worth noting that other studies in RMI also focus on mixing time. Vorobieff et al. ^[134] found changes in the mixing characteristics during shock-gas film interaction at dimensionless time $u_1't/\lambda \approx 6.6$ ($u_1't$ is the gas velocity after the shock wave, $\lambda$ is the perturbation wavelength of the gas film). Niederhaus et al. ^[135] observed the mixing process of light and heavy gases under various shock Mach numbers and discovered that dissimilar working conditions can be better normalized under dimensionless time $u_1't/R$. Recently, based on Richtmyer's model ^[109] and the numerical simulation outcomes of Lombardini, Pullin, and Meiron ^[136], Oggian et al. ^[137] observed a self-similar growth in multi-mode RMI at the characteristic time $At^+\Delta ut/\lambda_{min}\approx 250$. Table 1 summarizes the mixing time in the current RMI/SBI studies.

The mixing time can be modeled based on the first-order statistics theory, as in PS mixing. SBI is a single-vortex mixing problem, and the stretching indicator can be obtained using the single-vortex model, after which the mixing time can be defined. Marble ^[30] and Meunier ^[34] developed the mixing time model for the single-vortex system in PS mixing given by Eq 2.11. Though the mixing characteristics, including stretching and diffusion effects, are considered, normalizing the mixing time of SBI is challenging because of the unconsidered secondary baroclinic vorticity characteristics of VD flow. In contrast to the mixing time model shown in Eq 2.11 proposed by Meunier et al. ^[19], Liu et al. ^[138] proposed a mixing time model predicting the mixing time accurately in VD SBI by considering additional stretching from second baroclinic vorticity:

where $\varGamma$ is the circulation of the vortex, $\mathscr{D}$ is the diffusion coefficient, $R_0$ and $s_0$ are the shape factors, $\eta_0$ is the compression rate ^[139], and $At^+$ is the post shock Atwood number. Compared with Marble's mixing time model ^[30], Liu et al.'s model proposed in Eq 3.5 accurately predicts mixing time in a large $Ma$ interval. Moreover, the ratio between the VD mixing time in Eq 3.5 and PS mixing time in Eq 2.11 can be defined as the mixing enhancement number:

3.3.   Second-order statistics: new variable-density mixing indicators

When investigating VD mixing using second-order statistics, the governing equation of mixing changes from the advection-diffusion equation (Eq 2.1) to the component-transport equation (Eq 3.1). As a result, the form of the mixing indicators also changes.

The first kind of VD mixing indicators redefines the concept of mixedness. Zhou ^[108] introduced an indicator similar to mixedness used in PS mixing when studying VD mixing of RMI and RTI and named it mixed mass. The physical meaning of mixed mass is the mass of the mixture produced by the mixing process, and it is defined as:

Mixed mass can also be related to an indicator similar to the scalar dissipation rate used in PS mixing:

This equation shares similarities with the mixedness equation for PS mixing 2.15. Therefore, it provides an analogy that helps explain the relationship between mixedness and scalar dissipation rate for VD mixing.

To normalized the mixed mass indicator, Zhou ^[108] used the following equation:

The denominator represents the mixture mass generated when all the fluids in the area are mixed. Together with two other mixing indicators, mixing length and molecular mixing fraction, these three indicators were used to study the asymptotic phenomenon in RMI and RTI. The mixing length indicator is defined as:

And, the molecular mixing fraction is given by:

These three indicators were found to have limits in RMI and RTI and tended to be around 0.8 shown in Figure 8. This asymptotic phenomenon was also observed in the numerical simulations by Thornber et al. ^[140].

Tian et al. also used mixed mole fraction and mixed mass to study RMI and RTI mixing, and observed similar asymptotic phenomenon ^[141,142]. They established a prediction model for the asymptotic limit ^[143], which is given by:

The second type of VD mixing indicators retains the original mixedness definition but proposed a new definition of scalar dissipation rate. Yu et al. ^[144] studied this problem in SBI and proposed a new definition of scalar dissipation rate for VD mixing based on the mixedness definition proposed by Cetegen ^[1].

The new scalar dissipation rate for mixedness in VD mixing is denoted as $\chi^*$ and can be expressed as

along with the coefficient of density gradient accelerated dissipation term $\mathcal{K} _{1, f} (\sigma, Y)$ and the redistributed diffusion term $\mathcal{K} _{2, f} (\sigma, Y)$. It was found that this new scalar dissipation rate has a weak correlation with $Pe$ and $Re$ numbers in their research. Using this new defination of scalar dissipation rate, the optimal Atwood number interval for different SBI can be obtained as in Figure 9 ^[145].

Furthermore, as the scalar dissipation rate characterizes the mixing rate, the original definition of scalar dissipation rate provided by Eq 2.14 was employed to express the mixing rate in RMI/SBI experiments (Figure 10).

The stretching mechanisms involved in VD mixing have been investigated by several researchers. Tian et al. ^[146] explored VD mixing in shock-turbulence interaction using the Favre-averaged scalar variance $\widetilde{{\phi''}^2}$, another second-order statistic. Wong et al. ^[147] described the evolution of mixedness

and scalar variance $\overline{{X'^2_{SF_6}}}$ in RMI with reshock (Figure 11). Livescu et al. ^[148] studied the vortex stretch in shock-turbulence interaction. Aslangil et al. ^[149] investigated the joint distribution and evolution of the second and third invariants of velocity gradients, denoted as $Q$ and $R$, respectively, which is closely related to stretching mechanism, in VD buoyancy-driven turbulence. Tian et al. ^[150] also examined the joint probability density function (PDF) $\mathcal{P} (Q, R)$ in the post-shock turbulence, and discovered that joint PDF is almost completely symmetrical in the heavy-fluid regions, while the joint PDF in light-fluid regions has a similar distribution to that of isotropic turbulence.

3.4.   A brief summary on variable-density mixing

Figure 7 summarizes the important areas for VD mixing research. In summary, compared to the well-established understanding of PS mixing, current understanding of VD mixing is still at an early stage. VD mixing has two main characteristics: VD flow characteristics can be influenced by baroclinic vorticity (Eq 3.2), and the mixing indicators change due to the change of the mixing governing equation. VD flow characteristics have been described through linear stability theory ^{[105,109,110,121,122,123]} and vortex dynamics theory ^{[102,106,125,126,129,130,131,132]}. However, VD mixing indicators are affected by local flow indicators such as the velocity gradient tensors represented by vorticity and strain ^[49], which brought the challenge to relate mixing indicators and the flow characteristics described by traditional theory. In terms of mixing indicators, research on first-order statistics has mainly focused on the dimensionless mixing time in SBI and RMI problems ^{[134,135,137,138,153]}. Furthermore, research has been conducted on second-order statistics to define new VD mixing indicators and study VD mixing phenomena by these new VD indicators ^{[108,141,144]}. The evolution of second-order statistics in VD mixing problem have also been preliminarily discussed ^{[138,146,147]}. Therefore, coupling the velocity gradient tensor equations and the scalar dissipation rate equation, describing the flow characteristics of VD mixing from the perspective of velocity gradient tensor, and exploring the relationship between scalar dissipation rate and VD flow characteristics, are essential future works for VD mixing research.

4.   Conclusions and future work

Mixing is a fundamental problem in fluid mechanics that has been studied from many different perspectives. PS mixing has been studied in terms of chaotic mixing by Ottino ^[154], turbulent mixing by Warhaft ^[16], Dimotakis ^[23] and Screenivasan ^[155], and mixing versus stirring by Villermaux ^[31]. Meanwhile, VD mixing has been analyzed from the perspective of stability theory by Brouillette ^[99], vortex dynamics by Zabusky ^[106] and Ranjan ^[156], VD turbulence by Livescu ^[103] and Zhou ^[100,101], and turbulent stratified flow by Caulfield ^[10]. Despite much progress already made, the relationship between flow and mixing remains the central problem of mixing research and requires further investigation. Understanding this relationship offers the potential for more effective mixing enhancement in applications such as combustor design. Two types of mixing indicators, first-order statistics and second-order statistics, can reflect this relationship. This paper systematically reviews these mixing indicators and their dynamics in PS and VD mixing, as illustrated in Figure 12. Based on the relevant studies reviewed in this paper, the following conclusions can be drawn:

Regarding PS mixing:

1) The first-order statistics theory focuses on the concentration $c$ and reveals the mixing mechanisms on the lamella structures, as highlighted in Figure 12. The advection-diffusion equation (Eq. 2.1) governing the concentration on the moving lamellas can be simplified to the Fourier equation (Eq. 2.4) through the Ranz transform for analytical solution. Within this theoretical framework, key mixing indicators such as the probability distribution of concentration $p(c)$, mixing time $t_s$, and the Batchelor scale $\eta$ can be defined.

2) The second-order statistics theory focuses on the mixedness $f$ and scalar dissipation rate $\chi$ and reveals the mixing mechanisms on the general flow structures. An in-depth understanding of flow and mixing through scalar dissipation rate can be established by examining the stretching term $\nabla c \cdot \boldsymbol{S} \cdot \nabla c$ in the scalar dissipation rate equation (Eq 2.16). As highlighted by Figure 12, this stretching term is primarily determined by the alignment dynamics of scalar gradient, and by investigating these alignment dynamics, it is shown that scalar dissipation rate exhibits different dynamical properties in flows with different characteristics, which are determined by velocity gradient tensors. In two-dimensional flow, the stretching of the scalar dissipation rate is governed by dimensionless numbers $r$ (Eq. 2.32) and $n$ (Eq. 2.33), while in three-dimensional flow, it is determined by the invariants of the velocity gradient tensor (Eq. 2.34) and the rotation of the strain axis (Eq. 2.27).

Regarding VD mixing:

1) In comparison to PS mixing, baroclinic vorticity generated by the misalignment of the density gradient and pressure gradient can alter the flow in VD mixing. This VD flow characteristic is described from two perspectives: linear stability theory and vortex dynamics. These two methods describe the global characteristic of VD flow, which can hardly relate to molecular mixing.

2) Research on VD mixing under the framework of first-order statistics theory primarily focuses on the dimensionless mixing time in SBI and RMI. As highlighted in Figure 12, mixing time is primarily determined by the stretching rate of the bubble boundary contour given by Eq 2.22. During the modeling of mixing time, new phenomena were discovered, including the second baroclinic vorticity leading to additional stretching and accelerating the mixing process of SBI as in Eq 3.5.

3) Owing to change of the mixing governing equation, second-order indicators for VD mixing are re-defined, such as the mixed mass $\mathcal{M}$ in Eq 3.7 and new scalar dissipation rate $\chi^*$ in Eq 3.14. These mixing indicators have been applied to study VD mixing phenomena. Furthermore, the quantities related to the stretching mechanism, such as velocity gradient tensors invariants, have also been primarily investigated in VD mixing problems.

Based on these conclusions, future research in mixing has the following trends:

Regarding PS mixing:

1) Scalar dissipation rate will be key in mixing research and needs to be explored in some extreme conditions, such as turbulence with high $Re$ and Schmidt numbers $Sc$.

2) The relation between the scalar dissipation rate and "flow structures" needs further clarity.

3) An engineering method for designing flow to achieve optimal mixing is necessary.

Regarding VD mixing:

1) The flow characteristics caused by baroclinic vorticity should be described using local quantities such as velocity gradient tensors.

2) The new second-order statistics should also be explored in limiting conditions in VD turbulence with extremely high $Re$ and $Sc$ numbers.

3) The relationship between VD flow characteristics and VD mixing indicators deserves further study.

This review mainly focuses on non-reacting flow. As mentioned above, scalar dissipation rate is also crucial for studying combustion phenomena and should be used to improve our understanding of combustion and the design of combustion related devices.

Acknowledgement

This work was supported by the Natural Science Foundation of China (nos. 91941301).

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interests

There is no conflict of interest.