On the role of vortical structures for turbulent mixing using direct numerical simulation and wavelet-based coherent vorticity extraction

The influence of vortical structures on the transport and mixing of passive scalars is investigated. Initial conditions are taken from a direct numerical simulation database of forced homogeneous isotropic turbulence, with passive scalar fluctuations, driven by a uniform mean gradient, are performed for Taylor microscale Reynolds numbers (R λ) of 140 and 240, and Schmidt numbers 1/8 and 1. For each R λ, after reaching a fully developed turbulent regime, which is statistically steady, the Coherent Vorticity Extraction is applied to the flow. It is shown that the coherent part is able to preserve the vortical structures with only less than 4% of wavelet coefficients while retaining 99.9% of energy. In contrast, the incoherent part is structureless and contains negligible energy. By taking the total, coherent and incoherent velocity fields in turn as initial conditions, new simulations were performed without forcing while the uniform mean scalar gradient is maintained. It is found that the results from simulations with total and coherent velocity fields as initial conditions are very similar. In contrast, the time integration of the incoherent flow exhibits its primarily dissipative nature. The evolutions of passive scalars at Schmidt numbers 1/8 and 1 advected by the total, coherent or incoherent velocity suggest that the vortical structures retained in the coherent part play a dominant role in turbulent transport and mixing. Indeed, the total and coherent flows give almost the same time evolution of the scalar variance, scalar flux and mean scalar dissipation, while the incoherent flow only gives rise to weak scalar diffusion.


Introduction
Mixing of passive scalars in turbulent flows is ubiquitous in nature as well as in many technical devices, and has been studied for many decades. Mostly turbulent flow conditions are encountered that guarantee accelerated mixing of all transported quantities like mass, momentum and energy, and thus lead to an enhanced diffusion of all quantities, see, for example, [1]. Indeed, fundamental understanding from a theoretical point of view is a prerequisite to predict, optimize and control mixing, which is of huge practical importance.
There is a large variety of applications, like transport and mixing of pollutants in the atmosphere, or the micromixing of chemical species in reactors, just to mention a few. In addition, studying mixing of passive scalars yields new insights into the properties of turbulence itself, for review we refer to [2].
Small-scale turbulence is of major importance for turbulent flows [3] and in particular with small-scale structures of passive scalars that becomes more important with increasing Schmidt number [4]. The Kolmogorov phenomenology, extended to scalar mixing in [5,6], illustrates the role of Schmidt numbers for the inertial range scaling. Moreover, advanced numerical simulations have facilitated the study of mixing in turbulent flows with high Reynolds number, for example in [7].
Vortical structures are observed in many turbulent flows, see for example, [8][9][10]. These vortices are superimposed on a random background flow, which is structureless and noise-like. This observation motivated the development of methods to extract vortical structures. A prominent example is the Coherent Vorticity Extraction (CVE) method, based on denoising techniques in wavelet space and originally proposed by Farge et al. in [11] and [12] for homogeneous isotropic two-and three-dimensional turbulent flows, which decomposes the flow into two contributions, coherent vortices and an incoherent background flow, using an orthogonal wavelet decomposition of vorticity. Previous studies have shown that the coherent vortices are well represented with a small number of wavelet coefficients and the statistics of the remaining background flow are Gaussian-like. In [13], it was illustrated that the CVE method becomes more efficient with increasing Reynolds number, i.e. a smaller percentage of wavelet coefficients is necessary to reconstruct the vortical structures for higher Reynolds number. The CVE method was applied successfully to study vortical structures in different prototypes of turbulent flows, like sheared and rotating turbulence [14], temporally developing mixing layers [15], or channel flows [16]. For isotropic turbulence we checked the influence of forcing by considering different simulations with different forcing techniques, see for example, [17], [13], [18] or [4]. The results are robust and do not depend on the details of the forcing scheme applied to the velocity field. Based on the wavelet-filtered Navier-Stokes equations a new turbulence model, the Coherent Vorticity Simulation (CVS), has been introduced in [11,19]. The underlying idea is to compute the evolution of coherent vortices deterministically on an adaptive grid, while filtering out the incoherent background flow is sufficient to model turbulent dissipation. For a recent review on wavelet techniques in computational fluid dynamics we refer to [20]. In the context of mixing, wavelet techniques have been applied so far to two-dimensional flows only, where wavelet analyses of passive scalars have been presented in [21]. In [22][23][24][25] the predominant role of vortical structures for mixing of passive scalar has been investigated using CVE for two-dimensional turbulence, with and without chemical reactions. In [22], the coherent flow was represented by 0.8% of the wavelet modes, which retain 99% of energy and 85% of enstrophy. In contrast to the present work, in [22] the mixing of a decaying passive scalar for unit Schmidt number was studied and it was shown that the coherent flow is entirely responsible for the mixing dynamics while the incoherent flow only corresponds to diffusion. In the present work we consider for the first time threedimensional isotropic turbulence with passive scalars in the presence of a mean scalar gradient.
In the present study CVE is applied to investigate the mixing of passive scalars in three-dimensional isotropic turbulence by means of direct numerical simulation (DNS). The aim is to understand the role of vortical structures on turbulent transport and mixing, and to quantify their influence statistically. Indeed, the passive scalar is typically rolled up by the coherent vortex tubes and its fluctuations are produced by the interaction of the velocity field and the mean scalar gradient. Scalar fluctuations are therefore created and advected by the velocity field associated with the vortical structures. The dissipation of scalar fluctuations takes place mainly in regions where iso-scalar surfaces are stretched and folded by the velocity fields. In this paper, we assess the effects of the presence of a range of scales in the velocity and scalar fields using simulations over a limited range of Reynolds and Schmidt numbers (as noted in Table 1). The remainder of the paper is organized as follows. In the next section, the direct numerical simulations with the CVE method are described. The results section is decomposed into three subsections. First, the CVE is applied to the initial condition only (corresponding to one realization of a fully developed turbulent flow) to extract the vortical structures from the total field, which forms the coherent part, the remaining is called the incoherent part. Secondly, the dynamics of the total, coherent and incoherent flows are analyzed. Finally, the transport and mixing of passive scalars in the total, coherent and incoherent flows are investigated. Conclusions of the present work are presented in Section 4.

