Effect of Schmidt Number on Forced Isotropic Turbulence with Passive Scalars

Abstract

Traditionally, Fourier spectra have been employed to gain a deeper understanding of turbulence flow structures. The investigation of isotropic forced turbulence with passive scalars offers a straightforward means to examine the disparities between velocity and passive scalar spectra. This flow configuration has been extensively studied in the past, encompassing a range of Reynolds and Schmidt numbers. In this present study, direct numerical simulations (DNS) of this flow are conducted at sufficiently high Reynolds numbers, enabling the formation of a wide inertial range. The primary focus of this investigation is to quantitatively assess the variations in scalar spectra with the Schmidt number (Sc). The spectra exhibit a transition from a k^−5/3 scaling for low Sc to a k^−4/3 scaling for high Sc. The emergence of the latter power law becomes evident at Sc = 2, with its width expanding as Sc increases. To gain further insights into the underlying flow structures, a statistical analysis is performed by evaluating quantities aligned with the principal axes of the strain field. The study reveals that enstrophy is primarily influenced by the vorticity aligned with the intermediate principal strain axis, while the scalar gradient variance is predominantly controlled by the compressive strain. To provide a clearer understanding of the differences between enstrophy and scalar gradient variance, joint probability density functions (PDFs) and visualizations of the budget terms for both quantities are presented. These visualizations serve to elucidate the distinctions between the two and offer insights into their respective behaviors.

1. Introduction

The Schmidt number is a dimensionless quantity that characterizes the relative importance of momentum diffusion to the diffusivity of passive scalars, namely $\mathit{Sc}=\nu /{\nu }_{\theta }$. It is particularly relevant when dealing with the transport of passive scalars, such as temperature, concentration, or other non-momentum-transporting properties, in a turbulent flow. In forced isotropic turbulence, external energy is continuously injected into the system to maintain turbulence. Passive scalars in this context refer to quantities that do not actively contribute to the turbulent energy generation but are instead transported by the turbulent flow. The effect of the Schmidt number on forced isotropic turbulence with passive scalars strongly impacts several physical processes. At low Schmidt number ($\mathit{Sc}<<1$) molecular diffusivity is high compared to the turbulent mixing. In this case, passive scalars are mixed efficiently, and their gradients are rapidly flattened due to the dominant influence of molecular diffusion. As a result, sharp gradients in the scalar field are quickly smoothed out. At higher Schmidt numbers ($\mathit{Sc}>1$), molecular diffusion becomes less significant compared to turbulent mixing. This allows for concentration fluctuations to persist over longer-length scales before being smoothed out by diffusion. Consequently, the scalar field might exhibit more complex and smaller-scale structures. The scalar dissipation rate, which quantifies how quickly the scalar variance is dissipated, is affected by the Schmidt number, as lower Schmidt numbers generally result in higher scalar dissipation rates due to the dominance of diffusion, whereas higher Schmidt numbers lead to lower scalar dissipation rates due to enhanced mixing.

Direct numerical simulations (DNS) have long been employed as a powerful tool for investigating the intricate dynamics of passive scalars advected by turbulent flows. In the context of studying these phenomena, isotropic turbulence has often been simulated using DNS with a low wavenumber forcing. This setup enables the attainment of statistically stationary states within a few eddy turnover times, facilitating the collection of extensive datasets for an in-depth analysis of high-order statistical properties. It is worth noting that the extensive body of literature on forced isotropic turbulence with and without passive scalars makes it impractical to provide a comprehensive list of relevant papers in this brief discussion. However, the review by Ishihara et al. [1] provides a comprehensive compilation of papers on forced isotropic turbulence without passive scalars, while Donzis et al. [2] catalog the studies specifically focused on passive scalars, shedding light on the differences and similarities between velocity and passive scalar statistics, including their spectral properties. An alternative approach to studying passive scalars in turbulent flows involves DNS of decaying turbulence, initiated with prescribed spectra and random phases. Although this method deviates from realistic initial conditions, it provides insights into the effects of varying parameters, such as the Schmidt number. Orlandi and Antonia [3] conducted a study specifically focused on understanding the influence of Schmidt number variations on the decay of passive scalars. They compared DNS results with those obtained using eddy-damped quasi-normal Markovian (EDQNM) closures, which offer information about high Reynolds number regimes that would otherwise require computationally prohibitive efforts to achieve with DNS. By leveraging DNS and incorporating spectral closures, researchers have been able to delve into the complex dynamics of passive scalars in turbulent flows, exploring a wide range of scenarios and parameters. The combination of these techniques allows for a comprehensive investigation of the complex physics underlying the behavior of passive scalars, providing insights that contribute to our understanding of turbulent transport processes.

In a distinct approach, Orlandi et al. [4] employed a minimal flow unit (MFU) framework to investigate various aspects of isotropic turbulence and the turbulent transport of passive scalars. They initiated their simulations by introducing two orthogonal Lamb dipoles, allowing them to interact and evolve. Initially, the flow exhibited a vortical-dominated stage where vortex structures of smaller scales were generated through vortex stretching. During the initial phase, the number and shape of the vortex structures generated from the initial conditions were similar to those generated under inviscid conditions. However, as the simulation progressed and approached the time of the finite-time singularity, the vortical structures in the viscous case became dependent on the viscosity. The statistical properties, including the velocity and passive scalar spectra, coincided with those observed in forced isotropic turbulence.

The utilization of the MFU framework provided a detailed analysis of the dynamics of the components of the vorticity and scalar gradient vectors along with the principal axes of the strain rate tensor. These analyses corroborated the findings reported by Tsinober [5] regarding the importance of the intermediate strain in the evolution of enstrophy and the compressive strain in the scalar gradient variance. Additionally, flow visualizations were presented to enhance the understanding of the underlying physics. Notably, Orlandi et al. [4] limited their simulations to a passive scalar with unit Schmidt number.

Through the implementation of the MFU framework, Orlandi et al. [4] provided a comprehensive examination of the complex dynamics and statistical properties of isotropic turbulence and passive scalar transport. Their work shed light on the behavior of vortical structures and their relationship with viscosity, as well as the crucial roles played by different strain components in the evolution of enstrophy and scalar gradient variance.

In the present paper, simulations of forced turbulence were conducted at both high and low Reynolds numbers to investigate the transport of passive scalars, at various Schmidt numbers. Two different forcing methods were employed for the passive scalar, one similar to that used for the velocity field to maintain a constant passive scalar variance during the evolution towards the statistical steady state, and the other method employed by Donzis et al. [2]. The choice of forcing method had an impact on the dynamics, leading to oscillations in the scalar dissipation. The simulations using the second forcing method exhibited growing oscillations, rendering them unsuitable for analyzing the spectra and the terms in the balance equations of enstrophy and scalar gradient variance.

One unresolved question is whether the scalar spectrum, at any Schmidt number, exhibits an inertial range with a power-law decay of $k^{-5/3}$ similar to the kinetic energy spectrum, or maybe other decay laws. In this respect, Lohse [6] argued that in sheared turbulence, passive scalars should have a $k^{-4/3}$ spectral range. In fact, based on the assumption that in shear-dominated turbulence, the velocity spectrum should display a $E_u\sim k^{-7/3}$ spectral range [7], the respective Fourier coefficients should scale as $\widehat{u}\left(k\right)\sim k^{-2/3}$. Assuming the constancy of the passive scalar flux in wavenumber space, $T_{\theta }\left(k\right)\sim k\widehat{u}\left(k\right){\widehat{\theta }}^2\left(k\right)\sim \chi$ (with $\chi$ the scalar dissipation rate), it readily follows that $E_{\theta }\left(k\right)\sim k^{-4/3}$. Support for this scaling law came from measurements of temperature fluctuations in the wake of a heated cylinder [8]. A theory for the behavior of scalar spectra in stratified turbulence was proposed by Bolgiano [9] and Obukhov [10]. In analogy with Kolmogorov’s theory, and assuming that next to a given flow scale k, the only relevant parameters are the thermal dissipation rate and the product of the thermal expansion coefficient and gravity, a dimensional analysis readily yield the prediction $E_{\theta }\left(k\right)\sim k^{-4/7}$, which is commonly referred to as the Bolgiano–Obukhov scaling. Numerical experiments by Camussi and Verzicco [11] in cylindrical cells supported the reduced decay predicted by the Bolgiano–Obukhov scaling, as shown in Figure 5 of Lohse and Xia [12]. Orlandi and Pirozzoli [13] carried out DNS of natural convection at high Rayleigh numbers and a Prandtl number $\mathit{Pr}=1$ and observed a wide power-law spectral range with an exponent close to either $k^{-4/3}$ or $k^{-7/5}$.

To explore the dependence of the decay exponent on the Schmidt number and shed light on the differences in the dynamics of enstrophy and scalar gradient variance, DNS of forced isotropic turbulence transporting passive scalars at various $\mathit{Sc}$ was conducted in this study. These simulations provide information about the Schmidt number dependence of the decay exponent and contribute to understanding the distinct behaviors of enstrophy and scalar gradient variance.

2. Relevant Equations

The momentum and continuity equations for incompressible flows are as follows

$$\frac{\partial u_i}{\partial t}+\frac{\partial u_i{u}_j}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\frac{1}{\mathit{Re}}\frac{{\partial }^2{u}_i}{\partial x_j\partial x_j}+F_i;\frac{\partial u_j}{\partial x_j}=0,$$

(1)

where $u_i$ are the components of the velocity vector in the i directions, p is the pressure, $x_1$, $x_2$, and $x_3$ are the three orthogonal space directions, and $F_i$ is the external forcing to prevent decay of the turbulence kinetic energy $q=<u_i^2>/2$. The passive scalar $\mathit{T}$ is transported by the velocity field according to

$$\frac{\partial \mathit{T}}{\partial t}+\frac{\partial {\mathit{Tu}}_{\mathit{j}}}{\partial x_j}=\frac{1}{\mathit{Re}\mathit{Sc}}\frac{{\partial }^2\mathit{T}}{\partial x_j\partial x_j}+F_{\theta },$$

(2)