Numerical simulations and coherent vorticity extraction
In the following we describe the numerical approach and briefly summarize the CVE method. Direct numerical simulations of forced homogeneous isotropic turbulence, with passive scalar fluctuations driven by a uniform mean gradient, are performed for different Taylor microscale Reynolds numbers, R λ = u λ/ν, where ν is the kinematic viscosity, the Taylor microscale λ 2 = u 2 / (∂u/∂x) 2 , and u = u 2 is the root mean square (rms) velocity. Different Schmidt numbers, Sc = ν/D, where D is the molecular diffusivity, are considered using methodologies similar to those described in [4] and [26]. The numerical scheme is based on the Fourier pseudo-spectral algorithm of Rogallo [27] and solves the incompressible Navier-Stokes equations: where u = (u 1 , u 2 , u 3 ) is the velocity field, p the pressure, ρ the fluid density and f = (f 1 , f 2 , f 3 ) a forcing term (which vanishes for decaying turbulence). The domain is a periodic cube with side length 2π . The scalar fluctuation ϕ satisfies the following transport equation: where ∇ = (G, 0, 0) is the mean scalar gradient in Cartesian coordinates, which is responsible for the production of anisotropic fluctuating scalar field.

Coherent vorticity extraction
Here, we introduce the notation for the orthogonal wavelet decomposition of a threedimensional vector field and summarize the principle of the CVE method. For more details we refer the reader to the original papers [12,15]. We consider a vector field, the vorticity, 3 given at resolution N = 2 3J , where J is the number of octaves in each spatial direction. The decomposition of ω into an orthogonal wavelet series yields where the multi-index λ = (j, i x , i y , i z , d) denotes the scale j , the position i = (i x , i y , i z ) and the seven directions, d = 1, 2, ..., 7 of the wavelets. The corresponding index set is and d = 1, . . . , 7}.
Owing to the orthogonality, the wavelet coefficients are given byω λ = ω, ψ λ , where ., . denotes the L 2 -inner product, defined by f, g = [0,2π] 3 f (x)g(x)d x. The coefficients measure fluctuations of ω around scale 2 −j and around position 2πi/2 j in one of the seven possible directions. The N wavelet coefficientsω λ are efficiently computed from the N grid point values of ω using the fast wavelet transform, which has linear complexity [28]. In this study, we have chosen the Coiflet 30 wavelet, which has 10 vanishing moments ( x p ψ(x)d x = 0, p = 0, ..., 9) and is well adapted to represent the current flow simulations. The main steps of the coherent vorticity extraction method are explained in the following; the reader can find more details on the coherent vorticity extraction in [12]: r The wavelet coefficients of the vorticity ω λ are computed by applying the fast wavelet transform to each vorticity component. r A thresholding is applied to the wavelet coefficients,ω λ . Thus, we split the field into two contributions, the wavelet coefficients of the coherent part defined bỹ and those of the incoherent part as the remainder. Z = ω, ω /2 is the enstrophy and N = N x N y N z is the resolution. The threshold ε is motivated from the denoising theory and depends only on the resolution and the variance of the noise, which equals twice the incoherent enstrophy [29]. However, as the enstrophy of the incoherent field is unknown, we estimate a first threshold calculated from the total field and then we apply the thresholding. Subsequently the total field is split with a threshold calculated from the incoherent field. An iterative algorithm developed in [30] can be used for the choice of the optimal threshold. Here we decided to perform only one iteration. Thus, we prefer a good compression rate rather than a perfectly denoised contribution. We checked a posteriori that the incoherent flow is indeed noise-like. The probability density function (PDF) of the incoherent velocity is almost Gaussian and the energy spectrum shows a k 2 slope characteristic for energy equipartition. The noise-like character shows that the incoherent flow is in statistical equilibrium and neglecting its contribution may allow the modeling of turbulent kinetic dissipation rate. In [13] we have shown that modifying the number of iterations in the CVE algorithm has little influence on the results.
r The coherent vorticity ω C and the incoherent one ω I are reconstructed using the fast inverse wavelet transform. The index · C and · I denote the coherent and the incoherent part, respectively. Note that the two fields obtained, ω C and ω I , are orthogonal by construction, which insures a separation of the total enstrophy into Z = Z C + Z I because the interaction term ω C , ω I vanishes. By construction, we also have ω C + ω I = ω.
r The corresponding velocities u = u C + u I are reconstructed using the Biot-Savart law, Since the Biot-Savart operator is not diagonal in wavelet space, i.e. the wavelet decomposition is not orthogonal with respect to the velocity and thus the kinetic energy is split into κ = κ C + κ I + κ , where κ is below 0.09% of the total energy and can thus be neglected [11].
In the results section, the CVE is applied to fully developed turbulent flows, which are statistically stationary, computed by DNS at different Reynolds numbers.
By taking the total, coherent and incoherent velocity fields in turn as initial conditions, three simulations are performed for about five eddy turnover time scales t 0 . For convenience we refer to the flows generated by integrating either the total, coherent or incoherent velocity field as initial condition, the total, coherent and incoherent flows, respectively. An illustration of the three computations is given in Figure 1. The initial eddy turnover time scale is defined as t 0 = L/u , where L denotes the integral length scale of the initial total velocity field. The Kolmogorov microscale of the initial total field is η = (ν 3 / ) 1/4 , where is the mean dissipation rate. The simulations are performed without forcing (such that the turbulence decays, i.e. f ≡ 0 in Equation (1)) while the uniform mean scalar gradient is maintained. The simulation parameters are summarized in Table 1. In order to improve the statistics,  four realizations are taken at different times of the statistically stationary forced flow and used as initial conditions. Subsequently, the ensemble averaged statistics over these four independent realizations are computed.