where $F_{\theta }$ is introduced to prevent a decay of the passive scalar variance, $\mathsf{\Theta }=<{\mathit{T}}^2>/2$. The numerical simulations are initiated with velocity and scalar fields with random phases, representing flows without structures but with high energy at low wavenumbers and small energy at large wavenumbers. The initial kinetic energy spectrum is given by the equation $E\left(k\right)=64{(k/8)}^4{e}^{-(k/8)}/8$, with $q=\int E\left(k\right)=1.5$, and $\int k^{2}E\left(k\right)=180$. The reference velocity is $u_{\mathrm{rms}}=2q/3$, and the Taylor Reynolds number (later defined) is related to the initial spectrum $E\left(k\right)$ through

$$\mathit{Re}={\mathit{Re}}_0{\left(\frac{3}{20}\right)}^{1/2}\frac{{\left(2\int k^{2}E\left(k\right)\right)}^{1/2}}{\int E\left(k\right)}.$$

The value of $\mathit{Re}$ is determined by assigning a value to ${\mathit{Re}}_0$, as well as specifying the values of $\mathit{Sc}$. The same spectrum is used to generate the initial distribution of the passive scalar field $\mathsf{\Theta }$. The velocity components and the passive scalar in physical space are obtained using the method described by Rogallo [14].

The velocity components and the passive scalar in wavenumber space are given by

$$\begin{array}{c}\hfill {\widehat{u}}_1=\left(\alpha \right|k|k_2+\beta k_1{k}_3)/(\left|k\right|k_h)\\ \hfill {\widehat{u}}_2=(-\alpha |k|k_2+\beta k_3{k}_2)/(\left|k\right|k_h)\\ \hfill {\widehat{u}}_3=(-\beta k_h^2)/\left(\right|k|k_h)\\ \hfill \widehat{\mathit{T}}=(-\gamma k_h^2)/\left(\right|k|k_h),\end{array}$$

(3)

where the hat symbol is used to denote complex quantities. This expression satisfies the condition for a divergence-free velocity field in wavenumber space. Additionally, it can be demonstrated that if the following expressions are satisfied

$$\alpha =\sqrt{\frac{E\left(k\right)}{2\pi k^2}}{e}^{\iota {\theta }_1}cos\left(\varphi \right),\beta =\sqrt{\frac{E\left(k\right)}{2\pi k^2}}{e}^{\iota {\theta }_2}sin\left(\varphi \right),\gamma =\sqrt{\frac{E\left(k\right)}{2\pi k^2}}{e}^{\iota {\theta }_3},$$

(4)

the three-dimensional energy spectrum of the velocity field matches the given $E\left(k\right)$ spectrum. To distribute the angles ${\theta }_1,{\theta }_2,{\theta }_3$, and $\varphi$ uniformly over the interval $(0,2\pi )$, four sets of random numbers are used. By applying the fast Fourier transform (FFT), the initial velocity and thermal fields in physical space are obtained.

To achieve a statistically steady state in the transport equations, forcing terms are introduced. Various forcing methods were discussed by Eswaran and Pope [15], based on the concept that at high Reynolds numbers, the dynamics of the small scales should be decoupled from the large scales. However, manipulating the large scales can impact the rate of energy dissipation. A statistical steady state is reached when the dissipation rate oscillates around a constant value over time. Without forcing, the root mean square velocity $<u_i^2>=\int E_i(k,t)dk$, at each instant, has a total energy component $q_i\left(t\right)=\int E_i(k,t)dk$, that is less than the initial energy $q_0=\int E_i(k,0)dk$. By selecting a threshold wavenumber ${\left|k\right|}_F$, the energy $\Delta q_i\left(t\right)$ within the range $\left|k\right|\le {\left|k\right|}_F$ is evaluated to determine the quantity $F_i=(q_0-q_i\left(t\right))/(\Delta q_i\left(t\right))$, which represents the proportion of energy lost. By modifying the ${\widehat{u}}_i\left(n\right)$ components for $\left|k\right|\le {\left|k\right|}_F$ through ${\widehat{u}}_i=\sqrt{F_i}{\widehat{u}}_i\left(n\right)$, the total energy $q_i$ can be kept constant over time. In the simulations described in this paper, a numerical scheme was employed in physical space, and thus three-dimensional FFT operations were used to transition between physical and wavenumber spaces. Within the wavenumber space, the ${\widehat{u}}_i$ components were modified for $\left|k\right|\le {\left|k\right|}_F$. Subsequently, an inverse three-dimensional FFT was applied to obtain the modified velocity in physical space. To save computational time, the forcing was only applied during the last of the three Runge–Kutta steps used for time advancement. The same forcing approach was also employed in the solution of the passive scalar transport equation, resulting in a constant $\mathsf{\Theta }\left(t\right)=\int E_{\theta }(k,t)dk$, throughout the time evolution.

In contrast to the forcing scheme just described, Donzis et al. [2] used a forcing term for Equation (2), of the type $F_{\theta }=-u_{i}dT/d{x}_i$, mimicking the effect of a mean passive scalar gradient as found in natural convection flow. They assumed a constant gradient $dT/d{x}_i=1$. In this scenario, not only does the rate of scalar dissipation vary over time, but the quantity $\mathsf{\Theta }\left(t\right)$ also experiences temporal fluctuations. Both forcing methods, as described by Eswaran and Pope [15] and Donzis et al. [2], were utilized in the present paper, and their respective results are discussed later on.

The system of equations given by (1) and (2), subject to periodic boundary conditions in the Cartesian directions $x_i$, was numerically solved using a finite difference scheme as described in Orlandi [16]. To handle the temporal integration, an implicit Crank–Nicholson scheme was employed for the viscous terms, while a third-order Runge–Kutta scheme was utilized for the convective terms. By adopting a staggered arrangement of the velocity components, any odd–even decoupling phenomena were effectively eliminated, ensuring the discrete conservation of the total kinetic energy in the inviscid limit. The nondimensional temperature $\mathit{T}$ was placed at the same location as $u_2$, which proves particularly advantageous in the presence of mean stratification, as it preserves the sum of potential and kinetic energy. To facilitate the implementation of MPI (message passing interface) coding, the computational domain was divided into layers parallel to the $x_2$ direction.

3. Results

3.1. Global Flow Parameters and Validation

Numerical simulations were performed at several Reynolds and Schmidt numbers, with an appropriate resolution chosen to ensure the maximum resolved wavenumber $k_M^{\ast }\approx 1$, where the superscript * indicates a normalization in Kolmogorov units. The following definitions are required to obtain this normalization from nondimensional computational units

$$\begin{array}{c}\hfill u^{\prime }=\sqrt{\int Edk/3},D_v=\int E{k}^{2}dk,ϵ=2D_v/\mathit{Re},\eta ={\left(ϵ{\mathit{Re}}^3\right)}^{-1/4}\\ \hfill \lambda =\sqrt{\frac{10\int Edk}{ϵ\mathit{Re}}},R_{\lambda }=\mathit{Re}\lambda u^{\prime },S_E={\left(ϵ{\mathit{Re}}^5\right)}^{-1/4}{S}_U=S_E/\eta \\ \hfill T^{\prime }=\sqrt{\int E_{\theta }dk},D_{\theta }=\int E_{\theta }{k}^{2}dk,\chi =\frac{2D_{\theta }}{\mathit{Re}{\mathit{Sc}}^{3/2}},{\eta }_b=\eta /\sqrt{\mathit{Sc}},S_T=\frac{\chi {\eta }^{5/3}}{{ϵ}^{1/3}},S_R=\frac{S_T}{{\eta }_b}.\end{array}$$

(5)

The energy spectra in Kolmogorov units were obtained by scaling the velocity energy spectrum by $S_E$, which represents the normalization factor for energy. Similarly, the root-mean-square (rms) velocity was obtained by dividing the square root of the velocity variances in computational units by $S_U$, the normalization factor for velocity. The scalar energy spectrum was scaled by $S_T$, while the passive scalar rms was normalized by $S_R$.

The global data for the relevant simulations regarding the velocity are presented in

Table 1. Global quantities related to the solution of Equation (1).

Flow Case	${\mathit{N}}_{\mathit{i}}$	$\mathit{Re}$	${\mathit{R}}_{\mathit{\lambda }}$	${\mathit{K}}_{\mathit{M}}^{\ast }$	${\mathit{S}}_{\mathit{E}}{10}^3$	${\mathit{S}}_{\mathit{U}}10$	$\mathit{\lambda }/\mathit{\eta }$	${\mathit{D}}_{\mathit{v}}$	$\mathit{ϵ}$
A	1152	2694	438	1.89	4.19	1.27	34	581	0.44
B	1152	2021	374	2.36	5.97	1.47	32	433	0.43
C	768	367	164	5.72	49.7	3.33	21	74	0.41
D	768	367	164	5.72	49.7	3.33	21	74	0.41
E	768	2449	442	1.41	4.53	1.23	35	457	0.37

, and the data for the passive scalar are reported in

Table 2. Global quantities related to the solution of Equation (2), including the Schmidt number ($\mathit{Sc}$), Batchelor length scale (${\eta }_B$), and normalization factors for the scalar energy spectrum ($S_T$), and for the passive scalar rms ($S_R$), passive scalar variance ($\mathsf{\Theta }$), the integrated scalar spectrum dissipation ($D_{\theta }$), and the scalar gradient dissipation rate ($\chi$).