Initial condition: CVE at t = 0
This part of the paper shows the results of CVE used to obtain the initial conditions at t = 0 for the three runs, i.e. total, coherent and incoherent flows, for each R λ . Tables 2 and 3 show that more than 92% of enstrophy and almost all the energy are retained in the coherent part, while the enstrophy and the energy of the incoherent part are small. The compression rate, corresponding to the number of wavelet coefficients retained in the coherent part, is slightly better for R λ = 240 with 3.2% than for R λ = 140 with 3.6%. We also observe that for increasing R λ , the incoherent enstrophy increases while the incoherent energy decreases. This is important since, in the limit of large Reynolds number, keeping only the coherent vorticity and reconstructing from it the coherent velocity, using the Biot-Savart law, will tend to preserve the total energy. This confirms that CVE becomes more efficient for increasing R λ [13]. Figure 2 shows the probability density functions (PDFs) and the enstrophy spectra of the total, coherent and incoherent vorticity fields. The PDFs of the total and coherent vorticity nearly collapse and display a stretched exponential behavior, while the PDF of the incoherent vorticity is exponential with a width smaller than the ones we can find for the total and coherent vorticities. The enstrophy spectra of the total and coherent fields coincide up to wavenumbers kη ≈ 0.34 for R λ = 140 and kη ≈ 0.36 for R λ = 240. A faster decay of the spectrum is observed for the coherent field in the dissipative range with wavenumbers kη > 0.34 and kη > 0.36 for R λ = 140, 240, respectively. The spectra of the incoherent field contain contributions of all wave numbers, but it is significant only in the dissipative range. Indeed, the spectra of the incoherent field become more important than the ones of the coherent field for wavenumbers kη > 0.89 for both R λ , i.e. R λ = 140, 240, which are however in the dissipation range. These results are in agreement with [12]. Moreover, the total and coherent spectra exhibit an inertial range with a slope close to k 1/3 , which corresponds to the slope of the enstrophy spectrum   according to Kolmogorov's 1941 theory [31]. We can also note that the inertial range is slightly wider for R λ = 240. The spectra of the incoherent part are close to k 4 , which corresponds to the slope of a white noise for kinetic energy, since an energy spectrum proportional to k 2 corresponds to an energy equipartition in the three-dimensional Fourier space. Figure 3 shows the PDFs and the spectra of the total, coherent and incoherent velocity fields. The PDFs of the total and coherent velocities are close to the Gaussian behavior, while the incoherent velocity is Gaussian too but with a strongly reduced variance. The same quantitative behavior as above is observed for the spectra of the kinetic energy that correspond to the enstrophy spectrum divided by a factor k 2 . Indeed, the spectra of the incoherent part are close to k 2 , which corresponds to an energy equipartition.