Flow Case	$\mathit{Sc}$	${\mathit{\eta }}_{\mathit{b}}$$(\times {10}^2)$	${\mathit{K}}_{\mathit{b}}^{\ast }$	${\mathit{S}}_{\mathit{T}}$$(\times {10}^4)$	${\mathit{S}}_{\mathit{R}}$$(\times 10)$	$\mathbf{\Theta }$	${\mathit{D}}_{\mathit{\theta }}$	$\mathit{\chi }$
$A1$	0.5	4.645	2.676	1.624	5.878	3.0	1145	1.70
$A2$	0.75	3.840	2.212	1.686	5.448	3.0	1719	1.70
$A3$	1.0	3.285	1.892	1.615	4.917	3.0	2278	1.69
$A4$	1.5	2.716	1.564	1.450	4.359	3.0	2955	1.46
$A5$	2.0	2.352	1.355	1.741	5.234	3.0	4732	1.75
$B1$	2.5	2.593	1.494	2.283	5.567	3.0	4139	1.64
$B2$	3.0	2.367	1.364	2.273	5.544	3.0	4947	1.63
$B3$	4.0	2.050	1.181	2.254	5.498	3.0	6541	1.62
$C1$	2.0	10.55	4.053	15.75	10.56	3.0	474	1.29
$C2$	4.0	7.463	2.866	15.65	10.49	3.0	941	1.28
$C3$	6.0	6.093	2.340	15.24	10.21	3.0	1375	1.24
$C4$	8.0	5.272	2.026	14.61	9.790	3.0	1758	1.20
$C5$	10.0	4.720	1.812	14.86	9.958	3.0	2235	1.21
$C6$	16.0	3.731	1.433	13.88	9.302	3.0	3341	1.14
$D1$	2.0	10.54	4.046	24.76	16.62	3.47	749	2.04
$D2$	4.0	7.450	2.861	25.71	17.26	3.82	1555	2.11
$D3$	6.0	6.083	2.336	26.31	17.66	4.04	2386	2.16
$D4$	8.0	5.268	2.023	26.76	17.96	4.19	3236	2.20
$D5$	10.0	4.712	1.809	27.12	18.20	4.32	4100	2.23
$D6$	16.0	3.725	1.430	27.87	18.70	4.58	6740	2.29
$E1$	1.0	3.675	1.411	1.992	5.420	3.00	2006	1.64

. The results presented in Table 1 demonstrate that a series of simulations with an $N_i=1152$ grid resolution achieve a sufficiently high Taylor Reynolds number ($R_{\lambda }$) to capture a power-law inertial range characterized by a $k^{-5/3}$ scaling. Moreover, the maximum resolved wavenumber $k_M^{\ast }>1$ indicates that the simulations effectively resolve the exponential decay range of the spectrum. Simulations with an $N_i=768$ resolution were also conducted, employing the same random forcing as the $N_i=1152$ cases. The influence of the passive scalar forcing was investigated using two different forcing methods as previously mentioned. The simulations at $R_{\lambda }=164$ aimed to study the effect of high Schmidt numbers ($\mathit{Sc}>4$) on the passive scalar spectra. In Table 2, the global quantities related to the passive scalar are denoted by the same letter as in Table 1, with an additional index corresponding to the specific value of the Schmidt number.

The simulations conducted in this study exhibit a linear growth of the dissipation rate $D_v$ with respect to the Reynolds number ($\mathit{Re}$), approximately given by $D_v\approx 0.225\mathit{Re}$. Consequently, the energy dissipation rate $ϵ$ remains nearly independent of the Reynolds number. A similar linear growth of $D_v$ with $\mathit{Re}$ was observed in the interaction of two orthogonal Lamb dipoles by Orlandi et al. [4]. In that case, without any external forcing, $ϵ$ increased over time at a rate dependent on $\mathit{Re}$ and eventually reached the same value at a certain nondimensional time. Subsequently, $ϵ$ decreased over time independently of the Reynolds number. In the present simulations, the forcing was employed to achieve a statistical steady state within a certain number of eddy turnover times. This transient phase involved the formation of turbulent structures that were absent at the initial time. Depending on the specific forcing method employed, the energy may remain constant during this transient phase, while significant variations in $D_v$ and $ϵ$ can occur. At a certain point in the simulations, the statistical quantities started to exhibit oscillations. From this point onward, time averages were computed to obtain the mean values of $D_v$ and $ϵ$, which are presented in Table 1. Additionally, the effect of the forcing method on the spectra was investigated, and comparisons were made with pseudospectral simulations available in the literature.

To validate the accuracy of the present second-order finite difference scheme and examine the spectra, comparisons were made with the results from the pseudospectral simulations conducted by Jiménez et al. [17]. The pseudospectral simulations considered various values of the Taylor Reynolds number, with the highest reported at $R_{\lambda }=140$. In the simulations by Jiménez et al. [17], the forcing at wavenumbers ${k<\left|k\right|}_F$ was implemented by assuming a negative viscosity coefficient of appropriate magnitude, adjusted periodically to maintain $k_M^{\ast }\approx 2.5$. In the simulations conducted by Donzis et al. [2] at $R_{\lambda }=650$, the forcing in the momentum equation was implemented by fixing the $E\left(k\right)$ spectrum for wavenumbers ${k<\left|k\right|}_F$, based on prior simulations that utilized a stochastic forcing scheme. With this type of forcing, the total energy varies over time, and the spectrum at low wavenumbers appears relatively smooth. In contrast, in the present simulations, the total energy was kept constant, resulting in oscillations in the spectra at low wavenumbers. Figure 1 compares the spectra of cases A, B, and C with those from the simulations of Jiménez et al. [17] and Donzis et al. [2]. The longitudinal and transverse one-dimensional spectra (Figure 1a and Figure 1b, respectively) highlight the effects of the forcing method. In wavenumber space, the simulations exhibit a smoother trend at low wavenumbers. In contrast, the simulations in physical space result in spectra that show an increase in energy content at low $k^{\ast }$, potentially leading to a narrowing of the inertial range.

The differences in the spectra at low and intermediate wavenumbers can indeed be observed in cases C and D, which have a large $k_M^{\ast }$ and sufficient resolution to simulate the evolution of the thermal field at high Schmidt numbers. The compensated spectra in Figure 1 highlight these differences and emphasize the amplitude of the spectral “bottleneck” phenomenon, namely the energy pile-up at the start of the dissipative range. It is worth noting that the one-dimensional spectra of Donzis et al. [2] in Figure 1 were derived from the three-dimensional spectra using the relationships reported by Pope [18] (p.216). Comparing the three-dimensional spectra in Figure 1c, it can be observed that they exhibit a better agreement. Both types of simulations, using the stochastic forcing scheme and the present simulations in physical space, produce a spectral bottleneck with a significant amplitude. These findings support the conclusions of Bogucki et al. [19], Dobler et al. [20], who argued that a strong bottleneck is only observed in numerical simulations, in which three-dimensional spectra can be directly evaluated. In contrast, experimental spectra are generally derived from frequency spectra using the Taylor hypothesis, which explains why a clear spectral bottleneck is not observed.

In their work, Su et al. [21] conducted a comprehensive analysis of three-dimensional spectra obtained from direct numerical simulations. In Figure 2 of their study, they presented compensated spectra, highlighting their collapse around a peak that varied between a value of 2.2 at $R_{\lambda }=140$ and 2.0 at $R_{\lambda }=1000$, located at $k^{\ast }=0.128$. The present data for cases A and B, shown in Figure 1d, exhibit a good agreement with the spectra reported by Donzis et al. [2] at $R_{\lambda }=650$, using the same scaling as in Su et al. [21]. On the other hand, the spectra of cases C and D align well with those presented by Jiménez et al. [17] at $R_{\lambda }=140$. This agreement further supports the consistency between the present simulations and the findings reported in the literature.

The global data presented in Table 2 provide insights into the variations of quantities related to the passive scalar with the Schmidt number. It is noteworthy that for all cases, the maximum value of $K_b^{\ast }$ indicates that even at the highest $\mathit{Sc}$, Equation (2) is fully resolved. Furthermore, $D_{\theta }$ increases with the Peclet number ($\mathit{Pe}$). The data in Table 2 suggest an approximate relationship $D_{\theta }\approx 0.8\mathit{Pe}$, implying that the scalar dissipation rate $\chi$ remains nearly constant. It is important to note that this relationship only holds for the simulations with a constant $\mathsf{\Theta }$, leading to the exclusion of flow case D, where $\mathsf{\Theta }$ varied. This observation highlights the influence of the forcing method on the computed results.

To validate the passive scalar results at $\mathit{Sc}=1$, a comparison was made with the pseudospectral simulation by Donzis et al. [2], denoted as case E in the present study. Although the forcing method in Donzis et al. [2] did not maintain a constant $\mathsf{\Theta }$, their higher resolution diminished the impact of the forcing method. The compensated three-dimensional spectra of the passive scalar and energy in Figure 2a demonstrate a good agreement between the present simulations and the results of Donzis et al. [2]. In our simulations, the relatively short inertial range was primarily due to reaching a lower Reynolds number ($R_{\lambda }=442$) with a resolution of $N_i=768$, which was smaller than the resolution used by Donzis et al. [2] to achieve $R_{\lambda }=650$. Both simulations, however, indicate a larger amplitude of the passive scalar spectral bottleneck compared to that of the velocity field. This difference can be attributed to the different geometry of vortical structures and scalar gradient structures, as commented later. Figure 2b illustrates that similar to the velocity field, the amplitude of the spectral bottleneck in the three-dimensional spectrum of the passive scalar is larger than that in the one-dimensional spectrum. This observation suggests that our low-wavenumber forcing method leads to a less pronounced increase in the amplitude of $E_{\theta }\left(k_i\right)$, particularly in the three-dimensional spectrum.

3.2. Effects of Schmidt Number on Spectra

Efforts have been made in the past to provide a theoretical justification for the presence of a $k^{-1}$ range in three-dimensional spectra at high Schmidt numbers, following the inertial range and preceding the exponential decay range. The work of Batchelor [22] played a significant role in this regard, arguing that scalar fluctuations existed only in regions where the spatial variations of the velocity field were linear. These fluctuations can be expressed as the product of the principal strain rates and the position in that eigenframe. Furthermore, Batchelor [22] assumed that the strain rates were locally uniform in space for wavenumbers much larger than the Kolmogorov wavenumber ($k>>1/\eta$) and persistent in time, neglecting temporal fluctuations. In that scenario, the fluctuating scalar gradients align preferentially with the eigenvector corresponding to the most compressive principal strain rate, denoted as ${\tilde{S}}_3$. This alignment leads to the emergence of a $k^{-1}$ range in the scalar spectra. The work of Donzis et al. [2] revisited and expanded upon the arguments put forth by Batchelor [22] to understand the behavior of passive scalar fluctuations in turbulent flows. They considered the effects of strain rates, spatial uniformity, and temporal persistence on scalar gradients, particularly in the context of high Schmidt numbers.