Dynamics of the total, coherent and incoherent flows
Following the methodology proposed in [32] in the context of two-dimensional turbulence and later applied to three-dimensional turbulence [14], we study the time evolution of the three flows (total, coherent and incoherent) integrated up to five eddy turnover time scales  In isotropic turbulence, the longitudinal integral scale L can be defined as follows: where k is the wave number and E(k) is three-dimensional energy spectrum. Figure 4 shows the evolution of the Taylor microscale Reynolds number R λ (averaged over three directions) for the total, coherent and incoherent flows. We observe that the evolution of R λ for the coherent part coincides with the one of the total part (the curve for the total flow in green is hidden by the red one). Figure 5 shows the rms velocity as a function of time. Again, the coherent part follows the curve of the total flow. The velocity fluctuations in the incoherent part are negligible compared to those in the coherent part. Moreover, the decay of the incoherent flow is faster  than for the total and coherent flows; indeed the dominance of small scales in the incoherent flow is responsible for stronger dissipation and weaker velocities. This can be deduced from Figure 2, which shows that the maximum in the spectrum of the incoherent vorticity appears for higher wavenumbers, so the dissipation is stronger in the incoherent part. The longitudinal integral length scale is computed using one-dimensional energy spectra. The longitudinal integral length scale along direction α is given by where · denotes the space average. Figure 6 reports the time evolution of the average integral length scales for three flows. The total and coherent flow results confirm that the transverse scale is approximately half the longitudinal one. Also, the integral scales for the coherent field are almost the same as those for the total field at any given time. This can be   Figure 6. This is because the smallest scales are dissipated more rapidly, so that at later times there remain only the larger scales. Indeed, as explained in [33], the presence of vortical structures in the total flow prevents the auto-organization of the incoherent part. Conversely, the incoherent flow is able to organize itself into coherent structures if vortical structures are not present in the initial conditions. Furthermore, for N = 512 3 , at short times the incoherent flow contains smaller scales than for N = 256 3 , and larger scales are present for N = 512 3 in the coherent flow.
In order to get further insight into the higher order statistics, the velocity-derivative skewness and flatness as a function of time are plotted in Figure 7. The skewness and flatness of the total and coherent fields have the same tendency but we can remark that the values differ slightly. The incoherent flow exhibits symmetric behavior (skewness around   zero) and flatness which increases sharply with time. This shows that the incoherent flow obeys different dynamics compared to the total and coherent flows. The value of skewness is directly related to the energy transfer between scales. The small value of skewness of the incoherent part indicates that there is only weak nonlinear interaction. The flatness goes rapidly to a value around three, corresponding to the Gaussian-like behavior. At later times the flatness increases, since only sparsely distributed regions of the activity remain. Moreover, for long times the strong increase of flatness is due to the creation of coherent structures in the incoherent flow. Figure 8 shows the normalized vorticity variance ω 2 / ω 2 t=0 and Figure 9 shows the kinetic energy as a function of time. For all three flows, both quantities obey a power-law decay. Indeed for long times, the energy decay presents an algebraical behavior with an exponent larger than one. A more precise determination would require a longer simulation. It is however expected for flows in which the integral length scale is limited by the size of the domain that the decay exponent is equal to two [34]. The coherent part of the field is nearly the same as that of the total field with a similar power law behavior for both cases. Also, the vorticity variance of the incoherent part has an algebraic behavior. Furthermore, a zoom of Figure 9 (not shown here) reveals that the kinetic energy of the coherent part is slightly larger than that of the total field in the power-law region of the decay. At t = 0, the kinetic energy of the total field is larger than the kinetic energy of the coherent part, but the two curves collapse around t/t 0 ∼ 2.5. Indeed, this might be explained by the fact that the total initial field contains by construction more energy than the initial coherent field and the coherent flow dissipates less than the total flow.