The theoretical considerations of Batchelor [22] led to an expression for the three-dimensional scalar spectrum, which takes into account the effects of strain rates and spatial variations in the velocity field,

$$\frac{E_{\theta }}{\chi {\tau }_{\eta }{\eta }_B}=C_B{\left(k{\eta }_B\right)}^{-1}exp(-C_B{\left(k{\eta }_B\right)}^2),{\tau }_{\eta }={(\nu /ϵ)}^{1/2},$$

(6)

where $\chi$ is the scalar dissipation rate, ${\tau }_{\eta }$ is the Kolmogorov timescale, and ${\eta }_B$ is the Batchelor scale. The constant $C_B$ is approximately equal to two and is determined by the average value of the compressive principal strain rate, denoted as $<{\tilde{S}}_3>$. Surprisingly, in the present simulations, $C_B$ was found to be in the range of $1.56<C_B<2.2$. This suggests that the compressive strain ${\tilde{S}}_3$ plays a crucial role in determining the presence of the $k^{-1}$ range in the scalar spectra preceding the bottleneck. Furthermore, it can be speculated that the different actions of the three principal strain rates (${\tilde{S}}_i$) on the velocity and scalar fields may be responsible for the formation of scalar gradient structures with a different shape compared to vortical structures. This observation is supported by the absence of a $k^{-1}$ range in the energy spectrum shown in Figure 1, indicating that the dynamics of velocity fluctuations are not controlled by the same mechanisms as scalar fluctuations.

The results of the present simulations are divided into two groups based on the resolution used: one group with a resolution of $N=1154$ for low and intermediate Schmidt numbers, and another group with a resolution of $N=768$ for high Schmidt numbers. The compensated spectra for these groups are shown in Figure 3, with Figure 3a,b representing the $N=1154$ simulations, and Figure 3c,d representing the $N=768$ simulations. The figures demonstrate a significant collapse of the spectra in the dissipation range, starting at $k_b^{\ast }=0.2$. At low wavenumbers, the compensated spectra exhibit a plateau value ($C_{\theta }$) over a certain range of wavenumbers, which is not easily discernible due to the scalar forcing for $k<3$. The spectra from the $N=1154$ simulations (top panels) show a larger number of data points at low $k_b^{\ast }$ and match well with the constant value $C_{\theta }$. The dependence of $C_{\theta }$ on the Schmidt number can be observed in Figure 3b,d, which have scales similar to those in Figure 2 of Su et al. [21]. The transition from the inertial range to the exponential decay range is characterized by a reduction in the slope of $E_{\theta }^{\ast }$, which is required to have a maximum value which is independent of $\mathit{Sc}$.

The range of $k_b^{\ast }$ where the transition from the inertial range to the exponential decay range occurs increases with the Schmidt number. Within this range, the slope of ${\mathcal{E}}_{\theta }^{\ast }$ exhibits a variation from $n=1/3$ to $n=2/3$. These two limits correspond to power laws ($E\approx k^{-m}$) with $m=4/3$ and $m=1$, respectively. By analyzing the data plotted in Figure 3, the variations of $C_{\theta }$ and n with $\mathit{Sc}$ have been evaluated and are shown in the insets of Figure 3a,c. It is worth noting that the value of $m=1$ predicted by Batchelor [22] can be observed for $\mathit{Sc}>2$. Additionally, as discussed earlier, the range of $k_b$ where $E_{\theta }\approx k^{-1}$ increases with an increasing $\mathit{Sc}$.

3.3. Dissipation Spectra and PDFs in the Principal Strain Axes

Figure 4a provides a clear visualization of the amplitude of the spectral bottleneck by plotting the energy and passive scalar dissipation spectra. The collapse of the spectra in Kolmogorov units is shown for both the A and C cases, indicating the range of wavenumbers in which collapse occurs. In physical space, the shape of the spectra at high wavenumbers is determined by the vorticity field and by the scalar gradients. Figure 4a reveals that the largest content of energy and scalar dissipation resides at wavenumbers corresponding to the maximum of the spectral bottleneck. This observation suggests that the size of the vorticity structures is greater than that of the passive scalar gradients. Furthermore, it is expected that the probability distribution functions (PDFs) of the vorticity components (${\omega }_i$) and the passive scalar gradients ($g_{\theta i}=\partial \theta /\partial x_i$) in forced isotropic turbulence are approximately independent of their direction. However, these PDFs are not specifically shown in the figure.

To understand the flow physics in greater detail, the evolution of vorticity components can be analyzed by projecting them along the principal axes of the strain rate tensor. The strain rate tensor eigenvalues, denoted as ${\tilde{S}}_{\lambda }$, provide valuable information about the flow characteristics. The incompressibility condition requires the sum of these eigenvalues to be zero. Among the three eigenvalues, ${\tilde{S}}_1$ is the largest and positive, while ${\tilde{S}}_3$ is the smallest and negative. The sign of the intermediate eigenvalue, ${\tilde{S}}_2$, determines the sign of the determinant $R_{\tilde{S}}={\tilde{S}}_1{\tilde{S}}_2{\tilde{S}}_3$. This determinant, in turn, determines the nature of the vortex structures near a specific point in space. If $R_{\tilde{S}}>0$, the structures are rodlike, whereas if $R_{\tilde{S}}\le 0$, they are sheetlike. From various databases of turbulent flows, including forced and decaying flows, as well as wall-bounded and free-shear flows, it has been observed that sheetlike regions dominate. These sheetlike structures are more unstable and play a significant role in the energy cascade process [23].

Previous studies have also explored the flow topology associated with concentrations of scalar gradients. One notable finding by Kerr [24] was that while vorticity is concentrated in tube-like structures, the scalar gradient and the largest principal rate of strain align perpendicular to these tubes. Kerr [24] also observed that the highest values of the scalar gradient occurred in sheets that wrapped around the vortical tubes. However, the simulations in Kerr’s study were conducted at lower Reynolds numbers than in the present study. In addition to studying the flow topology, the present study also focused on examining the enstrophy and scalar gradient variance budget equations, specifically their alignment with the principal axes of the strain rate tensor. These investigations aimed to provide further insights into the flow dynamics and the role of different terms in the budget equations.

At each point in the computational domain, the strain rate tensor $S_{ij}$ was evaluated, and its eigenvalues ${\tilde{S}}_{\lambda }$ and eigenvectors were determined. These eigenvalues and eigenvectors were used to calculate the respective components of the vorticity vector (${\tilde{\omega }}_{\lambda }$) and of the scalar gradient vector (${\widetilde{g_{\theta }}}_{\lambda }$). To analyze the statistics of these quantities, the normalized PDFs of the generic quantity q were considered, with $q^{\prime }=(q-<q>)/<q^{\prime 2}{>}^{1/2}$, where $<>$ denotes the average value over the entire computational domain. In Figure 4b, the PDFs of the three components ${\tilde{S}}_{\lambda }^{\prime }$ are shown. The compressive strain component (${\tilde{S}}_3^{\prime }$) is depicted in red and exhibits a negative skewness, indicating a concentration of intense negative values. On the other hand, the extensional strain component (${\tilde{S}}_1^{\prime }$) is shown in black and exhibits a positive skewness, indicating a concentration of intense positive values. The magnitude of events in ${\tilde{S}}_3^{\prime }$ is generally stronger than that in ${\tilde{S}}_1^{\prime }$, and this magnitude tends to increase with the Reynolds number. It is worth noting that the intermediate component ${\tilde{S}}_2^{\prime }$ can take both positive and negative values, and consistent with previous findings, the determinant $R_{\tilde{S}}$ is negative for all flow cases, indicating a prevalence of sheetlike structures over rodlike structures.

Table 3. Skewness ($\sigma$) and flatness ($\kappa$) coefficients for the principal strain rates (${\tilde{S}}_{\lambda }$).

Flow Case	${\mathit{\sigma }}_{\mathit{S}1}$	${\mathit{\sigma }}_{\mathit{S}2}$	${\mathit{\sigma }}_{\mathit{S}3}$	${\mathit{R}}_{\tilde{\mathit{S}}}$$(\times {10}^3)$	${\mathit{\kappa }}_{\mathit{S}1}$	${\mathit{\kappa }}_{\mathit{S}2}$	${\mathit{\kappa }}_{\mathit{S}3}$
$A1$	1.657	0.800	−2.119	−1.588	9.201	6.549	10.233
$C6$	1.330	0.472	−1.758	−0.0538	7.259	4.648	7.7777

provides the values of the skewness ($\sigma$) and flatness ($\kappa$) coefficients for the strain fluctuations, as obtained from the PDFs shown in Figure 4b. These coefficients quantify the differences in the distributions of the strain components. The larger value of ${\kappa }_{S3}$ compared to ${\kappa }_{S1}$ indicates a higher amplitude of events for the compressive strain component (${\tilde{S}}_3^{\prime }$) compared to the extensional strain component (${\tilde{S}}_1^{\prime }$). This observation supports the formation of high-amplitude events for ${\tilde{S}}_3^{\prime }$. Additionally, the positive skewness coefficient ${\sigma }_{S2}$ indicates a prevalence of positive events for the intermediate strain component (${\tilde{S}}_2^{\prime }$), contributing to the overall negative value of $R_{\tilde{S}}$. The flatness coefficient ${\kappa }_{S2}$ suggests that the amplitude of events for ${\tilde{S}}_2^{\prime }$ is smaller compared to the compressed and extensional strain components. In Figure 4c, the PDFs of the vorticity components aligned with the ${\tilde{S}}_{\lambda }$ are shown. The symmetric distribution of these components is evident, as indicated by the low values of the skewness coefficients (${\sigma }_{Oi}$) presented in

Table 4. Skewness ($\sigma$) and flatness ($\kappa$) coefficients for the vorticity components in the principal strain axes (${\tilde{\omega }}_{\lambda }$).

Flow Case	${\mathit{\sigma }}_{\mathit{O}1}$	${\mathit{\sigma }}_{\mathit{O}2}$	${\mathit{\sigma }}_{\mathit{O}3}$	${\mathit{\kappa }}_{\mathit{O}1}$	${\mathit{\kappa }}_{\mathit{O}2}$	${\mathit{\kappa }}_{\mathit{O}3}$
$A1$	−0.223	−0.219	−0.207	9.053	8.293	14.128
$C6$	−0.210	−0.261	−0.153	7.157	6.600	10.885

. The negative values of the skewness coefficients, however, highlight a negative skewness of the vorticity components, which is not easily discernible in Figure 4c. It is noteworthy that the symmetry of the ${\tilde{\omega }}_{\lambda }^{\prime }$ events decreases with an increasing Reynolds number, as observed in the table and the figure, as well as in the other flow cases (B and E), which are not shown here.

Table 5. Skewness ($\sigma$) and flatness ($\kappa$) coefficients for the scalar gradient components in the principal strain axes (${\widetilde{g_{\theta }}}_{\lambda }$).

Flow Case	${\mathit{\sigma }}_{\mathit{g}1}$	${\mathit{\sigma }}_{\mathit{g}2}$	${\mathit{\sigma }}_{\mathit{g}3}$	${\mathit{\kappa }}_{\mathit{g}1}$	${\mathit{\kappa }}_{\mathit{g}2}$	${\mathit{\kappa }}_{\mathit{g}3}$
$A1$	−0.208	−0.227	−0.341	14.11	17.54	15.875
$C6$	−0.174	−0.223	−0.289	9.172	12.14	9.274

provides the skewness and flatness coefficients for the normalized scalar gradient components (${\widetilde{g_{\theta }}}_{\lambda }$), and these coefficients show relatively small differences compared to those of the vorticity components (${\tilde{\omega }}_{\lambda }$) presented in Table 4. However, the flatness coefficients for the scalar gradients are generally larger than those for the vorticity components. These findings, along with the spectra shown in Figure 4, suggest that the scalar gradient structures are thinner and more intense compared to the vorticity structures. To gain further insight into this behavior, the enstrophy and scalar gradient variance budgets of the components along the principal axes are analyzed in the subsequent sections.

3.4. Enstrophy Budgets

The budget equation for the time-averaged enstrophy, $\mathsf{\Omega }=<{\omega }_i{\omega }_i>/2$, expressed in terms of the quantities in the principal strain rate axes is given by

$$0=\underset{P_{O\lambda }}{\underbrace{<{\tilde{\omega }}_{\lambda }{\tilde{\omega }}_{\lambda }{\tilde{S}}_{\lambda }>}}+\underset{D_{O\lambda }}{\underbrace{\frac{1}{Re}<{\tilde{\omega }}_{\lambda }{\widetilde{{\nabla }^2\omega }}_{\lambda }>}}.$$

(7)

In this equation, the rate of entropy dissipation $D_{O\lambda }$ was evaluated by projecting on the principal axes of the strain tensor the three components of the Laplacian of the vorticity vector, and the quantities obtained were then multiplied by ${\tilde{\omega }}_{\lambda }$ and summed over $\lambda$. It is noteworthy that an accurate evaluation of the enstrophy dissipation rate requires a proper resolution of the $K^{4}E\left(k\right)$ compensated spectra, which we preliminarily checked.

The ${\tilde{S}}_{\lambda }$ and the ${\tilde{\omega }}_{\lambda }$ events are the contributors to $P_{O\lambda }$, which motivates an in-depth analysis of their correlation. The values of each component of $P_{O\lambda }$ are reported in

Table 6. Enstrophy production (${10}^3\times {\tilde{S}}_{\lambda }{\tilde{\omega }}_{\lambda }^2$), and scalar gradient variance production (${10}^3\times {\tilde{S}}_{\lambda }{\widetilde{g_{\theta }}}_{\lambda }^2$).

Flow Case	$<{\tilde{\mathit{S}}}_1{\tilde{\mathit{\omega }}}_1^2>$	$<{\tilde{\mathit{S}}}_2{\tilde{\mathit{\omega }}}_2^2>$	$<{\tilde{\mathit{S}}}_3{\tilde{\mathit{\omega }}}_3^2>$	$-<{\tilde{\mathit{S}}}_1{\widetilde{{\mathit{g}}_{\mathit{\theta }}}}_1^2>$	$-<{\tilde{\mathit{S}}}_2{\widetilde{{\mathit{g}}_{\mathit{\theta }}}}_2^2>$	$-<{\tilde{\mathit{S}}}_3{\widetilde{{\mathit{g}}_{\mathit{\theta }}}}_3^2>$
$A2$	4.8633	3.4859	−2.4440	−5.705	−0.719	27.159
$C6$	0.1541	0.1342	−0.0744	−2.230	−0.272	10.441

, which shows that the intermediate stress contribution is comparable to that of the extensional one. On the other hand, the component aligned with the compressive strain is the smallest, and its intrinsic negative sign reduces the overall enstrophy production. To understand whether and how events containing ${\tilde{S}}_{\lambda }$ and ${\tilde{\omega }}_{\lambda }$ are correlated, the joint PDF between ${\tilde{S}}_{\lambda }^{\prime }$ and ${\tilde{\omega }}_{\lambda }^{\prime }$ are shown in Figure 4. Small differences between the distributions at high and low Reynolds number have been found, although only the $A1$ flow case is shown. In the analysis of Figure 5, it is important to recall that the skewness of the vorticity components is small, hence the joint PDF are nearly symmetric with respect to the horizontal axis, and the correlation coefficients between ${\tilde{S}}_{\lambda }^{\prime }$ and ${\tilde{\omega }}_{\lambda }^{\prime }$ are thus small. The extensional strain (Figure 5a) and compressive strain (Figure 5c) have rather different distributions. The asymmetric distribution of ${\tilde{S}}_1^{\prime }$ and ${\tilde{S}}_3^{\prime }$ discussed in Figure 4 contributes to the different behavior observed in quadrants I and $II$ compared to quadrants $IV$ and $III$ (see Figure 5a,c). To gain more insight into the difference between the compressed and extensional stress, the contributions of each quadrant to the correlation coefficients are isolated in

Table 7. Contributions to $C_{{\tilde{S}}_l^{\prime }{\tilde{\omega }}_l^{\prime }}$ from each quadrant, evaluated from the joint PDF in Figure 5.

$\mathit{FLOW}$	${\mathit{Q}}_{\mathit{I}}$	${\mathit{Q}}_{\mathit{II}}$	${\mathit{Q}}_{\mathit{III}}$	${\mathit{Q}}_{\mathit{IV}}$	${\mathit{C}}_{{\tilde{\mathit{S}}}_{\mathit{\lambda }}^{\prime }{\tilde{\mathit{\omega }}}_{\mathit{\lambda }}^{\prime }}{10}^{-3}$
$A2$$(\lambda =1)$	$0.20515$	$-0.09259$	$0.09252$	$-0.20454$	$0.54140$
$A2$$(\lambda =2)$	$0.20714$	$-0.11684$	$0.11662$	$-0.20687$	$0.05558$
$A2$$(\lambda =3)$	$0.10140$	$-0.15281$	$0.15176$	$-0.10131$	$-0.97016$

. The first quadrant for the extensional stress contributes almost twice as much as the compressed stress, and the same occurs for the fourth quadrant. The values in the table quantify the higher contributions of quadrant I for the compressed stress compared to that of quadrant $III$ for the extensional stress. Similarly, there is a higher contribution from quadrant I for the extensional stress compared to quadrant $III$ for the compressive stress. This indicates that events with high magnitudes of the stress are associated with the generation of strong vorticity components aligned with them.

The joint PDF between the intermediate strain and the associated vorticity component in Figure 5b exhibits a different distribution compared to the previous cases. However, the contributions from each quadrant are similar to those of the extensional stress, as shown in Table 7. Despite the global lack of correlation between the strain and vorticity components, their interaction contributes to the positive enstrophy production, as seen in Table 6. To further analyze the enstrophy production, the joint PDF between ${\tilde{S}}_{\lambda }^{\prime }$ and ${\tilde{\omega }}_{\lambda }^{2\prime }$ is evaluated. The PDF of ${\tilde{\omega }}_{\lambda }^{2\prime }$ is positively skewed, so the range of ${\tilde{\omega }}_{\lambda }^{2\prime }$ was taken from $-1$ to 7, while the range of ${\tilde{S}}_{\lambda }^{\prime }$ was $-15$ to 15. The normalized joint PDF reveals that the four quadrants are no longer balanced, as shown in

Table 8. Contributions to $C_{{\tilde{S}}_l^{\prime }{\tilde{\omega }}_l^{2\prime }}$ from each quadrant, evaluated from the joint PDF between ${\tilde{S}}_l^{\prime }$ and ${\tilde{\omega }}_l^{2\prime }$.

$\mathit{FLOW}$	${\mathit{Q}}_{\mathit{I}}$	${\mathit{Q}}_{\mathit{II}}$	${\mathit{Q}}_{\mathit{III}}$	${\mathit{Q}}_{\mathit{IV}}$	${\mathit{C}}_{{\tilde{\mathit{S}}}_{\mathit{\lambda }}^{\prime }{\tilde{\mathit{\omega }}}_{\mathit{\lambda }}^{2\prime }}$
$A2$$(\lambda =1)$	$0.16470$	$-0.02388$	$0.08905$	$-0.05288$	$0.17699$
$A2$$(\lambda =2)$	$0.17575$	$-0.05729$	$0.08711$	$-0.05548$	$0.15009$
$A2$$(\lambda =3)$	$0.03248$	$-0.08558$	$0.05512$	$-0.06581$	$-0.06378$

.

Figure 5. Joint PDF between normalized principal strain fluctuations (${\tilde{S}}_{\lambda }^{\prime }$, horizontal axis) and corresponding vorticity fluctuations (${\tilde{\omega }}_{\lambda }^{\prime }$, vertical axis), for flow case $A1$; (a) $\lambda =1$ extensional strain, (b) $\lambda =2$ intermediate strain, and (c) $\lambda =3$ compressive strain. Contours are shown with spacing $\Delta =0.01$ Roman numbers denote the four quadrants.