Dynamics of passive scalars
The initial spectra of scalar fields are shown in Figure 10, where the Obukhov-Corrsin scale is η OC = ηSc −3/4 . The inertial range with a theoretical power law of k −5/3 is present for a small range of wavenumbers but it is slightly longer for R λ = 240. For small scales, the decay of scalar spectra is stronger for Sc = 1/8. The PDFs of scalar fields are not shown here but they are close to the Gaussian distribution.

The scalar variance equation in isotropic turbulence is
where D ϕ is the molecular diffusivity of the scalar ϕ and · denotes the spatial average. The scalar fluctuation production P is the first term on the right-hand side of Equation (9) and the second term is the mean scalar dissipation rate χ . Generally, large scales are responsible for the production of scalar fluctuations and the dissipation takes place at small scales. The scalar variance produced at large scales is transfered to the small scales. Figure 11 shows the time evolution of the scalar variance ϕ 2 for two Schmidt numbers Sc = 1/8 and 1, transported either by the total, coherent or incoherent velocities. Differences between the total, coherent and the incoherent flows can be observed for t/t 0 > 0.05. Indeed, the total and coherent flows decay slower after this time. We also observe that the incoherent flow is the one which dissipates most. The  coherent flow which contains the vortical structures is dominated by large scales, which are responsible for stronger production, while the structureless incoherent flow causes weak production. Moreover, the higher Schmidt number is responsible for the vertical shift obtained between Sc = 1/8 and Sc = 1 (Figure 11, left and right) corresponding to stronger values of scalar variance for larger Schmidt number. Nevertheless the behavior for these two Schmidt numbers is similar, which indicates a weak influence of Schmidt numbers for this range of values. Finally, the influence of Reynolds numbers can be observed in the decay of the incoherent flow. Faster decay is found for N = 256 3 , since the decay of kinetic energy in Figure 9 is slightly faster for N = 256 3 . Therefore less scalar variance is produced for N = 512 3 (R λ (0) = 240) than for N = 256 3 (R λ (0) = 140). Figure 12 shows the time evolution of the scalar flux − u 1 ϕ for two Schmidt numbers, Sc = 1/8 and 1, transported by either the total, coherent or incoherent flow. The scalar flux is the production term of scalar variance in Equation (9). No significant influence of Schmidt numbers on the scalar flux is found. We can consider that there are also no significant differences between the total and coherent parts. In contrast, the magnitude of the flux for the incoherent part is 10 4 or 10 5 times smaller, hence there is practically no flux in the incoherent flow because of small velocities in this flow. This shows that coherent structures have a predominant role on the scalar flux and there is almost no production in the incoherent flow. The time evolution of the mean scalar dissipation rate χ for the two Schmidt numbers Sc = 1/8 and 1 transported by either the total, coherent or incoherent flow is shown in Figure 13. By plotting Figure 13 in linear scales (not shown here), we find that there is more dissipation in the total part than in the coherent part and that the incoherent part decays faster. Furthermore, this figure confirms the observation and explanation for scalar variance. There is less scalar variance to dissipate in the incoherent flow and at small times the little variance, which must be dissipated at small scales, dissipates rapidly. The few remaining scales then produce scalar variance, hence the slow decay at long times. The small decay for the total and coherent flows is explained by the fact that the turbulence is decaying. Hence, there is less production of scalar variance, and consequently there is less scalar variance to dissipate. We can also note that the value of scalar dissipation for different contributions is insensitive to molecular diffusivity, which is in agreement with [35]. Figure 14 shows the time evolution of production-dissipation ratio of the scalar P / χ for two Schmidt numbers, Sc = 1/8 and 1, transported by either the total, coherent or incoherent flow. We observe in Figure 14 that for t/t 0 > 1 we obtain a quasi stationary state, as the transfer from large scales to small scales is well balanced. Furthermore, there is almost no production in the incoherent flow since the flux is weak and the dissipation decreases strongly with time.

Conclusions
In this study, we performed several direct numerical simulations to investigate the influence of vortical structures on transport and mixing of passive scalars in isotropic turbulent flows. In order to capture the vortical structures, we applied the CVE technique to fully developed turbulent vorticity fields before running the simulations. The effects of Reynolds and Schmidt numbers were studied by simulating scalar fields with Sc = 1/8 and 1 in decaying isotropic turbulence with well-developed initial states at R λ = 140 and 240.
Compared to [22], the compression rate differs between two and three dimensions and is stronger in the current three-dimensional study. We found that CVE becomes more efficient for increasing Reynolds number, which confirms the results obtained in [13]. Moreover, the total and coherent flows give almost the same time evolution, while time integration of the incoherent flow is mainly of dissipative nature as found in [14] and can thus be neglected for modeling turbulent dissipation. Furthermore, we observed that the behavior of the passive scalar advanced in time is almost the same for the total and coherent flows, which are thus the parts where the scalar is transported and mixed. In contrast, the incoherent flow is dominated by diffusion for small times because of the presence of structureless fields with small velocities, and reflected by the absence of production of scalar variance. Therefore, we conclude that the vortical structures are responsible for turbulent transport and mixing of passive scalars. Moreover, strong shear regions, which are important for scalar mixing, contain vorticity and as the CVE is based on vorticity, these regions are extracted and retained in the coherent part. In addition, for the Reynolds and Schmidt number range examined in this paper we have found that the observed mixing behavior has no significant dependence on either of these parameters.
These results motivate the use of CVS to model scalar mixing in turbulent flows as the vortical structures, computed deterministically in CVS, are responsible for the transport and mixing. The investigation of mixing at higher Schmidt numbers is an interesting problem and may show a change in the dynamics of scalar mixing in the incoherent flow.