Figure 5. Joint PDF between normalized principal strain fluctuations (${\tilde{S}}_{\lambda }^{\prime }$, horizontal axis) and corresponding vorticity fluctuations (${\tilde{\omega }}_{\lambda }^{\prime }$, vertical axis), for flow case $A1$; (a) $\lambda =1$ extensional strain, (b) $\lambda =2$ intermediate strain, and (c) $\lambda =3$ compressive strain. Contours are shown with spacing $\Delta =0.01$ Roman numbers denote the four quadrants.

The high positive values of $<{\tilde{S}}_1{\tilde{\omega }}_1^2>$ and $<{\tilde{S}}_2{\tilde{\omega }}_2^2>$ indicate a large contribution from quadrant I (see Table 8). The contribution of quadrant $II$ to $<{\tilde{S}}_1{\tilde{\omega }}_1^2>$ is smaller compared to $<{\tilde{S}}_2{\tilde{\omega }}_2^2>$, resulting in $<{\tilde{S}}_1{\tilde{\omega }}_1^2>><{\tilde{S}}_2{\tilde{\omega }}_2^2>$. On the other hand, the negative value of $<{\tilde{S}}_3{\tilde{\omega }}_3^2>$ is due to the negative contributions from quadrants $II$ and $IV$, which are larger than the positive contributions from the other quadrants. However, none of these contributions is comparable in magnitude to the positive ones for $\lambda =1$ and $\lambda =3$. This explains why $-<{\tilde{S}}_3{\tilde{\omega }}_3^2>$ is smaller than $<{\tilde{S}}_2{\tilde{\omega }}_2^2>$.

3.5. Budget of Scalar Gradient Variance

The budget equation for the scalar gradient variance $\chi =<{\widetilde{g_{\theta }}}_{\lambda }{\widetilde{g_{\theta }}}_{\lambda }>$, is expressed as

$$0=-\underset{P_{G\lambda }}{\underbrace{<{\widetilde{g_{\theta }}}_{\lambda }{\widetilde{g_{\theta }}}_{\lambda }{\tilde{S}}_{\lambda }>}}+\underset{D_{G\lambda }}{\underbrace{\frac{1}{\mathit{Re}\mathit{Sc}}<{\widetilde{g_{\theta }}}_{\lambda }{\widetilde{{\nabla }^2{g}_{\theta }}}_{\lambda }>}},$$

(8)

where the rate of scalar gradient dissipation is evaluated by projecting the three components of the Laplacian of the scalar gradient onto the principal axes of the strain rate tensor, multiplying by ${\widetilde{g_{\theta }}}_{\lambda }$, and then summing over $\lambda$. As for the enstrophy dissipation, we verified that the compensated $K^4{E}_{\theta }\left(k\right)$ spectra decayed at the highest k, hence supporting that the scalar gradient dissipation was properly resolved.

Differences between the terms in Equations (7) and (8) were discussed by Tsinober and Galanti [25], who emphasized the different signs of the production terms. The values presented in Table 6 indeed demonstrate that for the scalar gradient variance, the positive contribution from the compressive strain dominates over the other two components, and the contribution from the intermediate strain is the smallest. Notably, this dominance of the compressive strain contribution is independent of the Peclet number, as evident from the values reported in Table 6 for the high Reynolds number case ($A1$) and for the low Reynolds number case ($C6$).

3.5.1. Rate of Production of Scalar Gradient Variance

The joint PDF between ${\tilde{S}}_{\lambda }^{\prime }$ and ${\widetilde{g_{\theta }}}_{\lambda }^{\prime }$ was evaluated, demonstrating a good independence from the Peclet number. As a result, Figure 6 displays the distributions for flow case $A2$ only. Choosing a flow with $\mathit{Sc}=1$ makes the analysis of the differences between the interactions of strain with vorticity and scalar gradients more instructive, as $\mathit{Sc}$ does not affect the rate of scalar gradient dissipation, see Equation (8). In Figure 5, the three joint PDF differ significantly from each other. However, in Figure 6, the three joint PDF distributions are more similar, when considering the asymmetry of the PDFs of ${\tilde{S}}_1^{\prime }$ and ${\tilde{S}}_3^{\prime }$. The similarity between the joint PDFs shown in panels (c) of Figure 5 and Figure 6 is a preliminary indication of the importance of the compressive strain in the scalar gradient dynamics.

The symmetry of the joint PDFs with respect to the horizontal axis is also observed for the passive scalar gradient components. To quantify their magnitude, the contributions of each quadrant to the correlation coefficient $C_{{\tilde{S}}_{\lambda }^{\prime }{\widetilde{g_{\theta }}}_{\lambda }^{\prime }}$ are reported in

Table 9. Flow case A2: contributions to $C_{{\tilde{S}}_l^{\prime }{\widetilde{g_{\theta }}}_l^{\prime }}$ from each quadrant, as from the joint PDF reported in Figure 6.

$\mathit{\lambda }$	${\mathit{Q}}_{\mathit{I}}$	${\mathit{Q}}_{\mathit{II}}$	${\mathit{Q}}_{\mathit{III}}$	${\mathit{Q}}_{\mathit{IV}}$	${\mathit{C}}_{{\tilde{\mathit{S}}}_{\mathit{\lambda }}^{\prime }{\widetilde{{\mathit{g}}_{\mathit{\theta }}}}_{\mathit{\lambda }}^{\prime }}{10}^{-3}$
1	$0.14934$	$-0.09590$	$0.09574$	$-0.15013$	$-0.95457$
2	$0.11173$	$-0.12432$	$0.12429$	$-0.11172$	$-0.02500$
3	$0.08454$	$-0.16371$	$0.15490$	$-0.08724$	$-11.5133$

. A quick examination reveals that in this case as well, the quantities are uncorrelated. The negative contributions from quadrant $II$ for the $\lambda =3$ component are the largest, resulting in the sum of the three $C_{{\tilde{S}}_{\lambda }^{\prime }{\widetilde{g_{\theta }}}_{\lambda }^{\prime }}$ values being larger than that of $C_{{\tilde{S}}_{\lambda }^{\prime }{\tilde{\omega }}_{\lambda }^{\prime }}$. Table 9 indicates that the four contributions do not differ significantly for $\lambda =2$. Additionally, Figure 6b demonstrates that in this case, the distribution contours in the four quadrants are dissimilar, unlike in Figure 5b. The equilibration is attributed to the positively skewed distribution of ${\tilde{S}}_2^{\prime }$ and to the occurrence of strong ${\widetilde{g_{\theta }}}_2^{\prime }$ events.

The interaction between the ${\tilde{S}}_{\lambda }$ and ${\widetilde{g_{\theta }}}_{\lambda }$ events contributes significantly to the positive value of $-<{\tilde{S}}_3{\widetilde{g_{\theta }}}_3^2>$ in the scalar gradient production (see Table 6). To further understand the large differences among the three contributions to the scalar gradient production, the joint probability density function (PDF) between ${\tilde{S}}_{\lambda }^{\prime }$ and ${\widetilde{g_{\theta }}}_{\lambda }^{2\prime }$ was evaluated. The PDF of ${\widetilde{g_{\theta }}}_{\lambda }^{2\prime }$ exhibited a high positive skewness. Similar to ${\tilde{\omega }}_{\lambda }^{2\prime }$, the joint PDF of ${\tilde{S}}_{\lambda }^{\prime }$ and ${\widetilde{g_{\theta }}}_{\lambda }^{2\prime }$ was evaluated by varying the interval for ${\widetilde{g_{\theta }}}_{\lambda }^{2\prime }$ between $-1$ and 7, and the interval for ${\tilde{S}}_{\lambda }^{\prime }$ between $-15$ and 15. From the normalized joint PDF, the quadrant contributions to each component of $C_{{\tilde{S}}_{\lambda }^{\prime }{\widetilde{g_{\theta }}}_{\lambda }^{2\prime }}$ were evaluated.

Table 10. Flow case A2: contributions to $C_{{\tilde{S}}_l^{\prime }{\widetilde{g_{\theta }}}_l^{2\prime }}$ from each quadrant, evaluated from the joint PDF between ${\tilde{S}}_l^{\prime }$ and ${\widetilde{g_{\theta }}}_l^{2\prime }$, evaluated for $\mathit{Sc}=1$.

$\mathit{\lambda }$	${\mathit{Q}}_{\mathit{I}}$	${\mathit{Q}}_{\mathit{II}}$	${\mathit{Q}}_{\mathit{III}}$	${\mathit{Q}}_{\mathit{IV}}$	${\mathit{C}}_{{\tilde{\mathit{S}}}_{\mathit{\lambda }}^{\prime }{\widetilde{{\mathit{g}}_{\mathit{\theta }}}}_{\mathit{\lambda }}^{2\prime }}$
1	$0.09681$	$-0.04055$	$0.07985$	$-0.06420$	$0.07199$
2	$0.05210$	$-0.06244$	$0.06388$	$-0.06714$	$-0.01360$
3	$0.03026$	$-0.10503$	$0.05949$	$-0.07886$	$-0.09414$

reveals that the distributions in the four quadrants are no longer equilibrated.

The high positive values of $<-{\tilde{S}}_3{\widetilde{g_{\theta }}}_3^2>$ can be attributed to the negative contributions from quadrants $II$ and $IV$ (see Table 10), which outweigh the positive contributions from the other two quadrants. As for the $\lambda =2$ contribution, there is an increased equilibration among the quadrants, resulting in a small value of the correlation coefficient. This explains why $<{\tilde{S}}_2{\widetilde{g_{\theta }}}_2^2>$ reported in Table 7 is much smaller in magnitude compared to the other two components.

3.5.2. Rate of Dissipation of Enstrophy and Scalar Gradient Variance

In the study by Orlandi et al. [4], a “minimal flow unit” was employed to investigate the potential finite-time singularity of the Euler equation by tracking the evolution of relevant statistics. In the viscous case, the evolution of these statistics depended on the Reynolds number. The initial conditions for the “minimal flow unit” consisted of two orthogonal Lamb dipoles positioned at a certain distance from each other. The self-induction effect caused the dipoles to approach and eventually collide. During that time evolution, a range of small vortical structures formed, which were distinct from those observed in the inviscid evolution. The spectra of velocity and passive scalar at a specific time during the viscous collision were found to agree with those of the forced isotropic turbulence. This transition from a vorticity-dominated to a turbulent stage allowed for a connection between the evolution of spectra and the visualization of flow structures in physical space.

In the present paper, before presenting the visualizations of the quantities involved in Equations (7) and (8), it is important to present the averaged values for different combinations of Schmidt and Reynolds numbers. The production of enstrophy and scalar gradient variance, as shown in Table 7, should be balanced by the rate of dissipation of enstrophy and scalar gradient variance, reported in

Table 11. Rate of enstrophy dissipation ($-{10}^3\times D_{O\lambda }$), and scalar gradient dissipation ($-{10}^3\times D_{G\lambda }$).

Flow Case	$-{\mathit{D}}_{\mathit{O}1}$	$-{\mathit{D}}_{\mathit{O}2}$	$-{\mathit{D}}_{\mathit{O}3}$	$-{\mathit{D}}_{\mathit{G}1}$	$-{\mathit{D}}_{\mathit{G}2}$	$-{\mathit{D}}_{\mathit{G}3}$
$A2$	−1.3465	−3.1623	−0.2212	−2.4754	−1.9529	−11.589
$C6$	−0.0537	−0.1482	−0.0056	−0.8477	−0.6935	− 3.9475

. It is important to note that these terms are evaluated from individual fields, so the global balance observed when aggregating a large number of fields may not hold. The values reported in Table 11 highlight that the component aligned with the intermediate strain contributes the most to the rate of enstrophy dissipation, while the component aligned with the compressive strain contributes the least. On the other hand, the rate of scalar gradient variance dissipation is largely attributed to the component aligned with the compressive strain. For the simulations denoted as A and C in Table 1, the ratio $-<{\tilde{S}}_1{\widetilde{g_{\theta }}}_1^2>/<{\tilde{S}}_3{\widetilde{g_{\theta }}}_3^2>$ varied between 0.2 and 0.25 for $0.5<\mathit{Sc}<16$, whereas the ratio $-<{\tilde{S}}_2{\widetilde{g_{\theta }}}_2^2>/<{\tilde{S}}_3{\widetilde{g_{\theta }}}_3^2>$ was approximately 0.025. In contrast, the dissipation of scalar gradient variance shows a dependence on the Schmidt number. For case A, the ratio $D_{G1}/D_{G3}$ ranged from 0.17 to 0.27, and $D_{G2}/D_{G3}$ ranged from 0.12 to 0.20. In the case of C, the variation was minimal, with $D_{G1}/D_{G3}\approx 0.12$ and $D_{G2}/D_{G3}\approx 0.17$.

4. Flow Visualizations

4.1. Vorticity, Strain, and Scalar Gradient

The shape and size of the intense vorticity structures in the flow differ from those where scalar gradients are intensified. This distinction is qualitatively represented by the velocity and scalar spectra. By conducting flow visualizations, it is possible to investigate the association of zones with an intense enstrophy density ($\mathcal{O}={\omega }_i{\omega }_i/2$) and high scalar dissipation density ($\mathcal{C}=g_{\theta i}{g}_{\theta i}$), and regions of intense strain magnitude ($\mathcal{S}=s_{ij}{s}_{ij}$). In forced isotropic turbulence, several regions with intense turbulent activity are observed. However, due to the three-dimensional nature of the flow, the visualizations are not as clear as those presented in Orlandi et al. [4], where a hierarchy of turbulent structures formed in the small region where the Lamb vortices collided. In the present flows, two-dimensional visualizations were performed in a plane defined by the $x_1$ and $x_3$ coordinates, covering the entire computational domain ($-\pi \le x_i\le \pi$). These visualizations revealed several spots with a high enstrophy and scalar gradient variance. To obtain a clear view of these turbulent structures, visualizations were performed for the $C6$ case at a low Reynolds number and $\mathit{Sc}=16$, focusing on a limited area ($-1\le x_i\le 1$). For the $A1$ flow case at a high Reynolds number and $\mathit{Sc}=0.5$, the structures were smaller compared to those in the $C6$ case. Therefore, the visualization area was further reduced to $-0.5\le x_i\le 0.5$.

The normalized fluctuations of the above properties (${\mathcal{O}}^{\prime }$, ${\mathcal{C}}^{\prime }$, ${\mathcal{S}}^{\prime }$) can take both positive and negative values, indicating contributions greater or smaller than their mean values, respectively. These quantities are used to generate normalized PDFs and planar visualizations, as shown in Figure 7. The top row of Figure 7 shows that the three PDFs do not exhibit significant differences, and the probability of obtaining values that are five times greater than their root-mean-square value is relatively small. Furthermore, the positive tail of the PDFs is the same for both the $A1$ and $C6$ flows. The insets in Figure 7b,c indicate that there is only a Reynolds number effect on the negative values, which are associated with the shape and size of the large scales. From the PDFs alone, it can be concluded that they do not provide indications of the influence of the Reynolds number and Schmidt number on the shape and intensity of the structures.

However, such influences can be appreciated through the visualizations of the enstrophy, strain, and scalar gradient variance contours presented in Figure 7. In that figure, negative values (blue contours) occupy a large part of the flow, whereas positive fluctuations (red and green contours) are concentrated in small regions. At low Reynolds number, the contours of ${\mathcal{O}}^{\prime }$ (Figure 7d) and ${\mathcal{S}}^{\prime }$ (Figure 7e) are organized into larger patterns compared to those at high $\mathit{Re}$ (see Figure 7g,h). This difference is even more apparent when considering that the areas which are depicted is four times smaller. The visualizations at low $\mathit{Re}$ of contours of ${\mathcal{S}}^{\prime }$ superimposed to those of ${\mathcal{O}}^{\prime }$ highlight the strong probability that regions with a large ${\mathcal{O}}^{\prime }$ are located near regions with a strong ${\mathcal{S}}^{\prime }$.

To quantify the correlation between the various quantities, the corresponding joint PDFs were evaluated, and the quadrant contributions to the correlation coefficients $C_{a^{\prime },b^{\prime }}$ are shown in

Table 12. Contributions to correlation coefficients $C_{a,b}$ from each quadrant of the joint PDF, related to the visualizations reported in Figure 7.

Flow Case	${\mathit{Q}}_{\mathit{I}}$	${\mathit{Q}}_{\mathit{II}}$	${\mathit{Q}}_{\mathit{III}}$	${\mathit{Q}}_{\mathit{IV}}$	${\mathit{C}}_{\mathit{a},\mathit{b}}$
$A1$$(\mathsf{\Sigma },\mathsf{\Omega })$	$0.44643$	$-0.01391$	$0.08527$	$-0.02023$	$0.49755$
$C6$$(\mathsf{\Sigma },\mathsf{\Omega })$	$0.44226$	$-0.02354$	$0.10089$	$-0.02922$	$0.49038$
$A1$$(\mathsf{\Sigma },\chi )$	$0.13713$	$-0.02998$	$0.05963$	$-0.04336$	$0.12342$
$B3$$(\mathsf{\Sigma },\chi )$	$0.12070$	$-0.03855$	$0.06354$	$-0.05073$	$0.09501$
$C6$$(\mathsf{\Sigma },\chi )$	$0.11348$	$-0.05480$	$0.07497$	$-0.06394$	$0.06970$
$C2$$(\mathsf{\Sigma },\chi )$	$0.12510$	$-0.04297$	$0.07006$	$-0.05744$	$0.09475$

, where a and b denote any two generic quantities. The correlation coefficient between $\mathcal{S}$ and $\mathcal{O}$ is quite large and relatively independent of the Reynolds number. The first quadrant contributes the most, indicating that high and positive values of ${\mathcal{S}}^{\prime }$ are associated with points of intense ${\mathcal{O}}^{\prime }$. Points with small negative values of ${\mathcal{S}}^{\prime }$ and ${\mathcal{O}}^{\prime }$ contribute to $C_{{\mathcal{O}}^{\prime },{\mathcal{S}}^{\prime }}$ to a lesser extent, while the negative contribution to $C_{{\mathcal{O}}^{\prime },{\mathcal{S}}^{\prime }}$ from ${\mathcal{S}}^{\prime }$ and ${\mathcal{O}}^{\prime }$ of opposite sign is negligible. The influence of the compressive strain on the passive scalar leads to elongated patterns of ${\mathcal{C}}^{\prime }$ at low $\mathit{Re}$ and $\mathit{Sc}=16$ (see Figure 7f). This figure, together with Figure 7e, indicates that the correlation coefficient between the ${\mathcal{S}}^{\prime }$ and ${\mathcal{C}}^{\prime }$ events is small, as corroborated by the quadrant contributions in Table 12. Although $C_{{\mathcal{C}}^{\prime },{\mathcal{S}}^{\prime }}$ is reduced, Table 12 emphasizes that the first quadrant contributes the most. Even at high $\mathit{Re}$ and $\mathit{Sc}=0.5$ (see Figure 7i), the patterns of ${\mathcal{C}}^{\prime }$ are still rather elongated, with intense positive (red-colored) events poorly correlated with strong positive ${\mathcal{S}}^{\prime }$ events. Despite this occurrence, the $Q_I$ contribution is greater in the $A1$ flow case than in the $C6$ flow case. From the table, it can be concluded that the correlation coefficient $C_{{\mathcal{C}}^{\prime },{\mathcal{S}}^{\prime }}$ reduces with an increase in $\mathit{Sc}$.

4.2. Enstrophy and Scalar Gradient Variance Budgets Terms

To better understand the differences between the patterns of ${\mathcal{O}}^{\prime }$ and ${\mathcal{C}}^{\prime }$, it is useful to analyze the distributions of the respective terms in their budget equations, Equations (7) and (8), using the PDFs and planar visualizations of the production and dissipation terms. Figure 8 shows the visualizations corresponding to the $A1$ and $C6$ cases. The PDFs of the enstrophy dissipation density ($D_O^{\prime }$, Figure 8b) and of the scalar gradient dissipation density ($D_G^{\prime }$, Figure 8d) are similar and independent of the flow case. The insets emphasize that this similarity holds even for the low-amplitude events. The PDF of the enstrophy production density ($P_O^{\prime }$, Figure 8a) exhibits a good independence with respect to the Reynolds number. Since the positive events are stronger than the negative ones, their skewness coefficients are approximately equal to seven. On the other hand, regarding the scalar gradient production density ($-P_G^{\prime }$, Figure 8c) the positive events are much stronger than the negative ones, indicating a prevalence of the compressive stress over the extensional stress in the formation of strong scalar gradients. The inset in Figure 8c shows a dependence on the Schmidt number for the weak events, but this does not affect the high skewness coefficient, which is equal to 10.5.

At low $\mathit{Re}$, the positive events in Figure 8 are observed to occur within large-size regions, and there is a strong correspondence between negative values of enstrophy production in Figure 8e and zones depleted with enstrophy dissipation in Figure 8f. However, there is a good correlation between the positive events, as evidenced by the dominance of the first quadrant contribution ($Q_I=0.30344$) over the others. Thin and elongated regions with large positive scalar gradient production (Figure 8g) and scalar gradient dissipation (Figure 8h) are also observed, which reduce the correlation coefficient ($Q_I=0.22186$), consistent with the visualizations in Figure 7. At high $\mathit{Re}$, a similar behavior is observed in the bottom figures. The Reynolds number independence is demonstrated by $Q_I=0.30478$ for the correlation coefficient between the terms of Equation (7). There is a slight increase in the correlation between scalar gradient production and scalar gradient dissipation with the Schmidt number, as shown by $Q_I=0.28435$ for the $A1$ case.

4.3. Identification of Enstrophy- and Vorticity-Containing Structures

Figure 9a,b present isosurfaces of the enstrophy dissipation density ($D_O$), of the enstrophy production density ($P_O$), and of the Q criterion [26], all normalized by the mean enstrophy. These visualizations aim to explore the potential correspondence between the terms in Equation (7) and rod- or ribbonlike structures. Visualizations for flow case of $C2$ (low $\mathit{Re}$, $\mathit{Sc}=2$) are shown in Figure 9. The Q criterion, introduced by Jeong and Hussain [26], is a popular method for distinguishing between rodlike and ribbonlike structures in turbulent flows. It is defined as $Q=0.25({\omega }_{\lambda }{\omega }_{\lambda }-2S_{\lambda }{S}_{\lambda })$, evaluated with the components projected on the principal axes. Structures with $Q>0$ are considered rodlike, while those with $Q<0$ are classified as ribbonlike. Three-dimensional visualizations of enstrophy dissipation density ($D_O$), enstrophy production density ($P_O$), and Q, all normalized by the mean enstrophy, can help investigate the correspondence between these terms in Equation (7) and rod- or ribbonlike structures. Figure 9a,b show such visualizations for the case $C2$ at low Reynolds number and $Sc=2$. In these figures, the red surfaces, depicted as circular shapes, represent rodlike structures with $Q>0$, while the blue surfaces correspond to ribbonlike structures with $Q<0$. It is important to note that these visualizations were performed in a small region, and thus the observed volume may not be representative of the entire flow. However, it appears that regions dominated by vorticity occupy a larger volume compared to those dominated by strain. By measuring the volume occupied by the $Q<0$ (blue) and $Q>0$ (red) regions in the entire computational box, it was found that the former accounted for $60.2\%$ of the volume, while the latter occupied $39.8\%$. This observation can be explained by considering that under the constraint $<Q>=0$, the amplitude of the $Q<0$ patches is smaller than that of the $Q>0$ patches. In Figure 9a, the formation of yellow layers, indicating a positive $P_O$, can be observed primarily close to or within the rodlike structures. Figure 9b shows that the negative green regions representing $D_O$ are small, rare, and predominantly occurring in regions with a weak Q.

The distribution of the enstrophy production and dissipation densities can be better understood through contour plots in planar sections. From these plots and Figure 9b, it is evident that the dissipation density is mainly concentrated within the rodlike structures, while only a few spots are located in regions of weak Q. Additionally, locally, the enstrophy dissipation density may take small positive values, which are mostly located outside the rodlike structures. By examining the contour plots in the $x_1-x_3$ plane at the center of the computational box, one can investigate the contributions of each component to the enstrophy balance. In Figure 9c–e, the green regions representing the negative rate of enstrophy dissipation accumulate in elongated patches. The largest contribution comes from the $\lambda =2$ component in Figure 9d, while the contribution from the $\lambda =3$ component in Figure 9e is marginal. In Figure 9h, the component of enstrophy production along $\lambda =3$ is negative and smaller than the other two components. Also for enstrophy production, the $\lambda =2$ contribution (Figure 9g) reaches high values. It is worth noting that the $\lambda =2$ component, as observed in inviscid simulations by Orlandi et al. [4] of interacting Lamb dipoles, leads to an extreme growth of enstrophy, which is a key factor in the eventual formation of a singularity.

4.4. Structures Associated with Scalar Gradients

The visualizations of the terms in the scalar gradient balance equation (Equation (8)) provide insights into the differences with terms in the enstrophy balance equations. While previous studies [5] discussed these differences at low $\mathit{Re}$ and $\mathit{Sc}=1$, they did not emphasize the relative contributions of each component along the principal axes of the strain field. To obtain clearer visualizations, the results of the scalar gradient balance terms are shown for $\mathit{Sc}=2$ instead of $\mathit{Sc}=16$. Figure 9a,b and Figure 10a,b highlight the accumulation of the scalar gradient production term in sheets that are separate from regions with $Q>0$. At $\mathit{Sc}=16$, these sheets of high production become thinner and wider compared to those shown in Figure 10a. Figure 10b illustrates that the scalar gradient dissipation is concentrated in thin layers located where Q is weak. Furthermore, the contributions of the three contributions to the scalar gradient production and dissipation in the principal strain axes are found to be significantly different from the corresponding contributions to the enstrophy production and dissipation. Specifically, Figure 10d,g demonstrate that the contribution of the $\lambda =2$ components is relatively small compared to the other components. The strongest contributions in the scalar gradient balance equation come from the compressive strain, as quantified in Table 6 and Table 11. Figure 10e further illustrates that this component has greater positive contributions to the scalar gradient dissipation compared to the other components. The compressive strain is responsible for the strongest contribution to production, as depicted in Figure 10h and Table 6. It is important to note that in Equation (8), the production term has a negative sign, so the negative green lines in Figure 10h indicate a positive contribution. On the other hand, the extensional strain component gives a negative contribution to the formation of the strong thin layers of the scalar gradient. Figure 10f demonstrates that this component is located in the same regions as those for the enstrophy dissipation in Figure 10c. Indeed, the visualizations provide valuable insights and corroborate the different effects of the strain field on the vorticity and passive scalar gradient components. They allow for a deeper understanding of the spatial distribution and interactions of these quantities, highlighting the regions where they are concentrated and the relationships between them. By examining the visualizations, one can gain a more comprehensive knowledge of how the strain field influences the dynamics of vorticity and passive scalar gradients in the flow.

5. Conclusions

The study of passive scalar advection by turbulent velocity fields has been a subject of extensive research. Previous work, such as the review by Warhaft [27], examined this topic from experimental and theoretical perspectives. Donzis et al. [2] conducted well-resolved direct numerical simulations (DNS) to investigate the influence of Schmidt number variations on scalar statistics, which was also previously reviewed by Antonia and Orlandi [28]. The focus of the present study was to quantify the dependence of the slope of the transitional range in the nondimensional passive scalar spectra on the Schmidt number. This transitional range connects the inertial range (spectral slope $n=0$) of the compensated scalar spectra, at low wavenumbers, to the $n=0$ slope associated with the spectral bottleneck, which occurs at a fixed nondimensional wavenumber as described by Batchelor [22]. In the compensated scalar spectrum at low wavenumbers, there was a constant value $C_{\theta }$ which decreased with an increasing Schmidt number. However, the evaluation of the spectrum at extremely low and high $\mathit{Sc}$ values required a large number of grid points in the simulations to minimize the effect of forcing on the low wavenumber spectra. On the other hand, the dependence of the spectral slope n of the transitional range on the Schmidt number was measured, and it was found that for $\mathit{Sc}>2$, the slope remained fixed at $n=2/3$ [2]. This indicates that the scaling behavior in the transitional range is independent of the Schmidt number for $\mathit{Sc}$ values above a certain threshold.

A further investigation using refined DNS techniques based on codes optimized for GPU clusters, such as CUDA Fortran and OpenACC directives, along with the use of CUFFT libraries for efficient FFT execution, were used to provide valuable insights into the points mentioned. The comparison between the three-dimensional energy dissipation and passive scalar dissipation spectra highlighted the presence of spectral bottlenecks in both quantities. The compensated nondimensional passive scalar spectra exhibited a bottleneck with a maximum value and a wavenumber location that was independent of the Schmidt number. Notably, this nondimensional wavenumber was higher than that associated with the vortical structures. To further investigate the underlying physics, the strain rate tensor $S_{\lambda }$ was evaluated at each computational point. The vorticity and scalar gradient vectors were projected along the principal axes of $S_{\lambda }$, and the joint probability density function (PDF) among them and $S_{\lambda }$ provided insights into the distinct roles played by the components of $S_{\lambda }$ in the generation of vorticity and passive scalar gradient structures. This analysis allowed for a more comprehensive understanding of the dynamics and interactions between velocity and scalar fields in turbulent flows.

Investigating the influence of different vortical structures on passive scalar dynamics is an important research direction that can provide valuable insights. Adding solid body rotation to the flow presents an interesting case to explore using the same code and numerical framework as in the current study. In the presence of high rotation rates, the velocity spectra tend to exhibit a $k^{-3}$ inertial range, indicating the presence of more elongated vortical structures aligned with the rotation axis. Preliminary simulations suggest that for the passive scalar, the spectra maintain a $k^{-5/3}$ range over a wide range of wavenumbers. However, a more detailed analysis similar to the one performed in this paper is necessary to understand this behavior and its implications. Understanding why solid body rotation does not significantly affect the passive scalar spectra can have important implications, particularly in the field of combustion. The complex physics of passive scalar dynamics in rotating chambers can potentially lead to the design of more efficient engines. Therefore, further investigations into the interaction between rotation and passive scalar transport can contribute to advancements in combustion science and technology.