senior researcher, Francois G Schmitt

Enrico Calzavarini

Sergio Chibbaro

associate professor, Lipo Wang

A Wang

The present work focuses on the statistical properties of reactive scalars undergoing reversible chemical reactions in incompressible turbulence. Theoretical analysis about the statistical properties of scalars at different order of moments were carried out based on appropriately proposed approximations and models. The theoretically derived results were then compared with numerical results obtained by direct numerical simulation (DNS). In the direct numerical simulation, the spatial derivatives were mainly approximated by using a pseudo-spectral method, since the turbulent velocity and scalar fields are generally of periodic boundary conditions. For the special configurations in which the boundary condition is not periodic, a finite difference method with fine schemes was used to approximate the spatial derivatives. The numerical time integration was implemented by a third order Runge-Kutta scheme.

In statistically steady homogeneous isotropic turbulence, the chemical species retained close to a dynamical equilibrium state sustained by random large-scale reactant sources. A competition exists between the chemical reaction that tends to dump reactant concentration fluctuations and enhance their correlation intensity and the turbulent mixing that on the contrary increases fluctuations and remove relative correlations. A unique control parameter, the Damköhler number (Da θ ) that can be constructed from the scalar Taylor micro-scale, the reactant diffusivity and the reaction rate, was found to characterize the functional dependence of fluctuations and correlations. Such a dependence was validated in a variety of conditions, i.e., at changing the reaction order, the Reynolds and the Schmidt numbers. The larger is the Damköhler number, the more depleted are the scalar fluctuations as compared to the fluctuations of a passive scalar field in the same conditions, and the more intense are the correlations. A saturation in this behaviour was observed beyond Da θ O(10). We provide an analytical prediction for this phenomenon which is in excellent agreement with the direct numerical simulation results.

摘要

Résumé

Dans cette thèse, les propriétés statistiques des scalaires réactifs subissant des réactions chimiques réversibles en turbulence incompressible ont été étudiées au moyen de simulations numériques directes et d'analyses théoriques. Les cas étudiés incluent les réactions proches et fortement déviées des états d'équilibre chimique, en écoulement turbulent homogène et isotrope, et inhomogène anisotrope. Des analyses théoriques des propriétés statistiques des scalaires pour différents ordres de moments ont été effectuées sur la base d'approximations et de modèles. Les résultats théoriques ont ensuite été comparés aux résultats numériques obtenus par simulations numériques directes. LIST OF TABLES 6.1 Non-dimensionalized parameters for the simulations: Re = u H/ν is the Reynolds number based on large scale, where u is the single-component root-mean-square velocity, H is the length of the domain, ν is the viscosity; Re λ = u λ/ν is the Taylor scale λ based Reynolds number; Sc is the Schmidt number (ν/D); N 3 is the number of total grid points; |k| max •η is the resolution condition, where |k| max is the maximum amplitude of wave number kept by the dealiasing procedure, η if the Kolmogorov lenght scale; τ η is the Kolmogorov time scale; Γ = Da 1 /Da 2 , with Da 1 and Da 2 as the Damkholer numbers for forward and backward reactions respectively; L I is the integral length scale; T I is the integral time scale; ∆t is the numerical time step. . . . . . . . . . . .

In the non-homogeneous anisotropic turbulence: the turbulent Kolmogorov flow (TKF), it was found that the mean velocity profile has the same form with the forcing. The only non-zero shear stress term is proportional to the cosine function, and the normal stress components all involve a square cosine expression. The normal stresses are never equal, showing that as expected the turbulence is anisotropic. It was also shown that a quadratic nonlinear constitutive equation can be proposed for this flow, involving a linear term and two nonlinear terms in the form of traceless and symmetric tensors. For about half of the flow domain, the linear term is dominating.

Whereas for the vanishing mean velocity gradient regions, only one non-linear term remains non-zero and becomes constant. The reversible reactions discussed in the homogeneous isotropic case were also studied with the background flow as the TKF. The theoretical predictions about the dependence of the fluctuations and the correlations of the reactive scalars on the Da θ in homogeneous isotropic flow are found to work also for the TKF case. It indicates that, in the quasi-equilibrium state, the dependence of the statistical properties of reactive scalars on the scalar diffusion and reaction rate is weakly influenced by the background velocity field.

By adopting Dirichlet boundary conditions for the scalars, a combustor like configuration was also proposed.

In such a configuration, the reacting system is strongly deviated from the chemical equilibrium. The entire flow consists of two buffer layers and a bulk region, in which the scalar mean gradient is maintained by the boundary conditions to provide non-negligible source of the chemical reaction. A theoretical model based on the PDF of the passive scalar can satisfactorily predict the mean and fluctuation of reactants, if the associated Damköhler ii number (Da 1 ) is sufficiently large. The correlation coefficient between the scalar quantities are determined by two counteracting effects, the turbulent mixing and the chemical kinetics. Under the non-equilibrium condition with strong chemical sources, the chemical kinetics also plays important roles in determining the scalar energy spectra. For large Da 1 , the spatial distribution of the forward reaction and net reaction assume stripe like structures, making the scalar field more intermittent. Consistently, large Da 1 will shift the scalar energy from the small wave number range to the large wave number range, which is different from the quasi-equilibrium case.

All the works carried out in this thesis are devoted to the numerical and theoretical explorations about reactive scalars in incompressible turbulence of different configurations. Our findings also suggest new ideas for future studies, which are discussed in the conclusions.

Introduction

La turbulence, domaine complexe mais scientifiquement important, existe dans de nombreuses situations, telles que les écoulements dans les moteurs automobiles, les réacteurs chimiques, les écoulements environnementaux, y compris les océanographiques et météorologiques. Une caractéristique essentielle des écoulements turbulents est que le champ de vitesse du fluide varie de manière significative et irrégulière à la fois en position et en temps.

En conséquence, les vitesses sont apparemment aléatoires et chaotiques.

Cependant, du point de vue statistique, certaines propriétés intrinsèques et universelles de la turbulence peuvent être trouvées, par exemple la loi -5/3 décrivant le processus de cascade d'énergie régissant les tourbillons de turbulence. Dans la turbulence, il y a également des scalaires advectés.

Certains sont appelés scalaires passifs car ils n'apportent pas de rétroaction sur l'écoulement (généralement incompressible). Certains sont appelés scalaires réactifs car ils impliquent une réaction en plus du transport par l'écoulement. L'objectif principal de cette thèse est d'étudier les propriétés statistiques de scalaires réactifs en turbulence incompressible, possédant une réaction sans influencer l'écoulement. Dans le monde réel, des exemples de scalaire réactif dans une turbulence incompressible peuvent être les concentrations de réactifs chimiques (sels nutritifs) ou le phytoplancton dans l'océan.

Simulations numériques directes

Les équations aux dérivées partielles à résoudre numériquement dans cette thèse sont les équations de Navier-Stokes pour la vitesse et les équations de advection-diffusion pour les scalaires non réactifs, ajoutés à des temres de réaction pour les scalaires réactifs. Les simulations numériques directes ont été mises en oeuvre en résolvant ces équations aux dérivées partielles, dans différentes configurations, en utilisant la méthode pseudo-spectrale ou la méthode des différences finies. Lorsque les quantités ont des conditions aux limites périodiques, la méthode pseudo-spectrale est préférée en raison de sa grande précision dans l'approximation des dérivées spatiales. Les méthodes pseudo-spectrales utilisées ici ont adopté les séries de Fourier comme fonctions orthogonales. Pour éviter l'erreur de crénelage due à la convolution des modes à grand nombre d'ondes introduite par le terme non linéaire, la technique de désaliasing lisse a été utilisée. Pour les configurations dans lesquelles les conditions aux limites des scalaires ne sont pas périodiques, une méthode aux différences finies avec des schémas d'ordre élevé (schéma de huitième ordre à pas spatial vers l'amont et schéma de centre du dixième ordre pour approximer respectivement les dérivées spatiales du premier et du second ordre) a été utilisée.

L'intégration temporelle numérique a été implémentée de manière explicite par un schéma Runge-Kutta du troisième ordre. Les analyses statistiques des grandeurs d'intérêt ont été mises en oeuvre dans un état stationnaire statistique, qui a été maintenu par forçage numérique pour la vitesse et les scalaires. Dans la configuration anisotrope, le forçage des scalaires se faisait via un gradient moyen, tandis que dans la configuration isotrope, la vitesse et les fluctuations scalaires sont soutenues par un forçage portant sur l'ensemble des nombres d'onde.

Résultats et discussion

Nous avons d'abord étudié le cas des réactions chimiques réversibles en turbulence homogène et isotrope. La réaction fait intervenir trois réactifs R 1 , R 2 et P, dont les fluctuations ont été maintenues par un forçage à grande échelle avec une amplitude constante et une phase aléatoire, subissant une réaction réversible. En outre, un scalaire passif non réactif T est également considéré pour comparaison. Le système réactionnel se trouve dans un état de quasi-équilibre chimique, qui est maintenu par la compétition entre la réversibilité de la réaction, qui impose un équilibre chimique global à l'état statistiquement stationnaire, et la variabilité du forçage scalaire, qui introduit une stochasticité dans les champs scalaires. L'amplitude du forçage scalaire a été contrainte pour générer de petites fluctuations scalaires par rapport à la quantité moyenne globale de scalaire (~10% de la quantité moyenne), sinon la quantité scalaire totale pourrait être négative, ce qui n'a aucun sens pour la concentration de réactif. Dans ce cas, la vitesse de réaction nette globale était faible dans l'état d'équilibre dynamique, et donc trop faible pour modifier les lois d'échelle et la distribution statistique des scalaires. En conséquence, les scalaires réactifs ont montré la même loi d'échelle (spectre d'énergie en Fourier) que celui du scalaire non réactif. Les fluctuations des scalaires réactifs et non réactifs normalisés par leurs écarts types obéissaient à des distributions gaussiennes. A partir des spectres de cohérence entre scalaires réactifs, il a été trouvé que les coefficients de corrélation entre scalaires réactifs sont indépendants du nombre d'onde (échelle de longueur). Un tel résultat suggère que le coefficient de corrélation entre deux scalaires réactifs est le même que celui entre leurs gradients. L'effet de la réaction sur les propriétés statistiques des scalaires a été principalement mis en évidence dans les corrélations et l'amplitude des fluctuations, sur lesquelles nous avons réalisé non seulement des simulations numériques mais aussi des approches théoriques. Pour les analyses de modélisation des coefficients de corrélation des scalaires réactifs, quatre conditions ou hypothèses préalables sont nécessaires: 1. l'état stationnaire statistique; 2. la linéarisation des termes de la réaction, qui est viable à condition que les fluctuations des scalaires soient faibles par rapport aux quantités moyennes; 3. le fait que le coefficient de corrélation entre deux scalaires réactifs est le même que celui entre leurs gradients; 4. l'hypothèse que le forçage d'un scalaire ne peut pas être fortement corrélé avec celui d'un autre scalaire, ce qui est garanti par l'implémentation du forçage temporel aux scalaires. Puis dans le cadre de cette modélisation, à partir des équations gouvernant les scalaires réactifs, des dérivations mathématiques ont été effectuées pour obtenir des équations approchées sur les covariances et les auto-variances des scalaires réactifs, ce qui a finalement conduit à des expressions analytiques montrant que les coefficients de corrélation des scalaires réactifs sont des fonctions d'un paramètre de contrôle unique: le nombre de Damköhler (Da θ ). Da θ est construit comme le rapport entre l'échelle de temps de diffusion scalaire à travers un domaine de la taille de la micro-échelle scalaire de Taylor (λ θ 2 /D) et l'échelle de temps de réaction chimique (τ r ). De plus, Da θ s'est avéré être Le rôle clé de Da θ implique également l'importance de la micro-échelle scalaire de Taylor λ θ dans le mélange des espèces chimiques. Tout d'abord, il a été remarqué que λ θ ne varie pas significativement pour les différents champs scalaires R 1 , R 2 , P et le champ scalaire passif de référence T. De plus, nos résultats numériques ont montré que λ θ est à peu près inversement proportionnel à la racine carrée du nombre de Schmidt, surtout sous la condition d'un grand nombre de Reynolds. Une telle relation implique que Da θ peut être considéré comme le rapport de la plus grande échelle de temps de l'écoulement turbulent (grand temps de rotation des tourbillons) à l'échelle de temps typique associée au processus chimique (τ r ). Une telle idée a été témoignée dans une série de simulations où les champs scalaires sont advectés par l'écoulement turbulent spatialement filtré, ne conservant que de grands tourbillons. Il a été constaté que la micro-échelle de Taylor scalaire λ θ et les coefficients de corrélation des scalaires réactifs dans les champs de vitesse de seulement quelques modes les plus grands sont les mêmes que ceux convectés par tout le flux. Cette dernière observation peut être d'un intérêt pratique pour l'estimation des régimes atteints par les réactions biogéochimiques à petite échelle dans l'océan. Par exemple, si l'on considère un champ réactif comme la concentration locale de phytoplancton dans l'océan. Les échelles de temps utilisées pour construire le nombre de Da θ sont: le grand temps typique de rotation des tourbillons de l'écoulement turbulent tridimensionnel, qui est normalement de l'ordre de plusieurs heures, et le temps de croissance typique de la population, c'est-à-dire du ordre d'un jour. Cela conduit à une valeur de Da θ <1, ce qui signifie qu'à de petites échelles, plus précisément à des échelles où le flux océanique peut être approché comme un écoulement turbulent tridimensionnel avec une cascade d'énergie directe, la concentration de phytoplancton peut être considérée en toute sécurité comme un champ scalaire passif.

Dans le cas isotrope homogène, l'analyse théorique des propriétés statistiques des scalaires réactifs a été faiblement affectée par les propriétés à petite échelle de l'écoulement turbulent. Cela nous a incité à considérer la turbulence anisotrope non homogène, telle que l'écoulement turbulent de Kolmogorov (TKF). L'écoulement de Kolmogorov (KF) est un exemple d'écoulement de canal turbulent sans frontières, qui est statistiquement non homogène dans une direction et anisotrope. Généralement, le système KF correspond aux équations de Navier-Stokes étudiées dans un domaine périodique, à pression constante, et forcées par un forçage sinusoïdal. Audessus d'un nombre de Reynolds critique, l'écoulement devient turbulent et nous le désignons par écoulement turbulent de Kolmogorov (TKF).

Tout d'abord, nous avons acquis une compréhension plus approfondie du TKF forcé dans la direction x par une force sinusoïdale dépendant de la coordonnée z. On a constaté que le profil moyen de la composante de vitesse alignée sur le forçage a la même forme, avec un amortissement d'un facteur κ, par rapport à la valeur de vitesse moyenne calculée à partir du terme de forçage. On a trouvé que la valeur de κ était comparable à celle rapportée dans les travaux existants et qu'elle augmentait faiblement avec le nombre de Reynolds, indiquant une possible saturation asymptotique à très grand Re. Le seul terme de contrainte de cisaillement non nul est la covariance entre les composantes de vitesse le long des directions x et z.

Plus important encore, cette contrainte de cisaillement non nulle est proportionnelle à l'intégrale de la composante de vitesse moyenne le long de la direction z, qui est une fonction cosinus. Un tel résultat numérique a validé la convergence de la simulation, puisqu'il peut être théoriquement obtenu en faisant la moyenne de la composante x des équations de Navier-Stokes. Les composantes de contrainte normales impliquent toutes des fonctions cosinus carrées. Ces derniers ne sont jamais égaux, montrant que, comme prévu, la turbulence est anisotrope. En additionnant les trois contraintes normales, on obtient une expression en cosinus carré de l'énergie cinétique, dont les coefficients étaient en bon accord avec les valeurs rapportées par les travaux déjà publiés.

Ensuite, nous avons également étudié les fermetures pour l'écoulement turbulent de Kolmogorov, qui proposent des équations constitutives entre le tenseur de contrainte anisotrope et le tenseur de vitesse de déformation moyenne. Premièrement, selon l'hypothèse de viscosité turbulente de Boussinesq, un modèle de fermeture linéaire, qui suggère une relation linéaire entre les deux tenseurs. Un tel modèle linéaire suggère que l'alignement entre les deux tenseurs est de 1. Les résultats numériques ont montré que l'alignement est supérieur à 0,9 pour environ 46% du volume, où le gradient moyen de vitesse est grand. Cela indique que pour environ la moitié du volume avec un fort gradient de vitesse moyenne, la relation Les amplitudes des taux de réaction directe et inverse ont d'abord été examinées. On a trouvé que l'amplitude de la vitesse de réaction directe pouvait être significativement plus grande que celle de la vitesse de réaction inverse, en particulier lorsque Γ est grand. Ceci indique que le système réactif dans une telle configuration était en effet fortement dévié de l'état d'équilibre chimique. Pour une meilleure compréhension de cette configuration nouvellement développée composée de deux couches tampons et d'une région de mélange, les moments du scalaire passif ont été explorés à la fois numériquement et théoriquement. Sous la condition d'état stationnaire statistique et en utilisant l'hypothèse de longueur de mélange, l'expression analytique du profil moyen du scalaire passif peut être dérivée de l'équation d'advection-diffusion. Ces analyses théoriques suggèrent que le transport scalaire à la fois dans les couches tampons et dans la région centrale de mélange peut être considéré comme des diffusions scalaires turbulentes avec différentes diffusivités turbulentes. Surtout, les prévisions théoriques étaient en excellent accord avec les résultats numériques.

Avec un Γ grand, la fonction de densité de probabilité (PDF) de la vitesse de réaction nette (R net ) culmine pour la valeur R net = 0 et devient asymétrique et plus étendue pour les grandes valeurs de R net . Alors que lorsque Γ est proche de 1, la PDF culmine à une valeur modérée de R net .

Une telle différence doit être causée par le mécanisme chimique, ce qui peut être vu plus clairement à partir de la distribution spatiale des vitesses

Contributions scientifiques

Les contributions scientifiques de cette thèse peuvent être résumées comme suit: (1) le sujet de cette thèse, qui est une combinaison de réaction et de scalaire passif dans des écoulements turbulents, est nouveau; ( 2 [START_REF] Jayesh | Probability distribution of a passive scalar in grid-generated turbulence[END_REF]). For legibility purpose, the upper curve has been shifted by 2 decades with respect to the lower one.

Note that in the labels, the θ and u are the fluctuations and θ and u are the corresponding rms. The deviation of PDF from Gaussian curve for the scalar is larger than that of the velocity, 3.3 The DNS results of the energy spectra of the reaction rates (γ 1 R 1 R n 2 -γ 2 P ) of the first order (solid line) and second order (dashed line) reactions in the form of (3.40). k = |k| is the amplitude of wave number vector. The dealiasing is implemented as smooth dealiasing, with a filtering 6.11 The mean (solid lines) and r.m.s. (dashed lines) of (a) R 1 compared with T ; (b) P under the conditions of different Γ as functions of position (z). The main panel of (c) shows the mean profile of P normalized by its maximum, whose function as Γ is plotted in the inset plot. There is a prefect superposition for all Γ values. In all the plots, the vertical dotted lines mark the interfaces between the buffer layers and the bulk region. . . . . . . . . . . . . . . . . . . . . . . . 113 6.12 The scalar mean: (a) R 1 z and (b) P z , as a function of z obtained from theoretical analysis (dashed lines) based on Eq. (6.25) and DNS (solid lines with the same colors). The grey dashed lines are from the theoretical prediction at infinitely large Da 1 according to Eq. (6.20). We see that the prediction for Γ = ∞ is close to the curves for Γ = 100 and also the predictions for large Γ are close to the DNS results when Γ=10, 30, 100. The vertical dotted lines mark the interfaces between the buffer layers and the bulk region.

List of Tables

λ 0 = 15ν/ 0 u 0 , ν is the kinematic viscosity, u 0 = 1 3 i u 2 i is the global root-mean square of single component velocity, 0 = ν 2 i j (∂ i u j + ∂ j u i ) 2
is the global energy dissipation rate, • denotes the average in time and all over the domain. N 3 is the grid size; η 0 = (ν 3 / 0 ) 1/4 is the global Kolmogorov scale; |k| max • η 0 is the resolution condition; T total is the total simulation time and T l is the large eddy turnover time, i.e. T total /T l denotes the number of large

Objectives of the thesis

In this thesis, direct numerical simulations are implemented to numerically study reversible chemical reactions in incompressible turbulence. The main focus is on the fundamental properties of reactive scalar mixing, in comparison with a passive scalar, under a statistically stationary state. The flow statistics, from global to scale-dependent features, are studied in details and theoretically modelled.

The objectives of this thesis are the following:

• Obtain numerically statistical laws characterizing the fluctuations of reactive quantities for second order reaction or higher orders;

• Estimate the scaling properties of the reactive and passive scalars, using Fourier power spectra and structure functions for the intermittent properties;

• Perform a parametric study by varying e.g. the reaction order or the Damkholer number and study their effects in a statistical sense;

• Check different configurations to try to obtain the universal laws.

Organization of the thesis

The rest of this thesis is organized as follows. Chapter 1 introduces some basic concepts and classical theories about the incompressible turbulence and passive scalar fields convected. Chapter 2 presents some fundamentals and models about the chemical reaction, and it reviews the literature discussing the dynamics of reactive scalars in incompressible turbulence. Chapter 3 introduces some basic ideas of the numerical methods in computational fluid mechincs, in particular the pseudo-spectral and the finite difference approaches. It also discusses how they were implemented in this thesis. The subsequent chapters are about the original works of this thesis. Chapter

Introduction

Turbulence, scientifically important and challenging, exists widely in various situations, such as flows in automotive engines, chemical reactors, environmental flows including oceanographical and meteorological flows. It can be observed from the clouds on the sky, a turbulent jet ejected by an aircraft generating wakes, plumes spiraling up and spreading, or the whirls on the sea surface visualized from a satellite. An essential feature of turbulent flows is that the fluid velocity field varies significantly and irregularly in both position and time. Thus they are seemingly random and chaotic [START_REF] Pope | Turbulent Flows[END_REF]. From the statistical point of view, physicists have been focusing on finding some intrinsic and universal properties of turbulence, although turbulent motions can be of different configurations, different boundary conditions, and generated by various types of forcing (such as string to a cup of water). In turbulence studies, turbulent flows are considered to be composed of eddies of different sizes, which can span a very large range of scales. For example, the largest eddies in the ocean are of the order of kilometers or even larger, but the smallest eddies are invisible to human eyes. Typically, the external forcing is exerted at the large scales. Because of the energy cascade, the turbulent kinetic energy flows hierarchically towards the smaller scales till the smallest end where the kinetic energy is dissipated to heat. According to a picture originally proposed by Kolmogorov (1941a,b), at the scales much smaller than the domain or forcing scale but larger than the smallest eddies, the information of boundary or forcing from largest scale is lost, and the effect of dissipation is negligible. Thus universality is believed to exist in the properties of turbulent motions of these scales because they are fully determined by the inertial interaction between eddies. This range of scales is called inertial range and it is the main focus of most of fundamental studies in turbulence theories. This chapter is an introduction about some fundamentals concerning incompressible turbulence and nonreactive scalars transported in it. It is organized as follows. In section 1.1, the mathematical description of the Introduction turbulent motion, Navier-Stokes equations, is first introduced. Then the idea of energy cascade in turbulent and the classical K41 theory quantifying the scaling behaviors in turbulent flow are introduced. Section 1.2 is about passive scalars in incompressible turbulence, focusing on the governing equations, statistical properties, scaling behaviors and intermittency of passive scalar. In next chapter, Chapter 2, the fundamentals about the reactive scalars in turbulence, which is the main topic of this thesis, are presented.

Turbulent flow

Navier-Stokes equations

Lagrangian and Eulerian description

To derive the equations fully describing the motion of turbulent flows, the clarification of the framework under which the motion of fluids is described is needed. There are two ways to describe fluid motion. One is the Lagrangian description, in which fluid particles are followed as they move through the flow field. The other is the Eulerian description, in which the properties of the flow field are monitored using a fixed reference frame. In fact, most studies about fluids mechanics are based on the Eulerian description, because the global properties of the flow field or the properties at specific region are of physicists' interests, instead of any particular fluid particles. However, the understanding of the Lagrangian description is also necessary because the derivation of the momentum equation about fluid motion is exactly the application of Newton's second law to the fluid particles.

The Lagrangian description is based on tracing the motion of fluid particles. The position of a fluid particle is the function of current time t, the reference time t 0 and the position of this particle at the reference time r 0 : r(t; r 0 , t 0 ). And the velocity of this particle is the temporal derivative of r(t; r 0 , t 0 ): u = dr(t; r 0 , t 0 )/dt. In such a track-dependent description, a property of the fluid field F is described as function of the position and time:

F = F [r(t; r 0 , t 0 ), t].
The Eulerian description focuses on properties of flow field at the locations of interest, in which the fluid particles are not discriminated. A property of the fluid field F is the function of position and time: F = F (x, t). F [r(t; r 0 , t 0 ), t] and F (x, t) are two descriptions of one quantity. When r = x, F [r(t; r 0 , t 0 ), t] = F (x, t).

(1.1) dr

3 dt + ∂F ∂t = ∂F ∂x 1 u 1 + ∂F ∂x 2 u 2 + ∂F ∂x 3 u 3 + ∂F ∂t = (u • ∇)F + ∂F ∂t = D Dt F (x, t). (1.2)
Eq. ( 1.2) builds the connection between Lagrangian and Eulerian description of fluids, in which D Dt is called material derivative.

Continuum approximation and Reynolds' transport theorem

A fluid is composed of a huge number of fluid molecules. From the chemical or atomically physical point of view, the properties of these molecules may be taken into consideration. However, in fluid mechanics, the fluids are studied by considering a control volume much larger than the fluid molecule, and averaging the behavior of the fluid molecules inside. For example, the pressure of the flow field instead of the collision among fluid molecules is considered. These control volumes are microscopically infinitely large but macroscopically infinitely small. They are valid only when the length scale of physical interest is much larger than the mean free path of fluid molecules, which means the fluid is of continuous distribution. This is one of the fundamental hypotheses in fluid mechanics, called continuum approximation. The continuum approximation is valid when the Knudsen number Kn = ι/L is much smaller than unity, where ι is the mean free path of the molecules and L is the length scale of interest. Generally, Kn 1 is satisfied in the nature, except in very few cases, for example in the quasi-vacuum environment or at the thinner reaches of the tenuous gases.

Using the concepts of continuum and control volume, the governing equations of flow motions can be derived by considering the conservation laws of mass, momentum or energy. The physical laws applied to an control volume are always in Lagrangian description. As aforementioned, the operation of material derivative ( D Dt ) is needed, and it is necessary to transform such operations into the volume integrals of Eulerian description. The theorem for this transformation is called Reynolds' transport theorem. Let V (t) be a control volume of arbitrary shape and F (t) be a quantity of interest as functions of time t. The Reynolds' transport theorem states

D Dt V (t) F (t)dV = S(t) F (t)u • dS + V (t) ∂F ∂t
dV.

(1.

3)

The Reynolds' transport theorem can be mathematically derived by considering the definition of derivative from the perspective of limit [START_REF] Kundu | Fluid Mechanics[END_REF][START_REF] Currie | Fundamental Mechanics of Fluids[END_REF]. The right hand side of Eq. (1.3) are in Eulerian Introduction description and the left hand side ia applicable to refer to specific physical law.

Mass conservation and momentum conservation

By considering the mass conservation and the momentum conservation of fluid motion, the continuity equation and the equation of momentum conservation are readily obtained [START_REF] Kundu | Fluid Mechanics[END_REF][START_REF] Currie | Fundamental Mechanics of Fluids[END_REF]:

1 ρ Dρ Dt + ∇ • u = 0, (1.4) ρ( ∂u ∂t + (u • ∇)u) = ρf + ∇(τ ij ), (1.5)
where ρ is number density of fluid particles, f is the acceleration by the net body force and τ ij is the stress tensor. Mathematically, Eq. (1.4) is derived by considering the F in Eq. ( 1.3) as ρ, and Eq. (1.5) is derived by considering the F as ρu and using Newton's second law.

The set of equations of Eq. (1.4) and Eq. (1.5) has more unknowns than the number of equations. In order to produce a solvable system, a constitutive equation about the stress tensor τ ij (for newtonian fluid) is needed:

τ ij = -pδ ij + 2µ(S ij - 1 3 (∇ • u)δ ij ) + µ v (∇ • u)δ ij , (1.6)
where p is the pressure; µ is the coefficient of dynamic viscosity; µ v is the coefficient of bulk viscosity, which can be considered as 0 by Stokes assumption [START_REF] Kundu | Fluid Mechanics[END_REF]; δ ij is Kronecker delta tensor, whose component is 1 when i = j and 0 when i = j; S ij = 1 2 ( ui xj + uj xi ) is the strain rate tensor. Substituting Eq. (1.6) into Eq. (1.5) produces the Navier-Stokes Equations:

ρ( ∂u ∂t + (u • ∇)u) = ρf -∇p + µ u + (µ v + 1 3 µ)∇ • u. (1.7)
In this thesis, only the incompressible turbulent flows are considered, which are of constant density. For incompressible fluids, Eq. (1.4) and Eq. (1.7) are simplified to:

∇ • u = 0, (1.8) and ∂u ∂t + (u • ∇)u = -∇p/ρ + ν u + f , (1.9)
where ν = µ/ρ is the coefficient of the kinematic viscosity.

1.1 -Turbulent flow 11

Eq. (1.9) can be made dimensionless by choosing characteristic length scale L and characteristic velocity amplitude U as reference scales. This leads to the following dimensionless Navier-Stokes Equations for incompressible fluids:

∂u ∂t + (u • ∇)u = Re -1 u -∇p + f , (1.10)
where

Re = U L ν (1.11)
is the Reynolds number. Turbulent flow is characterized by large Reynolds number, and the flow with small

Reynolds number is called laminar flow. When Re = LU/ν 1, the non-linear convection term ((u • ∇)u) is dominant, and the inertial force of the fluid motion is much stronger than the viscous force which is to slow down the fluid motion. Thus the flow is chaotic and unstable. More importantly, the flow with large Re is of most interest of physicists, because the flow motion is of a large span of length scales. Then the inertial range of length scales, which are between the largest scales and the smallest scales and believed to contain some universal properties, is more distinguishable.

Vortex stretching

The turbulent motion can be considered as a collection of numerous eddies. Thus an important feature of turbulence is the rotational motion [START_REF] Pope | Turbulent Flows[END_REF]. By taking the curl of Eq. (1.10) without external body forcing, it yields

Dω Dt = Re -1 ω + ω • ∇u, (1.12)
where ω = ∇ × u is called as the vector of vorticity, which is twice of the rotation rate of the fluid. In the circumstance of large Re, the viscous term Re -1 ω can be neglected.

Considering an infinitesimally small material line l = δx, it evolves as

Dl Dt = l • ∇u, (1.13)
which is exactly the same as Eq. (1.12) without the viscous term. Thus the vorticity component in the direction, along which the fluid material line element is stretched by the velocity gradient, tends to increase. In other word, the vector of the vorticity tends to align with the principle direction of the mean stain rate. Such a phenomenon is called as vortex stretching. It is very important in turbulent motion because it is the process about how the eddies extract energy from the mean shear flow. It is noteworthy that the process of vortex Introduction stretching exists not only in two dimensional turbulence since the angular momentum of eddy is concerned.

This indicates the three dimensional and two dimensional turbulence are fundamentally different.

Energy cascade

In about 500 years ago, the turbulent flow was recognized to be composed by eddies of different sizes, e. The energy from outside of the turbulent flow is injected, by external forces or boundary conditions, into the flow at the large scales. The flow motion at these scales are of largest Reynolds number, and most unstable. As a result, they eventually will break up into smaller and smaller eddies. Such a process stops when the eddies of small enough scales are reached, in which the viscosity is important. The breaking-up of eddies is accompanied with the transfer of energy, which will eventually be dissipated into heat by viscosity. This whole process is called turbulent cascade. An illustration about this process is given in figure 1.2. The concept of energy cascade in turbulence was first introduced by [START_REF] Richardson | Weather Prediction by Numerical Process[END_REF] with a famous verse:

Big whirls have little whirls that feed on their velocity, And little whirls have lesser whirls and so on to viscosity.

In Navier-Stokes equations (Eq. (1.10)), the convection term ((u • ∇)u) is responsible for transferring the energy to smaller scales and the dissipation term (Re -1 u) is responsible for dissipating the energy into heat at the smallest scales. In order to have stationary turbulence, energy have to be added at the large scales to maintain the motion of eddies, otherwise the entire kinetic energy will gradually vanish, through the cascade process.

1.1 -Turbulent flow 13 turbulent motions of eddies much smaller than energy-containing or forcing scales, are statistically isotropic, in other words are independent of the motions of large eddies. In reality the motions of large scale eddies (L), can be inhomogeneous and anisotropic but in the cascade the information is lost so that the motion of very small scale can be locally homogeneous and isotropic.

Eulerian point of view 29

! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! "##$!%&'(! )*+(!,-./(0! " , ( 0 1 $ !2 * % 2 * # ( ! 3&%24-%!#&%%&5*6&4,! ",(01$!&,7(26&4,! 8!
The second hypothesis of Kolmogorov was that at very high Re number the statistics of these small scale turbulent motions are independent from large scales and characterized by the kinematic viscosity (n) and the average rate of dissipation of turbulence kinetic energy per unit mass (e). By

Energy dissipation, e

The average rate of dissipation of turbulence kinetic energy per unit mass means of dimensionless analysis, Kolmogorov scales were introduced as follows:

h = ✓ n 3 e ◆ 1/4 (3.11) u h = (ne) 1/4
(3.12) [START_REF] Ardeshiri | Dynamics of Copepods in Turbulent Flows[END_REF]. L is the length scale of domain and η is the smallest length scale, kolmogorov microscale. Energy is injected from outside of the flow at the largest scales.

t h = ⇣ n e ⌘ 1/2 (3.13) Figure 1.2: Illustration of turbulent cascade
Then the largest eddies break up into eddies of smaller and smaller size and energy is meanwhile transferred. This process continues until the smallest eddies and the energy is dissipated by viscosity.

Kolmogorov's 1941 theory

In general, the N-S equations can not be solved analytically. However based on the N-S equations, physicists have obtained some theories describing the fundamentals of turbulent motions, such as those proposed by [START_REF] Richardson | Weather Prediction by Numerical Process[END_REF], [START_REF] Taylor | Statistical theory of turbulence: Parts i-ii[END_REF] and Kolmogorov (1941b,a). A milestone in the studies of turbulence is the Kolmogorov theory in 1941 (Kolmogorov, 1941a,b), which is the first quantitative description to the energy cascade process in turbulence. The Kolmogorov's 1941 (K41) theory is based on the following assumpotions:

Kolmogorov's hypothesis of local isotropy. At sufficiently high Reynolds number, the smallscale turbulent motions of the scale (l L) are statistically isotropic.

Kolmogorov's first similarity hypothesis. In every turbulent flow at sufficiently high Reynolds number, the statistics of the small-scale turbulent motions of the scale (l L) have a universal form that is uniquely determined by ν and (energy dissipation rate, the average rate of dissipation of turbulent kinetic energy per unit mass).

Kolmogorov's second similarity hypothesis. In every turbulent flow at sufficiently high Reynolds number, the statistics of the small-scale turbulent motions of the scale in the range (η l L) have a universal form that is uniquely determined by , independent of ν .

The hypothesis of local isotropy states that the turbulent motions of small scales are isotropic, although the forcing or boundary conditions at large scale are generally anisotropic. When the energy cascade reaches small Introduction enough scales, the information of largest scales is lost and the isotropy is restored.

The first similarity hypothesis states that when the information of largest scales is lost, the statistical properties of turbulent flows are fully dependent on the viscosity ν and dissipation rate ( ν 2 ( ∂ui ∂xj + ∂uj ∂xi ) 2 ). Therefore, by dimensional analysis, the smallest scales (Kolmogorov scales) are estimated as:

η = ν 3 1/4 , u η = (ν ) 1/4 , τ η = η u η = ν 1/2 , (1.14)
where η, u η and τ η are the Kolmogorov length scale, Kolmogorov velocity scale and Kolmogorov time scale respectively. The Reynolds number at the Kolmogorov scale is Re η = (ηu η )/ν = 1.

The second similarity hypothesis is about the inertial range of scale, where the effect of energy dissipation is also negligible. Thus the dissipation rate , which is also the energy injection rate at largest scales and energy transfer rate at the intermediate scales because of energy conservation, is the only key parameter to determine the properties of flow motion.

The motion of the largest eddies is of the integral length scale L and characteristic velocity U . By dimensional analysis, L can be estimated as .15) According to Eq. (1.11) and (1.14), the ratio of the largest scale to the smallest scale is .16) This explains why the flow motion with larger Reynolds number spans a larger range of scales.

L = U 3 . ( 1
L η = Re 3/4 . ( 1
Kolmogorov assumes that the flow is self-similar at the scale much smaller than the largest scales. An important consequence of this assumption is that there is a unique scaling exponent h in the scaling relation of .17) This scaling relation implies the scaling behavior of the structure functions

the velocity difference δu(l) = u(x + l) -u(x): δu(αl) = α h δu(l). ( 1
S q (l) = |δu| q = |u(x + l) -u(x)| q ∼ l ζ(q) , (1.18)
where ζ(q) is the scaling exponent of qth order structure function.

-Turbulent flow 15

By considering the famous Kármán-Howarth equation in terms of the structure function and neglecting the viscous term in the inertial range [START_REF] Monin | Statistical Fluid Mechanics: Mechanics of Turbulence[END_REF][START_REF] Pope | Turbulent Flows[END_REF], the Kolmogorov 4/5th law can be obtained as

(u(x + l) -u(x)) 3 = - 4 5 l. (1.19)
This is an exact result about the scaling behavior of third order longitudinal structure function derived from the Navier-Stokes equations, indicating the scaling exponent of third order structure function is

ζ(3) = 1. K41
theory suggests that ζ(q) is linearly dependant on q, i.e. ζ(q) = q/3.

However, the linear dependance between ζ(q) and q was later found incorrect, because of the internal intermittency. The phenomenon of the internal intermittency in turbulence was first found in the turbulent motion by [START_REF] Batchelor | The nature of turbulent motion at large wave-numbers[END_REF], when they revealed that the instantaneous energy dissipation rate could evolve to large value intermittently in their experiment. In order to take the internal intermittency into account, the refined similarity hypotheses were proposed by [START_REF] Kolmogorov | A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high reynolds number[END_REF] and [START_REF] Obukhov | Some specific features of atmospheric turbulence[END_REF] for the velocity field. This hypothesis leads to a non-linear expression of ζ(q) as

ζ(q) = q 3 - µ 18 (q 2 -3q), (1.20)
where µ is a constant, called as intermittency exponent. Later, in the the multifractal framework [START_REF] Yaglom | The influence on the fluctuation in energy dissipation on the shape of turbulent characteristics in the inertial interval[END_REF][START_REF] Mandelbrot | Intermittent turbulence in self-similar cascades: divergence of high moments and dimension of the carrier[END_REF][START_REF] Frisch | Turbulence[END_REF], which suggests that the energy cascade has multifractal statistics, the scaling exponent is estimated as

ζ(q) = q 3 -K( q 3
), (1.21) where K( q 3 ) is model-dependant (for example β model [START_REF] Frisch | A simple dynamical model of intermittent fully developed turbulence[END_REF] or lognormal model [START_REF] Obukhov | Some specific features of atmospheric turbulence[END_REF][START_REF] Kolmogorov | A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high reynolds number[END_REF]).

Although K41 theory was found questionable, for second order structure function, ζ(2) = 2/3 is nevertheless a good estimation. In the work by Kolmogorov (1941b), the second order structure function was considered.

And by dimensional analysis, it was derived that

S 2 (l) = C 2/3 l 2/3 , (1.22)
where C is the Kolmogorov constant and believed to be universal. Such a result confirms the estimation of

ζ(2) = 2/3.
The scaling law of the energy cascade can also be derived based on the results of structure functions.

Introduction

Taking the advantage of Fourier Transform, for the eddies of length scale of l, k = 2π/l is the corresponding wave number. To fully understand the energy cascade, the velocity u is decomposed into mean part and fluctuating part: .23) in which • means the ensemble average, u is the mean velocity, i.e. the motion of 0 wave number (largest scale), and u is the sum of the velocity of other modes (smaller scales). An important feature of turbulent flow is that u is comparable with or dominant over u . The tensor of autocovariance is

u = u + u . ( 1
R ij (r) = u i (x)u j (x + r) . (1.24)
Taking the Fourier transform of R ij (r) produces

Φ ij (k) = 1 8π 3 R ij (r)e -ik•r dr, (1.25)
where k is the vector of wave number. Then the energy spectrum is defined in k space, as the integral of 1 2 Φ ii in the sphere with radius of |k|:

E(k) = 1 2 Φ ii (k)δ(k -|k|)dk, (1.26)
where the repeating index of Φ ii is Einstein summation convention, i.e. 1 2 Φ ii = 1 2 (Φ 11 + Φ 22 + Φ 33 ) in three dimensional space.

The autocovariance R ij has the same dimension with the second order structure function S 2 (l). Thus by dimensional analysis, a -5/3 scaling relation can be found in the energy spectrum [START_REF] Obukhov | Spectral energy distribution in a turbulent flow[END_REF]:

E(k) = C 2/3 k -5/3 , (1.27)
where C is a constant. This result is known as Kolomogorov-Obukhov's -5/3 law (Kolmogorov, 1941a,b;[START_REF] Obukhov | Spectral energy distribution in a turbulent flow[END_REF], quantitatively describing the energy cascade in the inertial range.

Passive Scalars

The concentration of a substance or intensity of a property convected by a turbulent flow can exhibit a very complex structure, evolving broadly in both space and time. Such substances can be dye diluted in the water or smoke dispersing in air. The property can be the temperature of air flow when cooling a hot fin. The [START_REF] Monin | Statistical Fluid Mechanics: Mechanics of Turbulence[END_REF][START_REF] Libby | Turbulent flows involving chemical reactions[END_REF][START_REF] Shraiman | Scalar turbulence[END_REF][START_REF] Warhaft | Passive scalars in turbulent flows[END_REF][START_REF] Mitrovic | Effects of a first-order chemical reaction on turbulent mass transfer[END_REF]).

Advection-diffusion equation

Let θ be a passive scalar field, for example representing the molecular concentration of any species. By considering the F in Eq. ( 1.3) as the volumetric density (θρ) and using the Fick's law of diffusion [START_REF] Bergman | Fundamentals of Heat and Mass Transfer[END_REF], the advection-diffusion equation of passive scalar can be obtained. In incompressible fluid, it reads

∂θ ∂t + (u • ∇)θ = D θ + q, (1.28)
where θ is the passive scalar, D is the species diffusivity for species or thermal conductivity for temperature, q

is the external source or sink. In the dimensionless from, the governing equation of passive scalar is .29) where Sc = ν/D is the Schmidt number.

∂θ ∂t + (u • ∇)θ = (Re • Sc) -1 θ + q, ( 1
Eq. (1.29) is in a form similar with Eq. (1.10). Thus physicists are interested in analogizing the properties of passive scalar to that of velocity, especially questing for the universality in the inertial range of scales and under the condition of large enough Reynolds (Re) and Peclét (P e = Re • Sc) numbers.

Introduction

Moments of passive scalar

Similar with the Reynolds decomposition of velocity u (Eq. (1.23)), θ can be decomposed into mean part and fluctuating part:

θ = θ + θ . (1.30)
The θ is the mean quantity of scalar, i.e. the amplitude of the 0 wave number mode (largest scale), and θ is the sum of the rest modes (smaller scales). Then the scalar variance is defined as the mean square of the fluctuating part: θ 2 . The scalar variance describing the magnitude of the fluctuation of θ is also called as scalar energy, because it is an analogy of the kinematic energy for velocity (u i u i /2) . Considering the ensemble average ( • ) of Eq. ( 1.28), it yields .31) The covariance between passive scalar and velocity u θ at the right-hand-side of the Eq. ( 1.31) is the turbulent flux of passive scalar.

∂ θ ∂t + ∇ • ( u θ ) = D θ -∇ • u θ + q . ( 1
The passive scalar flux is generally unknown. Its determination belongs to the closure of turbulence. As an analogy with molecular motion and the introduction of the diffusivity constant, a simplifying assumption is sometimes used, which sates that the passive scalar flux is linked with the mean gradient, using a diffusivity tensor. This was proposed by [START_REF] Batchelor | Diffusion in a field of homogeneous turbulence. i. eulerian analysis[END_REF], as

u θ = -D T,ij ∇ θ , (1.32)
where D T,ij is the turbulent diffusivity tensor. A further simplification is to assume that this diffusivity tensor is diagonal [START_REF] Monin | Statistical Fluid Mechanics: Mechanics of Turbulence[END_REF]. A last simplification is to assume that this tensor is proportional to the unit tensor, and thus depends on only one parameter D T , providing a proportionality between vectors [START_REF] Pope | Turbulent Flows[END_REF][START_REF] Tennekes | A First Course in Turbulence[END_REF]:

u θ = -D T ∇ θ , (1.33)
where D T is the turbulent diffusivity of the passive scalar. This is called the mixing length hypothesis, which is useful for modelling purpose, but is not compatible with intermittency and long-range correlations in turbulent fluctuations. As mentioned in the book by [START_REF] Tennekes | A First Course in Turbulence[END_REF], the mixing length hypothesis is rather well satisfied only when the turbulent flow is mean gradient driven. But in other cases, such as buoyancy driven flow or homogeneous flow, the mixing length hypothesis is not valid.

-Passive Scalars 19

The energy of the passive scalar is defined as the variance of the fluctuation θ 2 . In a canonical case where the mean gradient of passive scalar ∇θ is constant and the external source is of zero mean q = 0, the governing equation of θ can be derived from the Reynolds decomposition of Eq. (1.28), as

Dθ Dt = D θ + q -u • ∇ θ . (1.34)
By multiplying Eq. ( 1.34) with θ and taking ensemble average ( • ), it can be obtained that

1 2 D θ 2 Dt = -D (∇θ ) 2 + θ q -u θ • ∇ θ . (1.35)
The minus sign of the dissipation term D |∇θ 2 | indicates it is responsible for the energy lose of the fluctuation of passive scalar. On the other hand, in addition to the energy injection from the external source, the fluctuation of passive scalar can also be sustained by the passive scalar mean gradient, whose energy injection power is proportional to the flux of passive scalar u θ .

Obukhov-Corrsin theory

Similarily with Kolmogorov's and Obukhov's theories for the scaling properties of the velocity field, results for passive scalars can also be obtained.

The autocovariance of scalar is defined as

R θ (r) = θ (x)θ (x + r) . (1.36)
Taking the Fourier transform of R θ (r) produces

Φ θ (k) = 1 8π 3 R θ (r)e -ik•r dr. (1.37)
Then the energy spectrum is defined in k space, as the integral of 1 2 Φ θ in the sphere with radius of |k|:

E θ (k) = 1 2 Φ θ (k)δ(k -|k|)dk. (1.38)
The global scalar variance is the integral of E θ (k) all over the amplitude of wave numbers:

θ 2 = E θ (k)dk.
(1.39)

Introduction

Here the the energy spectrum of the scalar E θ (k) describes the fluctuation of θ at the length scale of (2π/k).

Following Kolmogorov (1941b), [START_REF] Obukhov | Structure of temperature field in turbulent flow[END_REF] and [START_REF] Corrsin | On the spectrum of isotropic temperature fluctuations in an isotropic turbulence[END_REF] argued that, under the condition of high enough Re and P e numbers, there is also an energy cascade process of passive scalars from large to small scales, at which the fields of passive scalars are locally isotropic. Together with the scaling of the velocity field, this forms the classical KOC (Kolmogorov-Obukhov-Corrsin) theoretical framework.

In real space, Obukhov and Corrsin's results [START_REF] Obukhov | Structure of temperature field in turbulent flow[END_REF][START_REF] Corrsin | On the spectrum of isotropic temperature fluctuations in an isotropic turbulence[END_REF], stated independently, are that the scaling relation second order structure function of a passive scalar follow: .40) and in Fourier space:

S 2 θ (l) = |δθ| 2 = C θ -1/3 θ l 2/3 , ( 1
E θ (k) = C θ -1/3 θ k -5/3 , (1.41)
where C θ and C θ are constants and θ (2D ∂θ ∂xi ∂θ ∂xi ) is the dissipation rate of scalar variance. The energy spectrum of scalar variance obeys a similar scaling law with the kinematic energy spectrum in the inertial range.

Effect of Schmidt number

An extension of Obukhov-Corrsin theory was done by [START_REF] Batchelor | Small-scale variation of convected quantities like temperature in turbulent fluid. part 1. general discussion and the case of small conductivity[END_REF][START_REF] Batchelor | Small-scale variation of convected quantities like temperature in turbulent fluid. part 1. general discussion and the case of small conductivity[END_REF][START_REF] Batchelor | Small-scale variation of convected quantities like temperature in turbulent fluid. part 2. the case of large conductivity[END_REF].

In the study about turbulent passive scalar field, the smallest length scale is characterized by the Batchelor scale η B [START_REF] Batchelor | Small-scale variation of convected quantities like temperature in turbulent fluid. part 1. general discussion and the case of small conductivity[END_REF] which is different from the smallest length scale of velocity field (Kolmogorov scale η) with a factor of the square root of Schmidt number Sc (ν/D), as:

η B = η Sc 1/2 .
(1.42)

When Sc ∼ 1, the -5/3 law for passive scalar refers to the scaling behavior in the inertial range (between large scale L and η). For the case of very large Sc (η η B ), the length scales between L and η is referred as inertial-convective range, in which the -5/3 scaling law holds. Moreover, in this circumstance, the range of length scales even smaller than η but larger than η B (called as the viscous-convective range) is also of interest, in which the energy of passive scalar variation scales as k -1 [START_REF] Batchelor | Small-scale variation of convected quantities like temperature in turbulent fluid. part 1. general discussion and the case of small conductivity[END_REF]. For a passive scalar field with Sc 1 (η B η), the inertial range can be divided into inertial-convective range (between L and η B ) and inertial-diffusive range (between η B and η). The -5/3 law holds in the inertial-convective range, while the 1.2 -Passive Scalars 21 scaling coefficient in the inertial-diffusive range is -17/3 [START_REF] Batchelor | Small-scale variation of convected quantities like temperature in turbulent fluid. part 2. the case of large conductivity[END_REF].

Numerical studies of the Obukhov-Corrsin theory

Soon after the use of computers, a lot of numerical studies based on the N-S equations and the convectiondiffusion equation of passive scalar were implemented for examining the Obukhov-Corrsin theory [START_REF] Sreenivasan | On local isotropy of passive scalars in turbulent shear flows[END_REF][START_REF] Sreenivasan | The phenomenology of small-scale turbulence[END_REF], whose prerequisite is high enough Re and P e numbers. With the development of super computers, the numerical spatial resolution becomes finer and finer. Thus the direct numerical simulations about the passive scalar in turbulence with high Re and P e numbers (the Kolmogorov scale is smaller and finer grid size is required to resolve it) is possible.

The -5/3 scaling behavior in the energy spectrum of the passive scalar variance is one main conclusion of the Obukhov-Corrsin theory. This have been well observed in some numerical studies [START_REF] Warhaft | Passive scalars in turbulent flows[END_REF][START_REF] Yeung | High-reynolds-number simulation of turbulent mixing[END_REF]. In the past few decades, several numerical studies about the passive scalar spectrum were interested in the effect of Sc number. Some works focused on the cases of small Sc [START_REF] Yeung | Schmidt number effects on turbulent transport with uniform mean scalar gradient[END_REF][START_REF] Yeung | Direct numerical simulation of turbulent mixing at very low schmidt number with a uniform mean gradient[END_REF][START_REF] Hill | Spectra of turbulently advected scalars that have small schmidt number[END_REF]. In the work by [START_REF] Briard | Passive scalar convective-diffusive subrange for low prandtl numbers in isotropic turbulence[END_REF], the energy spectrum of passive scalar variance in isotropic turbulence, with Re = 160 and P r between 10 -6 and 1, was numerically studied * . For the case of P r = 1, a -5/3 slope was found for the spectrum between P r 3/4 k η and k η , where k η is wave number corresponding to the Kolmogorov scale. With the decrease of P r, a more detailed scaling behavior than that stated by [START_REF] Batchelor | Small-scale variation of convected quantities like temperature in turbulent fluid. part 2. the case of large conductivity[END_REF] was detected. The spectrum gradually evolves into a combination of a -17/3 slope between P r 3/4 k η and √ P rk η and a -11/3 slope between √ P rk η and k η . Meanwhile, high Sc number, for which even better resolution is required to resolve the Batchelor scale (η/Sc 1/2 ), is also of high interest [START_REF] Warhaft | Passive scalars in turbulent flows[END_REF]Gotoh et al., 2015). In the work by [START_REF] Donzis | The batchelor spectrum for mixing of passive scalars in isotropic turbulence[END_REF], direct numerical simulation data was used to examine the spectrum of passive scalar fluctuation maintained by uniform mean gradient in forced stationary isotropic turbulence. The Sc number varies from 1/8 to 1024, and Taylor-scale Reynolds number (Re λ ) varies from 8 to 650 on the periodic domains of 64 3 to 4096 3 grid points. It was clearly observed that, for all Reynolds numbers, there was a general trend for the high wave number spectrum to flatten towards k -1

as Sc increases.

Probability density function, structure function and intermittency

The probability density function (PDF) of the passive scalar did not received attention as much as the energy spectrum until the 1980s, because the fluctuation of passive scalar normalized by its standard derivation had [START_REF] Jayesh | Probability distribution of a passive scalar in grid-generated turbulence[END_REF]). For legibility purpose, the upper curve has been shifted by 2 decades with respect to the lower one. Note that in the labels, the θ and u are the fluctuations and θ and u are the corresponding rms. The deviation of PDF from Gaussian curve for the scalar is larger than that of the velocity, because the intermittency of the scalar field is relatively stronger.

been assumed of Gaussian distribution in homogeneous turbulence, which was a satisfactory model. However, it was then found in experimental and theoretical studies that the PDF of the passive scalar has exponential tails (see figure 1.3). It indicates that the probability of the scalar to reach extremely high or low values is larger than what was believed before. Such a phenomenon is of practical interest, for example when the scalar is the concentration of a toxic chemical. The exponential tails in the profile of scalar PDF is a manifestation of the internal intermittency of passive scalars, associated with the sharp ramp-cliff structures in the small scale scalar field [START_REF] Shraiman | Scalar turbulence[END_REF][START_REF] Warhaft | Passive scalars in turbulent flows[END_REF]. The most common measure of intermittency is the scaling exponents of the structure functions, hence not for the field itself, but for its increments. For passive scalars, similarily with Eq. (1.18), the KOC theory implies a scaling relation of the structure function of the passive scalar as

S q θ (l) = |δθ| q = |θ(x + l) -θ(x)| q ∼ l ζ θ (q) , (1.43)
where ζ θ (q) is the scaling exponent of qth order passive scalar structure function. Without intermittency, the scaling behaviors of passive scalar and velocity are supposed to be identical, as

ζ(q) = ζ θ (q) = q/3. (1.44)
Taking the internal intermittency into consideration, similar with velocity [START_REF] Kolmogorov | A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high reynolds number[END_REF][START_REF] Obukhov | Some specific features of atmospheric turbulence[END_REF], the refined similarity hypotheses for passive scalars were also proposed [START_REF] Stolovitzky | Refined similarity hypotheses for passive scalars mixed by turbulence[END_REF]. Then it was soon confirmed that, with intermittency, the ζ-q and ζ θ -q curves are nonlinear instead of straight lines (see figure 1.5).

With the help of super computers, more precise inspection of the intermittency is an interesting topic for the numerical simulations of passive scalar in high Reynolds number turbulent flow [START_REF] Warhaft | Passive scalars in turbulent flows[END_REF][START_REF] Borgas | High schmidt number scalars in turbulence: Structure functions and lagrangian theory[END_REF][START_REF] Yeung | High-reynolds-number simulation of turbulent mixing[END_REF]. [START_REF] Iyer | Steep cliffs and saturated exponents in three-dimensional scalar turbulence[END_REF] studied the intermittency of a passive scalar in three-dimensional turbulence at Taylor-scale Reynolds number of 650, using direct numerical simulations on a 4096 3 grid. The Schmidt number is unity. The passive scalar fluctuation was forced by a uniform mean gradient. The scaling exponents of structure function of order higher than about 12 was found to saturate to about 1.2. Interestingly, the fractal dimension of the steepest fronts, associated with very large scalar gradients, was found as about 1.8, which adds up with the saturation value of scaling exponent to be the number of spatial dimension 3. In the work by Gotoh and Watanabe (2015), direct numerical simulation with 4096 [START_REF] Schmitt | Relating lagrangian passive scalar scaling exponents to eulerian scaling exponents in turbulence[END_REF][START_REF] Schmitt | Linking eulerian and lagrangian structure functions' scaling exponents in turbulence[END_REF], compared with the theoretical predictions of K41 and KOC theories [START_REF] Huang | Arbitrary Order Hilbert Spectral Analysis Definition and Application to Fully Developed Turbulence and Environmental Time Series[END_REF].

Chapter 2

Reactive scalars in incompressible turbulence

The problem of scalar transport and mixing becomes physically rich when reactions are taken into consideration, such as biological or chemical reactions. For example, in the marine system, many micro-biological species are under the influence of both the turbulent dynamics of their living fluid environment and the reactive interactions between the species, or species and nutrition distribution. Specifically, species interaction can be modeled as chemical or biological reactions, e.g. preys and predators or species and nutrition, involving chemical or biological quantities (such as the concentration of phytho-or zooplanktons) as reactive scalars. Different from the passive scalars, which are transported without reaction, reactive scalars are more complex because the reactions almost always make the system non-linear (Eq. (2.7)). Because of the motion of turbulent eddies in the ocean, the plankton and nutrition as reactants and products of reactions are efficiently mixed and redistributed. When transported in the turbulent flow, the turbulence is responsible for not only bringing reactants together so that reaction can occur within a finite time, but also for dispersing the products of reaction. These species and nutrients are at the bottom of the food chain in the marine world. Thus the dynamics of such system is important in a number of aspects, e.g. larval recruitment, fishing and environmental monitoring [START_REF] Hill | Homogeneous turbulent mixing with chemical reaction[END_REF][START_REF] Warhaft | Passive scalars in turbulent flows[END_REF][START_REF] Neufeld | Chemical and Biological Processes in Fluid Flows: A Dynamical Systems Approach[END_REF].

From a fundamental point of view, such a reactive scalar mixing problem in incompressible turbulence is especially meaningful in the following senses. Firstly, the nonreactive passive scalar turbulence alone has been the subject of many studies [START_REF] Monin | Statistical Fluid Mechanics: Mechanics of Turbulence[END_REF][START_REF] Libby | Turbulent flows involving chemical reactions[END_REF][START_REF] Mitrovic | Effects of a first-order chemical reaction on turbulent mass transfer[END_REF]. When the chemical source is added, how scalars interact needs to be taken into consideration. Secondly, reactive turbulence has been typically considered as the combustion case, with large variation of the fluid density and very fast reaction rates. There exist abundant literature on turbulence combustion [START_REF] Peters | Turbulent Combustion[END_REF][START_REF] Pope | Turbulent Flows[END_REF][START_REF] Poinsot | Prediction and control of combustion instabilities in real engines[END_REF]Zhao et al., 2018a;[START_REF] Zhang | Non-premixed turbulent combustion modeling based on the filtered turbulent flamelet equation[END_REF]. Relatively, the theoretical knowledge of statistical laws, including scaling behavior, of the reactive scalars and the effects of chemical reaction on the passive scalar mixing in incompressible turbulence is still very sparse.

In this chapter the basic concepts, fundamental theories and some numerical studies about reactive scalars in incompressible turbulence are introduced. It is organized as follows. Section 2.1 introduces the basic chemical kinetics for formulating the reaction source terms in the governing equations of reactive scalars, and some chemical and biological models, which simplify the real reacting system from a theoretical point of view. In section 2.2, the statistical properties of reactive scalars, including the mean, variance and correlation are discussed based on the governing equations of reactive scalars. Last but not least, in section 2.3, a literature review about the analytical, numerical and experimental studies of turbulent reactive scalar is presented.

Chemical kinetics

General descriptions

The basic concepts about chemical kinetics [START_REF] Law | Combustion Physics[END_REF][START_REF] Neufeld | Chemical and Biological Processes in Fluid Flows: A Dynamical Systems Approach[END_REF] are introduced here, before being used in the fluid mechanics framework.

The general form of irreversible reaction between N species:

M 1 , M 2 , . . . , M N is N i=1 n i M i γ 1 G GGGGG A N i=1 n i M i , (2.1)
where M i can be a reactant or a product; γ 1 is the reaction rate coefficient; n i and n i are the reaction order of M i . The order of the reaction is

N i=1 n i . Considering the mole concentration of M i : θ i , the rate of change of θ i is ω i = dθ i dt . (2.2)
For any two reacting species M i and M j , a species independent quantity ω can be introduced:

ω = ω i n i -n i = ω j n j -n j . (2.3)
Thus ω is defined as the reaction rate. The phenomenological law of chemical reaction states that ω is propor-2.1 -Chemical kinetics 27 tional to the product of the concentrations of the reactants:

ω = γ 1 N i=1 θ n i i . (2.4)
It is noteworthy that, according to Eq. ( 2.4), the dimension of the reaction rate coefficient γ 1 is not universal, but dependent on which reaction it refers to.

In reality, there is no absolutely irreversible reaction, every reaction of (2.1) is associated with the reverse reaction of it. The overall reaction is a reversible reaction:

N i=1 n i M i γ 1 G GGGGG B F GGGGG G γ 2 N i=1 n i M i , (2.5)
where γ 2 is the reaction rate coefficient of backward reaction. The reaction rate of the backward reaction is defined similar with the forward reaction. Thus the net reaction is

ω f orward -ω backward = γ 1 N i=1 θ n i i -γ 2 N i=1 θ n i i . (2.6)
For specific species M j , the net chemical source is

ω j = (n j -n j )(γ 1 N i=1 θ n i i -γ 2 N i=1 θ n i i ).
(2.7)

After enough long time and without external supply of reactant, the reversible reaction will finally reach its chemical equilibrium state, at which both the forward and backward reactions are still ongoing but of the same rate, i.e. the net reaction rate is 0.

Models in chemical and biological systems

In a incompressible turbulence environment, such as in the ocean, various reactions exist. For example the decomposition of carbonic acid into H 2 O and CO 2 and the process of photosynthesis. However, chemical reactions in reality are always highly multi-steps and of complex chemical dynamics. From a theoretical point of view, it is appreciated to propose some chemical models summarizing the main mechanism and behaviors in complex chemical reactions. More generally, some biological processes can also be modeled as reaction. For example, the cluster of plankton blooms where the reactants are various nutrients transported by the flow.

Theoretically, the interaction between plankton and their nutrition can be considered as two reactants in one

Reactive scalars in incompressible turbulence reaction with one's number grows at the cost of the consumption of the other. In this section, some mathematical models about the chemical and biological process are reviewed [START_REF] Neufeld | Chemical and Biological Processes in Fluid Flows: A Dynamical Systems Approach[END_REF].

Chemical models

As an example of the chemical kinematics described in section 2.1, an elementary second order irreversible reaction is in the from of:

A + B γ -→ C. (2.8)
If we consider it as a 0-dimensional problem, the concentrations of A, B and C evolve according to

dA dt = dB dt = - dC dt = -γAB, (2.9)
from which the difference between the concentration of A and B is found as invariant, because

d(A -B) dt = 0.
Thus we can define constant Q as

Q = A(t) -B(t) = A(0) -B(0),
which then produces the temporal ordinary partial equation of A as

dA dt = γQA -γA 2 .
(2.10)

Eq. ( 2.10) is an example of how a nonlinear term arise in a single chemical reaction.

In the reaction of (2.8), if the concentration of one reactant is hold as constant (for example the reactant of B), the equation about the temporal evolutions of A and C become (2.11) where γ = γB can be taken as another reaction rate coefficient. Then the reaction is which is called as zeroth reaction.

dA dt = - dC dt = -γ A,
A γ -→ C, ( 2

Biological models

When dealing with a single species, the evolution of the population can be modelled in the framework of birth- Another model for biological reaction is the prey-predator model, which can be used for studying the evolution of fish population due to fishing behavior. Let A and B represent the concentrations of prey and

Reactive scalars in incompressible turbulence predator. In addition to the aforementioned birth and death for both prey and predator, the number of preys is supposed to decrease with large number of predators, and the predators can flourish with plenty of preys. The behavior of preying is the key in the prey-predator model. Thus for the prey, it is assumed that the decrease in population is mostly due to being preyed, i.e. the natural death of prey is less important. Similarly, the increase of predator population is assumed as strongly dependant on the feeding and relatively weakly on proliferation.

Therefore, the evolutions of A and B follow:

dA dt = (γ 1 -γ 2 B)A, dB dt = (γ 2 A -γ 3 )B.
(2.17)

And the corresponding reactions are

A γ1 -→ 2A, A + B γ2 -→ 2B, B γ3 -→ ∅. (2.18)
Consequently, such a chemical scheme produces persistent chemical oscillations, in which A and B increase and decrease alternatively and periodically.

Advection-diffusion-reaction equation

The governing equation of a reactive scalar θ j is the advection-diffusion equation of passive scalar (Eq. (1.28))

with the addition of the chemical source term ω j :

∂θ j ∂t + (u • ∇)θ j = D θ j + ω j + q. (2.19)
Let Θ j (can be the global average of θ j ) be the characteristic quantity of θ j . The characteristic time scale for the reaction is Θ j /ω j . Let us note τ mix the characteristic time scale of scalar transport, which is generally dependant on the scalar convection and diffusion. The Damköhler number Da is defined as the ratio between the time scales of scalar transport and reaction, as:

Da = τ mix Θ j /ω j . (2.20)
Such a dimensionless parameter is a key parameter, quantifying the effect of reaction on the evolution of reactive scalar in comparison with the effect of turbulent convection-diffusion.

In the advection-diffusion-reaction equation (Eq. (2.19)), the reaction term ω j is generally non-linear. As 2.2 -Advection-diffusion-reaction equation a result, it is generally impossible to derive the analytical expression about the statistical properties of the reactive scalar (the mean quantity or the variance etc.), unless the reaction is simple and of low order, such as the first-order reaction [START_REF] Corrsin | Statistical behavior of a reacting mixture in isotropic turbulence[END_REF]. Consider θ as a reactive scalar undergoing first order reaction of

θ γ -→ ∅. (2.21)
in homogeneous isotropic turbulent flow. Then the evolution equation for θ(x, t) reads

∂θ ∂t + (u • ∇)θ = D θ -γθ + q. (2.22)
By writing the left-hand-side of Eq. ( 2.22) in the form of material derivative and doing the ensemble average, it yields

D θ Dt = D θ -γ θ + q . (2.23)
The diffusion term D θ is 0 because of the homogeneity. Thus if there is no external supply ( q = 0), the mean quantity θ is supposed to decay exponentially.

Considering the variance of θ, by multiplying the equation of fluctuating scalar (obtained from (2.22)-(2.23))

and doing ensemble average, it yields

1 2 D θ 2 Dt = -D |∇θ| 2 -γ 1 2 θ 2 + q θ . (2.24)
Similarly with Eq. (1.35), the fluctuations of reactive scalars in homogeneous isotropic turbulence are sustained by the external energy source q θ , and dissipated by the dissipation term -D |∇θ| 2 . In addition, the decrease of scalar fluctuation energy due to reaction is more rapid because the term of -γ 1 2 θ 2 produces an exponential decay.

Consider a second order reaction in the most simple form :

θ 1 + θ 2 γ -→ ∅.
(2.25)

With isotropy and homogeneity, similar with Eq. ( 2.24), the evolution equations of the mean and variance of

θ 1 are D θ 1 Dt = D θ 1 -γ θ 1 θ 2 -γ θ 1 θ 2 + q1 (2.26) and 1 2 D θ 2 1 Dt = -D |∇θ 1 | 2 -γ θ 2 1 θ 2 -γ θ 2 1 θ 2 -γ θ 1 θ 2 θ 1 + q 1 θ 1 . (2.27)
Compared with first order reaction (Eq. ( 2.23) and (2.24)), the evolutions of the mean and variance of reactive scalars in the second reaction are much more complicated, because the non-linear reaction term involves the second and even third order covariance of the fluctuations of reactants. The covariance normalized by the product of the r.m.s of scalar fluctuations is the correlation coefficient between the reactive scalars:

r(θ 1 , θ 2 ) = θ 1 θ 2 θ 2 1 1/2 θ 2 2 1/2 .
(2.28)

The correlation coefficient is a classical mathematical quantity characterizing how strong two zero-mean vari-

-1 ≤ r(θ 1 , θ 2 ) ≤ 1.
(2.29) r(θ 1 , θ 2 ) = 0 means θ 1 and θ 2 are not correlated. For example, in figure 2.1(b), a distribution of θ 1 is independent with θ 2 . A positive r(θ 1 , θ 2 ) corresponds to the case of θ 1 and θ 2 being positively correlated, in which θ 1 is more likely to be the same sign with θ 2 . Figure 2.1(c) shows the joint distribution of θ 1 and θ 2 in the case of r(θ 1 , θ 2 ) = 0.9, which means θ 1 and θ 2 are strongly correlated. Similarly, a negative r(θ 1 , θ 2 ) indicates that θ 1 is more likely to be the opposite sign with θ 2 (figure 2.1(a)). Extremely, when r(θ 1 , θ 2 ) is 1 or -1, θ 1 and θ 2 are perfectly positively or negatively correlated.

2.3 -Interactions between turbulence and reactions 33

Interactions between turbulence and reactions

Taking advantage of computer science and the governing equations, there were some analytical and numerical studies aiming at better understanding and exploring the statistical properties of the reactive scalars in the interaction with incompressible turbulence. On the other hand, a well known example for the turbulent reactive scalars in the real world is the ocean. Thus some experimental studies carried out measurement of the statistical properties of reactive scalars in the ocean or aquatic environment. In this section, some of these studies are reviewed.

Analytical studies on reactive scalars

In fully developed turbulence, universal scaling relations are classically discussed in the framework of Kolmogorov-Obukhov phenomenology for the velocity field, as well as the Obukhov-Corrsin approach for passive scalars advected by the turbulent velocity. This now forms the classical KOC (Kolmogorov-Obukhov-Corrsin) theoretical framework for fluid turbulence with advected passive scalars [START_REF] Warhaft | Passive scalars in turbulent flows[END_REF]. Chemical reactions have been considered quite early in such framework, for some specific cases. The earliest works on the statistical theory of turbulent chemical reactions were done by [START_REF] Corrsin | Statistical behavior of a reacting mixture in isotropic turbulence[END_REF][START_REF] Corrsin | The reactant concentration spectrum in turbulent mixing with a first-order reaction[END_REF]Corrsin ( , 1964a)), [START_REF] Pao | Statistical behavior of a turbulent multicomponent mixture with first-order reactions[END_REF] and [START_REF] O'brien | Decaying second-order isothermal reaction in final period turbulence[END_REF][START_REF] O'brien | Very rapid, isothermal, two-species reactions in final period turbulence[END_REF][START_REF] O'brien | Turbulent diffusion of rapidly reacting chemical species[END_REF].

In the work by [START_REF] Corrsin | Statistical behavior of a reacting mixture in isotropic turbulence[END_REF], the reactive scalar in first order reaction (Γ C -→ ∅) and homogeneous turbulence was considered. The equation of the global scalar variance (Eq. (2.24)), which indicates a exponential decay of scalar variance, was derived. [START_REF] Corrsin | The reactant concentration spectrum in turbulent mixing with a first-order reaction[END_REF] studied the mixing of a scalar contaminant undergoing a first-order chemical reaction in isotropic turbulence. Theoretically he deduced the energy spectrum of the reactive scalar in different wave number ranges. A -5/3 scaling relation was found for the reactive scalar in the inertial range (figure 2.2(a)). Later, [START_REF] Corrsin | Some statistical properties of the product of a turbulent first-order reaction[END_REF] extended the problem to the first order reaction with product (Γ C -→ P), and derived that the mean, variance and energy spectrum of P were similar with that of Γ, only with slight amendment. It is noteworthy that, in the process of deriving the variance of P, the covariance term between Γ and P was involved. [START_REF] Corrsin | Some statistical properties of the product of a turbulent first-order reaction[END_REF] simply assumed that Γ and P are perfectly correlated.

Then Corrsin (1964a) extended the theoretical analysis about the first order reaction to a second order reaction.

A -5/3 scaling relation was also found in the inertial range of the energy spectrum (figure 2.2(b)). Moreover, the nonlinear reaction term causes some additional spectral transfer and more importantly a loss in the spectral energy.

In the work by [START_REF] Pao | Statistical behavior of a turbulent multicomponent mixture with first-order reactions[END_REF], the theories about the evolution of mean and variance of reactive scalar in a first order reaction was developed into a general form, which theoretically quantifies the decaying process of

Having recorded these extreme doubts, but having no clearly better choice, we use (24) to estimate B. (26) is more tedious and not necessarily more reliable. Thus we write

- O2 w B/::: k-8 exp {3Ce-)k-f) dk. (27)
By integration, we get the final result

- (28) 02 ___ Beg [exp {3Ce-+ek*-f} -exp {~C&-~JQ}I. 2 c A log G
Since @k* 9 &L, we neglect the tirst exponential, getting the estj,mate

B NN 2C@e-)exp { -3C~-),kZf}. (29) A FIGURE 1. Qualitative sketch of reactant spectrum in inertial-convective range.
To gain some appreciation of (28) or (29), consider its limiting form for C --f 0. This is obtained most easily by putting (14) into (24), or by expanding (28) as a series, to give

B -+ D M $Pgk%.
Equating this estimate to the better one [equation (17)], we find that

(30) - €0 w O ~J & , E ~. ( 31
)
A similar approach to the inertial turbulent-energy spectrum gives the wellknown estimate for dissipation rate in terms of large-scale parameters,

- E M ( U 2 ) 4 k L , ( 32
)
€8 W e 2 ( U 2 ) + 8 k j J i .

(

) 33
omitting a factor of order unity. Putting

lim t→∞ A 2 ∼ t -3/2 , lim t→∞ A 2 ∼ t -3/2 , lim t→∞ A 2 / A 2 = π -1.
(2.30)

-Interactions between turbulence and reactions 35

While for the non-stoichiometric case, the asymptotical decay of the covariance A B follows t -9/4 e -γt 3/2 , where γ is the reaction rate coefficient.

Numerical studies on reactive scalars

The covariance between reactants and the development of the introduction of models for the covariance terms were important early topics of turbulent mixing analyses [START_REF] Lamb | A model of second-order chemical reactions in turbulent fluid. i. formulation and validation[END_REF][START_REF] Heeb | Turbulent mixing with multiple second-order chemical reactions[END_REF]. [START_REF] Heeb | Turbulent mixing with multiple second-order chemical reactions[END_REF] did numerical simulations about the irreversible reactions during turbulent mixing, and compared with experimental results to examine 14 closure theories. For the two non-premixed reactants case, it was found [START_REF] Toor | Turbulent mixing of two species with and without chemical reactions[END_REF] Among the existing results, the effects of chemical reaction on the turbulent mixing are also of particular interest [START_REF] Hill | Homogeneous turbulent mixing with chemical reaction[END_REF]. For instance, the second-order chemical reactions in a reactor [START_REF] Mao | Second-order chemical reactions with turbulent mixing[END_REF], [START_REF] Chakrabarti | Direct numerical simulation of chemical selectivity in homogenous turbulence[END_REF]. [START_REF] Sykes | Turbulent mixing with chemical reaction in the planetary boundary layer[END_REF] used

or series-parallel reactions (A + B - → R, R + B - → S)
Large Eddy Simulation (LES) to study the chemically reactive mixing of two species, in a system with one reactant injected into an uniformly mixed layer containing the other one. The segregation coefficient, defined as the covariance of the two reactants normalized by the product of their mean quantities, is the measure of the turbulent mixing. It was shown that the two reactants were significantly segregated when the reaction was fast. [START_REF] Komori | The effects of turbulent mixing on the correlation between two species and on concentration fluctuations in non-premixed reacting flows[END_REF] simulated two reactants of second-order irreversible reaction introduced through different parts of the bounding surface of turbulent flow. They developed a model with the Damköhler number based on the integral timescale to estimate the segregation parameter. [START_REF] Leonard | Direct numerical simulation of turbulent flows with chemical reaction[END_REF] studied the decaying reactive scalars undergoing a second-order irreversible reaction in homogeneous turbulence, with pseudospectral method and in up to 128 3 domains. Two interesting results were found. The first one is that the regions in the flow field where reaction rates are the highest are correlated with locations of high strain rates. The other one is that the chemical reaction rate constant does not appreciably affect the microscales of the dissipation of concentration fluctuations, which is supported by a closure theory. Numerical simulations about the decay moments of reactive scalars in irreversible reactions was also carried out, which are one main focus of early theoretical studies [START_REF] Corrsin | Statistical behavior of a reacting mixture in isotropic turbulence[END_REF](Corrsin, , 1964a;;[START_REF] O'brien | Decaying second-order isothermal reaction in final period turbulence[END_REF]. [START_REF] Neufeld | Chaotic mixing induced transitions in reaction-diffusion systems[END_REF] Following these early studies and numerical results, we developed in this thesis numerical simulations of reactive scalars, in different configurations, which are the topic of the different remaining chapters.

Turbulence and reactive scalars in biological oceanography

With ocean turbulence in mind, [START_REF] Lopez | Chaotic advection of reacting substances: Plankton dynamics on a meandering jet[END_REF] numerically studied the nutrient-phytoplankton-zooplankton interactions in chaotic flows. The simulation results displayed a smooth-filamental transition in the concentration patterns. Hernandez-Garcia and [START_REF] Hernandez-Garcia | Sustained plankton blooms under open chaotic flows[END_REF] studied the planktonic population living in an open and chaotic fluid flow, using a predator-prey model. This study showed that a strong chaotic flow is beneficial for sustaining plankton blooms by deforming the filament structures of the flow under the action of stretching and dilution. [START_REF] Groselj | How turbulence regulates biodiversity in systems with cyclic competition[END_REF] carried out a two dimensional numerical study about a cyclic competition between three biological species. The reacting system consists of

A + B - → A, A - → 2A, B + C - → B, B - → 2B, C + A - → C, C - → 2C.
The effect of reaction, in comparison with the turbulent convection-diffusion, on the evolution of the global statistical property of reactants was investigated. A transition from rotating spiral waves to collective oscillations as Damköhler number decreases from very large to very small was observed for the evolutions of the average densities of the species. The work by [START_REF] Powell | Turbulence, diffusion and patchiness in the sea[END_REF] explored the effect the biological interaction on the variance spectrum of the plankton populations in the sea, using a predator-prey model. They concluded that there is no general result in this phenomenon, since the biological interaction may either flatten or steepen the spectrum, which is dependant on the process of diffusion.

The scaling behavior of reactive scalar is also of high interest. Some studies focus on possible applications in the field of chemical and biological oceanography, where the typical times of biogeochemical reactions may be large. [START_REF] Seuront | Multifractal intermittency of eulerian and lagrangian turbulence of ocean temperature and plankton fields[END_REF][START_REF] Seuront | Universal multifractal analysis as a tool to characterize multiscale intermittent patterns: Example of phytoplankton distribution in turbulent coastal waters[END_REF] studied the intermittency of phytoplankton and temperature fields from ocean observations. It has been found that phytoplankton was nearly dynamically passive (similar to temperature) at small scales but biologically active at large scales. The phytoplankton statistics, considered via the proxy of fluorescence measurements, have found some scaling relations with -1.2 spectral slopes, interpreted as signature of biological activity [START_REF] Seuront | Multifractal intermittency of eulerian and lagrangian turbulence of ocean temperature and plankton fields[END_REF][START_REF] Lovejoy | Universal multifractals and ocean patchiness: phytoplankton, physical fields and coastal heterogeneity[END_REF][START_REF] Derot | Long-term high frequency phytoplankton dynamics, recorded from a coastal water autonomous measurement system in the eastern english channel[END_REF]. Such -1.2 slope was not always detected: [START_REF] Yamazaki | Phytoplankton microstructure in fully developed oceanic turbulence[END_REF] measured the microstructure of phytoplankton in fully developed oceanic turbulence by using LED (Light Emitting Diode) fluorescence probe. The spectra of both phytoplankton and velocity exhibited -5/3 slope in the inertial subrange.

An important aspect in understanding ocean turbulence is the species transport process, either in the interior or at the interface of the ocean. This issue is particularly acute for biogeochemical studies. For instance, the mixing process at the ocean and atmosphere interface drives the pelagic food because of the light available for photosynthesis at the surface. Till now, there is no unanimous agreement on the effects of turbulence on species transportation, including the feeding, growth and mortality rates of the species [START_REF] Mackenzie | Turbulence, larval fish ecology and fisheries recruitment: a review of field studies[END_REF]. Under some conditions, the swimming marine species, together with high dissipation rates produced by the aggregations of Reactive scalars in incompressible turbulence living organisms of different sizes, might generate intense turbulence comparable with those by the strong winds near the ocean surface [START_REF] Kunze | Biologically generated mixing in the ocean[END_REF]. In these studies, a major question was how efficient the mixing will be in the ocean. From the important multiscale feature in turbulence, it is reasonable to expect that the mixing efficiency hinges on length scales. Very small whorls introduced into a fluid will be quickly damped by friction, and thus will not mix the fluid. Probably, zooplankton schooling introduces larger scales and increases mixing efficiency [START_REF] Dewar | A fishy mix[END_REF].

So far, there are no considerable number of studies devoted into the turbulence-reaction interaction concerning oceanology and biology. The main reasons are that such a topic is highly interdisciplinary and the conditions in the ocean are not controllable. Thus the problems are too complicated to have universal and unanimous conclusions. In the contrast, more studies about turbulent reactive scalar are interested in the numerical simulations of the canonical cases.

Chapter 3

Numerical implementation of direct numerical simulations

The partial differential equations in turbulence studies (Eqs. (1.9), (1.28), (2.19)) are generally not analytically solvable, i.e. the analytical expressions of the velocity and scalars as functions of continuous spatial and temporal variables (x and t) generally can not be derived from the governing equations. Alternatively, the numerical solutions expressed as functions of discrete variables can be used to describe the flow and scalar fields. These numerical solutions are obtained by solving the set of algebraic equations derived from the governing equations in the discrete domain. To have a detailed and reliable numerical solution, considerable computational time and storage for obtaining and recording the data are needed. Thus the numerical simulations of turbulence problems strongly rely on the availability of computers.

A lot of methods have been developed for the numerical simulations of the turbulence problems. Some numerical methods include models about the turbulent motions, for example the Smagorinsky model [START_REF] Smagorinsky | General circulation experiments with the primitive equations: I. the basic equations[END_REF] and dynamic model [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF][START_REF] Lilly | A proposed modification of the germano subgrid-scale closure method[END_REF][START_REF] Meneveau | A lagrangian dynamic subgrid-scale model of turbulence[END_REF] in Large Eddy Simulation (LES), k -model [START_REF] Jones | The prediction of laminarization with a two-equation model of turbulence[END_REF] and k -ω model [START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF] in Reynolds-Averaged Navier-Stokes (RANS). The general idea of these models is to build the linkage between the small-scale quantities and the large-scale quantities based on the modeling of turbulent fluxes. Then the large-scale quantities are to be numerically solved, while the small-scale statistics are provided by the model. The benefit of using models in numerical simulations is lower computational cost, since the fineness of the spatial resolution is only required for resolving the large scale. However, the drawback is also evident. The models introduce some artificial hypotheses and simplifications to the problem, at the risk of reducing the accuracy and credibility of the results.

In contrast, the Direct Numerical Simulation (DNS) implements numerical simulation based directly on the governing equations without modeling assumptions. The spatial resolution of DNS is required to be fine enough to resolve the smallest scale (e.g. Kolmogorov scale for velocity (1.14) and Batchelor scale for scalar (1.42)).

The results obtained by DNS are much more reliable than using models. Thus DNS is also called numerical experiment, as the preferred numerical method for examining theoretical analysis. However, The computational cost is high [START_REF] Pope | Turbulent Flows[END_REF][START_REF] Laizet | High-order compact schemes for incompressible flows: A simple and efficient method with quasi-spectral accuracy[END_REF][START_REF] Koblitz | Direct numerical simulation of particulate flows with an overset grid method[END_REF], with a computational time roughly proportional to Re 3 [START_REF] Davidson | Turbulence: An Introduction for Scientists and Engineers[END_REF]. The numerical simulations for this thesis are implemented as DNS.

Spatial discretization and approximation of derivatives

The discretization of the spatial domain is the first and very important step for numerical simulations of turbulent flows. It not only defines the framework in which the numerical solutions are expressed, but also allows the approximation of the spatial derivative with an algebraic expression of the discrete quantities. Essentially, the process of spatial discretization means to build a mesh of discrete grid points in the spatial domain. Such a mesh is supposed to be all over the positions of interest, orthogonal or at least near-orthogonal, and fine enough or specially refined to capture the physically interesting structures [START_REF] Liseikin | Grid Generation Methods[END_REF][START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF][START_REF] Zikanov | Essential Computational Fluid Dynamics[END_REF]. The spatial discretization can be very complicated in irregular or nonuniform geometries, such as the wing flow or engine.

In this thesis, the canonical case of turbulence and reaction have been studied. The computational domain is a three dimensional cube. The spatial discretization is implemented by building uniformly structured discrete grids along the straight lines of Cartesian coordinate systems. The spatial coordinate vector x in the discrete domain is

x = i L x N x e x + j L y N y e y + k L z N z e z , (3.1)
where L x , L y , L z and N x , N y , N z are the length of the domain and the number of discrete grids in each direction;

i, j and k are integers ranging from 1 to N x , N y and N z respectively. Because of the cubic computational domain,

we have N x = N y = N z = N and L x = L y = L z = L in this thesis.
Larger N provides better spatial resolution, but meanwhile increases the cost of computation.

Pseudo-spectral methods

As proposed in K41 theory (Kolmogorov, 1941a,b) and validated by numerous works [START_REF] Sreenivasan | On local isotropy of passive scalars in turbulent shear flows[END_REF][START_REF] Sreenivasan | The phenomenology of small-scale turbulence[END_REF], for the turbulent flow at the scale much smaller than the geometric size of the domain, the information of large scale or boundary conditions are lost. Thus in the studies focusing on the universal properties of turbulent motion, the boundary condition can be as simple as possible, generally periodic. In this circumstance, pseudo-spectral methods [START_REF] Bibliography Orszag | Numerical simulation of incompressible flows within simple boundaries. i. galerkin (spectral) representations[END_REF]Eswaran and Pope, 1988b;[START_REF] Mansour | Decay of isotropic turbulence at low reynolds number[END_REF][START_REF] Sripakagorn | Extinction and reignition in a diffusion flame: a direct numerical simulation study[END_REF][START_REF] Hou | Computing nearly singular solutions using pseudo-spectral methods[END_REF] are the preferred numerical approach, because of their high accuracy in approximating the spatial derivatives.

The general idea of pseudo-spectral methods

The idea of pseudo-spectral methods is similar with the method of separation of variables used for solving partial differential equations analytically. When the method of separation of variables is used for solving PDEs, the solutions are expressed as a series of infinite number of continuous eigenfunctions with coefficients. These eigenfunctions are linearly independent and orthogonal with each other, for example cos, sin, Bessel functions and Chebyshev polynomials [START_REF] Zikanov | Essential Computational Fluid Dynamics[END_REF][START_REF] Moin | Fundamentals of Engineering Numerical Analysis[END_REF]. Then in the numerical simulation with pseudo-spectral methods, the solutions are expressed as a series of finite number of discrete eigenfunctions with coefficients, where the number of the eigenfunctions are generally equal or proportional to the number of discrete grids. For both methods, the final solutions are obtained by determining the coefficients of the eigenfunctions.

In the numerical simulation of turbulent motions in a periodic domain with a pseudo-spectral method, the eigenfunctions mostly used are the Fourier series (other functions such as Chebychev polynomials for bounded domain). Specifically for solving the Navier-Stokes equations (Eq. (1.9)) and the diffusion-convection equations of the scalars (Eq. (1.28) and (2.19)), the velocity and scalar in real space u(x, t) and θ(x, t) are represented in spectral space (also called Fourier space) by finite Fourier series:

u(x, t) = k û(k, t)e ik•x , θ(x, t) = k θ(k, t)e ik•x . (3.2)
The discrete wave number vector k is

k = 2π L (ie kx + je ky + ke kz ), (3.3)
where i, j and k are integers ranging from -N 2 + 1 to N 2 (N is generally even integer). N π L represents the computational maximum resolved wave number in Fourier space. Actually, when non-linear term is concerned, for dealiasing, the maximum resolved wave number should be smaller [START_REF] Hou | Computing nearly singular solutions using pseudo-spectral methods[END_REF]). e ik•x is the discrete complex eigenfunction at the mode k. û(k, t) and θ(k, t) are complex numbers, called the Fourier coefficients of velocity and scalar at the mode k. The objective of pseudo-spectral methods is to determine the û(k, t) and Numerical implementation of direct numerical simulations θ(k, t), which can then be transformed into real space quantities (u(x, t) and θ(x, t)) via Eq. (3.2). Eq. (3.2) is called as the discrete Fourier inverse transform for obtaining real space quantities with Fourier space quantities.

The reverse of Eq. (3.2) is called the discrete Fourier transform:

û(k, t) = x u(x, t)e -ik•x , θ(k, t) = x θ(x, t)e -ik•x , (3.4)
where x goes over all the discrete grid points. It is very convenient to approximate the spatial derivatives by using Fourier series in a pseudo-spectral method, because the operation of spatial derivation in real space is a arithmetic operation of multiplication in Fourier space. For example from Eq. ( 3.2), it can be easily found that

∇ • u(x, t) = k ik • û(k, t)e ik•x , ∂θ ∂x (x, t) = k ik x θ(k, t)e ik•x . (3.5)
It is noteworthy that the pseudo-spectral method is not called spectral method because not every term in the Navier-Stokes equations (Eq. (1.9)) and the diffusion-convection equations of the scalars (Eq. ( 1.28) and

(2.19)) are computed in Fourier space [START_REF] Bibliography Orszag | Numerical simulation of incompressible flows within simple boundaries. i. galerkin (spectral) representations[END_REF][START_REF] Orszag | Comparison of pseudospectral and spectral approximation[END_REF]. In Fourier space, it is efficient and accurate in computing the real space derivatives but not the non-linear terms (convection term and probably also reaction term), which become convolutions. In a pseudo-spectral method, it is in real space that the non-linear terms are computed, whose Fourier coefficients are then obtained by Fourier transform. By computing the non-linear term in real space, the number of operations at each time step can be proportional to N log(N ), where N represents the number of discrete grids in real space and also the maximum resolved wave number in Fourier space, instead of N 2 by computing the convolution [START_REF] Canuto | Spectral Methods in Fluid Dynamics[END_REF].

Smooth dealiasing

In the discrete Fourier transform of non-linear terms, there are high wave number modes beyond the resolved range generated. These extra high wave number modes are called aliasing error. To maintain the numerical stability, dealiasing is needed to get rid of the aliasing error. Otherwise, as shown in figure 3.1, the aliasing error can lead to the abnormal peak of the energy spectrum at the resolved high wave number modes. Moreover, the numerical results at the the low wave number modes can be also distorted. The implementation of dealiasing leads to neglect some of the largest resolvable wave number components, generally performed for each numerical time step. The most commonly method is the 2/3 spherical truncation [START_REF] Bibliography Orszag | Numerical simulation of incompressible flows within simple boundaries. i. galerkin (spectral) representations[END_REF], which means to enforce the Fourier coefficients of the modes with |k| > 2

3

N π L as 0. In this thesis, the method of smooth dealiasing was used. Instead of the sudden cut-off, a filter of the high wave number modes with a relative smooth filtering

kL N E(k) ZLWKGHDOLDVLQJ ZLWKRXWGHDOLDVLQJ
)LOWHULQJ)XQFWLRQ Specifically:

F (k) = e -α( kL N π ) β , (3.6)
where α = β = 36 is generally suggested (the dotted curve in figure 3.1). Compared with the conventional 2/3 rule approach, the smooth dealiasing is capable to reduce numerical high frequency instabilities [START_REF] Hou | Computing nearly singular solutions using pseudo-spectral methods[END_REF].

The implementation of pseudo-spectral method in this thesis

In this thesis, the pseudo-spectral method is used for solving the Navier-Stokes equations (Eq. (1.9)), also the diffusion-convection equation of the scalars (Eq. (1.28) and (2.19)) when the scalar fields are of periodic Numerical implementation of direct numerical simulations boundary conditions. Specifically, Eq. (1.9) and (2.19) are first rewritten as:

∂u ∂t = ν u + L(u), ∂θ ∂t = D θ + L(θ), (3.7)
with

L(u) = -(u • ∇)u - ∇p ρ + f , L(θ) = -(u • ∇)θ + ω + q.
In Fourier space, Eq. (3.7) are expressed as

∂ û ∂t = -ν|k| 2 û + L(û), ∂ θ ∂t = -ν|k| 2 θ + L( θ), (3.8)
where û, L(û), θ and L( θ) as the Fourier coefficients obtained from the discrete Fourier transform of u, L(u), θ and L(θ) respectively. The pressure term -∇p/ρ in L(u) is dealt with projection method (page 211 in [START_REF] Pope | Turbulent Flows[END_REF]). Then new variables ũ and θ are defined as:

ũ = û • exp(ν|k| 2 t), θ = θ • exp(ν|k| 2 t), (3.9)
which lead to

∂ ũ ∂t = ( ∂ û ∂t + ν|k| 2 û)exp(ν|k| 2 t), ∂ θ ∂t = ( ∂ θ ∂t + ν|k| 2 θ)exp(ν|k| 2 t). (3.10)
Thus the compact expressions about ũ and θ can be derived by substituting Eq. (3.9) into Eq. (3.8), as .11) Eq. (3.11) are transient PDEs, which can be solved by performing discrete time marching strating from the initial conditions. More details about the time marching are presented in section 3.2.

∂ ũ ∂t = L(ũ • exp(-ν|k| 2 t))exp(ν|k| 2 t), ∂ θ ∂t = L( θ • exp(-ν|k| 2 t))exp(ν|k| 2 t). ( 3

Finite difference methods

In addition to the pseudo-spectral method, the finite difference method (Zhao et al., 2018a,b;[START_REF] Demosthenous | Direct numerical simulations of premixed methane flame initiation by pilot n-heptane spray autoignition[END_REF] and finite volume method [START_REF] Versteeg | An Introduction to Computational Fluid Dynamics: the Finite Volume Method[END_REF] are also widely used in approximating the spatial derivatives. In this thesis, the cases with the scalar fields of Dirichlet boundary conditions, by which the scalar quantity at the boundary are prescribed, are also investigated. In such circumstance, the finite difference method was used for solving the scalar equations.

Some basic concepts about finite difference methods

As commonly used methods in the direct numerical simulation of turbulent motions [START_REF] Kristoffersen | Direct simulations of low-reynolds-number turbulent flow in a rotating channel[END_REF][START_REF] Desjardins | High order conservative finite difference scheme for variable density low mach number turbulent flows[END_REF][START_REF] Demosthenous | Direct numerical simulations of premixed methane flame initiation by pilot n-heptane spray autoignition[END_REF], the general idea of the approximation of spatial derivative using finite difference method is originated from the definition of derivatives. In a two-dimensional continuous domain with x and y as coordinate variables, the derivative of a function F (x, y) with respect to x is defined as

∂F ∂x (x, y) = lim ∆x→0 F (x + ∆x, y) -F (x, y) ∆x . (3.
12)

It is very straightforward to show that in the discrete domain, the ∂F ∂x at the grid point of (i, j) can be approximated as

∂F ∂x i,j ≈ F i+1,j -F i,j ∆x , (3.13)
where ∆x = x i+1 -x i is the grid step. When the spatial resolution tends to be infinitely fine, i.e. ∆x approaches 0, such an approximation approaches the exact derivative. The approximation of (3.13) is called first order forward scheme. It can also be derived by considering the Taylor series of F (x, y) at the grid of (x i+1 , y i ) with respect to (x i , y i ), which is

F i+1,j = F i,j + ∂F ∂x i,j ∆x + ∂ 2 F ∂x 2 i,j (∆x) 2 2 + ∂ 3 F ∂x 3 i,j (∆x) 3 6 + . . . . (3.14)
Then a precise expression of ∂F ∂x i,j is obtained as

∂F ∂x i,j = F i+1,j -F i,j ∆x + ∂ 2 F ∂x 2 i,j ∆x 2 + ∂ 3 F ∂x 3 i,j (∆x) 2 6 + . . . . (3.15)
It is clearly seen that the approximation of (3.13) can be obtained by neglecting the polynomial terms with the order of O(∆x) and higher in (3.15). These neglected terms are called truncation error. The scheme of (3.13) is of first order accuracy because its truncation error is dominated by the first order polynomial.

Construction of finite difference schemes of high order of accuracy

One advantage of finite difference methods is the convenience in constructing difference schemes of high order of accuracy, which is necessary for capturing the intermittent structure, for example in turbulence or rapid reactions.

One most commonly used method for constructing a difference scheme of high order of accuracy is by considering the Taylor expansions at more grid points, i.e. taking the information from more grid points into consideration. For example, to approximate ∂F ∂x i,j , the Taylor expansion at the grid point of (x i-1 , y j ):

F i-1,j = F i,j - ∂F ∂x i,j ∆x + ∂ 2 F ∂x 2 i,j (∆x) 2 2 - ∂ 3 F ∂x 3 i,j
(∆x) 3 6 + . . . , (3.16) can also be considered.

By subtracting Eq. (3.14) by Eq. (3.16) and neglecting the polynomial terms with order higher than O(∆x 2 ),

it yields ∂F ∂x i,j ≈ F i+1,j -F i-1,j 2∆x .
(3.17)

The approximation by Eq. ( 3.17) is of second order accuracy, called second order center scheme.

Similarly, difference schemes with higher order of accuracy or for higher order derivatives can be constructed by considering more grid points. For example, the fourth order center scheme for ∂F ∂x i,j reads

∂F ∂x i,j ≈ -F i+2,j + 8F i+1,j -8F i-1,j + F i-2,j 12∆x , (3.18)
the second order center scheme for ∂ 2 F ∂x 2 i,j reads .19) and the second order center scheme for mixed derivative ∂ 2 F ∂x∂y i,j reads

∂ 2 F ∂x 2 i,j ≈ F i+1,j -2F i,j + F i-1,j ∆x 2 , ( 3
∂ 2 F ∂x∂y i,j = 1 2∆x F i+1,j+1 -F i+1,j-1 2∆y - F i-1,j+1 -F i-1,j-1 2∆y + O((∆x) 2 , (∆y) 2 ). (3.20)

Boundary conditions

In finite difference method, for the grid points in the inner domain, the difference scheme discussed above can be directly applied to the governing equations to obtain the discrete algebraic equations. However, for the grids at the boundaries, special treatments are generally needed, because the boundary conditions are generally not periodic. Otherwise the boundary grids are essentially the same as those in the inner domain, like in pseudo-spectral methods.

There is one important requirement for the scheme used for the boundary grids, that is to have the order of accuracy not lower than that of the inner domain grids, because the accuracy of the entire numerical simulation is determined by the calculation on all the grid points. A big truncation in one grid will be propagated to other grids, and gradually to the entire domain. Thus a lower order of accuracy at the boundary will compromise

-3F 0 + 4F 1 -F 2 2∆x = β and 3F N -2 -4F N -1 + F N 2∆x = β. (3.21)
An important technique for implementing the boundary condition is using the ghost grids, i.e. defining the grids beyond the boundary, such as the grids of x -2 = -2∆x, x -1 = -∆x and x N +1 = L+∆x, x N +2 = L+2∆x.

For example with the Dirichlet boundary of F (0) = α, if the approximation of ∂F ∂x in the inner domain is implemented with the fourth order center scheme (Eq. (3.18)), we can have not only F 0 = α but also F -1 = α at the ghost grid. The benefit of doing this is that there is no special treatment needed in the approximation of ∂F ∂x 1 , which requires the value of F -1 according to Eq. (3.18). With the Neumann boundary condition ∂F ∂x 0 = β, F -1 can be prescribed as F 0 -∆xβ, which suggests a slope of the Neumann boundary condition value between the boundary grid and ghost grid.

The implementation of finite difference method in this thesis

In some configurations explored in this thesis, the scalar fields are of Dirichlet boundary conditions (Chapter 6), which requires a finite difference method for solving the scalar equations. In order to ensure the accurate approximation of the intermittent structures in the turbulent motion or fast reaction, very fine finite difference schemes are used. Specially, the first-order derivative is approximated by using eighth order upwind difference scheme, which is dependent on the local velocity. For example, at the grid point of (x i , y j , z k ), when u x;i,j,k > 0, ∂θ ∂x i,j,k is approximated as

∂θ ∂x i,j,k ≈ 1 280 ( θ i,j,k -θ i-5,j,k ∆x ) - 1 28 ( θ i,j,k -θ i-4,j,k ∆x ) + 1 6 ( θ i,j,k -θ i-3,j,k ∆x ) - 1 2 ( θ i,j,k -θ i-2,j,k ∆x ) + 5 4 ( θ i,j,k -θ i-1,j,k ∆x ) - 1 2 ( θ i,j,k -θ i+1,j,k ∆x ) + 1 14 ( θ i,j,k -θ i+2,j,k ∆x ) - 1 168 ( θ i,j,k -θ i+3,j,k ∆x ),
and when u x;i,j,k < 0, ∂θ ∂x i,j,k is approximated as

∂θ ∂x i,j,k ≈ 1 280 ( θ i+5,j,k -θ i,j,k ∆x ) - 1 28 ( θ i+4,j,k -θ i,j,k ∆x ) + 1 6 ( θ i+3,j,k -θ i,j,k ∆x ) - 1 2 ( θ i+2,j,k -θ i,j,k ∆x ) + 5 4 ( θ i+1,j,k -θ i,j,k ∆x ) - 1 2 ( θ i-1,j,k -θ i,j,k ∆x ) + 1 14 ( θ i-2,j,k -θ i,j,k ∆x ) - 1 168 ( θ i-3,j,k -θ i,j,k ∆x ).
As to the second-order derivative ∂ 2 θ ∂x 2 i,j,k , tenth order center difference scheme is used:

∂ 2 θ ∂x 2 i,j,k ≈ 5 3 ( θ i+1,j,k -2θ i,j,k + θ i-1,j,k ∆x 2 ) - 20 21 ( θ i+2,j,k -2θ i,j,k + θ i-2,j,k 4∆x 2 ) + 5 14 ( θ i+3,j,k -2θ i,j,k + θ i-3,j,k 9∆x 2 ) - 5 63 ( θ i+4,j,k -2θ i,j,k + θ i-4,j,k 16∆x 2 ) + 1 126 ( θ i+5,j,k -2θ i,j,k + θ i-5,j,k 25∆x 2
).

The polynomial coefficients in the above schemes are determined by considering the Taylor expansions of the neighbor grids, by which 8th and 10th order of accuracy are ensured for the approximations of the first and second order derivative respectively.

The Dirichlet boundary conditions are employed by prescribing the scalar quantities at the boundary grids and using 4 extra ghost grids beyond the boundary. Thereby, the above schemes can be applied to both the grids on the boundaries and the inner domain.

Temporal discretization and time integration

For a time-dependent PDE, the temporal discretization is needed for approximating the time-derivative terms.

The discretization of the temporal domain is generally implemented by discretizing the temporal coordinate t into a series of uniform time steps: t n = t 0 +n∆t with n as an integer rangeing from 0 to infinity. t 0 corresponds to the initial time, at which the quantities of the entire spatial domain are prescribed by the initial conditions. Then a specific marching scheme is applied to the governing equations for doing the numerical temporal integration, in which the known information of one time step (t n ) can be used to compute the unknown quantities at next time step (t n+1 ).

3.2 -Temporal discretization and time integration 49

Numerical marching schemes

Implicit and explicit schemes

In turbulent studies, the temporal derivative is generally linear and only of first order (Eq. ( 1.9) and (2.19)).

Therefore, the PDE about of F (representing velocity components or scalar), which is a function of spatial and temporal coordinates (x and t), can be written in the form of

∂F ∂t = L(x, F, ∂F ∂x , ∂ 2 F ∂x 2 , • • • ), (3.22)
where L includes every term of the PDE except the temporal derivative, like Eq. (3.11).

A simple way to approximate ∂F ∂t can be a scheme similar with Eq. (3.13). Moreover, after the spatial discretization of the spatial derivatives, the right hand side of Eq. ( 3.22) are approximated with algebraic expressions of the discrete F . Therefore, in a discrete form, Eq. ( 3.22) can be approximated as

F n+1 i -F n i ∆t ≈ L(• • • , F n i-1 , F n i , F n i+1 , • • • ). (3.23)
In Eq. (3.23), the quantities with superscript of n are of the time step of t n , which are known; and the superscript of n + 1 indicates the time step of t n+1 , which are unknown and to be computed. Thus Eq. ( 3.23) provides a straightforward scheme for numerical temporal integration, as

F n+1 i ≈ L(• • • , F n i-1 , F n i , F n i+1 , • • • )∆t + F n i . (3.24)
Eq. ( 3.24) is one of the explicit schemes, which means the every unknown quantity at the time step of t n+1 can be explicitly expressed as function of quantities at the time step of t n . In contrast, if the quantities used in the right hand side of Eq. (3.23) are all of the unknown time step, an implicit scheme can be derived, as

F n+1 i -L(• • • , F n+1 i-1 , F n+1 i , F n+1 i+1 , • • • )∆t ≈ F n i . (3.25)
With such a scheme, the numerical marching from t n to t n+1 is generally implemented by considering Eq. (3.25) for all the spatial grid points simultaneously and doing the computation in the form of vector and matrix.

The explicit scheme is easier for implementation than the implicit scheme, however the implicit scheme is generally more stable [START_REF] Magoulès | Computational Fluid Dynamics[END_REF], which allows a larger ∆t. Sometimes, a compromise can be made between the explicit and implicit schemes, i.e. considering the discrete quantities at both the current and next 50 Numerical implementation of direct numerical simulations time steps in the approximation of spatial derivatives. For example Eq. ( 3.22) can be approximated as

F n+1 i -F n i ∆t ≈ L(• • • , F n i-2 , F n i-1 , F n i , F n+1 i+1 , F n+1 i+2 , • • • ), (3.26)
which provides a semi-implicit temporal marching scheme.

The third order Runge-Kutta scheme

The numerical temporal integration in this thesis is implemented by using third order Runge-Kutta scheme [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes[END_REF]. For a canonical temporal differential equation

∂F ∂t = L(F ), (3.27)
the third order Runge-Kutta scheme reads

F 1 =F 0 + ∆t • L(F 0 ), F 2 = 3 4 F 0 + 1 4 F 1 + 1 4 ∆t • L(F 1 ), (3.28)
F 3 = 1 3 F 0 + 2 3 F 2 + 2 3 ∆t • L(F 2 ).
where F 3 are the unknown quantities at the time of (t + ∆t), and F 0 are the known quantities at the time t. Let us take the pseudo-spectral method as example (it is more straightforward in finite difference method).

Applying (3.28) to Eq. ( 3.11), it yields the expressions of obtaining û3 and θ3 at the time of t + ∆t with û0 and θ0 at the time of t as

û1 = (û 0 + ∆t L(û 0 )) • exp(-ν|k| 2 ∆t), û2 = 3 4 û0 • exp(-ν|k| 2 ∆t 2 ) + 1 4 û1 • exp(ν|k| 2 ∆t 2 ) + 1 4 ∆t L(û 1 ) • exp(ν|k| 2 ∆t 2 ), (3.29) û3 = 1 3 û0 • exp(-ν|k| 2 ∆t) + 2 3 û2 • exp(-ν|k| 2 ∆t 2 ) + 2 3 ∆t L(û 2 ) • exp(-ν|k| 2 ∆t 2 ), 3

Choice of the time step

Generally, a larger numerical time step ∆t is favorable for the efficiency of numerical simulation, because fewer time steps of computation are required for the simulation to converge or reach a statistically stationary state.

However, similarly with ∆x, which needs to be smaller enough for a good resolution, the ∆t also cannot be too large for the sake of numerical stability.

The stability of a numerical simulation mainly concerns the round-off error. After approximating the spatial and temporal derivatives with certain difference schemes, the continuous PDEs are transformed into some algebraic equations. These algebraic equations are solved by computers. Since the data stored in the computers are of finite digits, the final solutions we can get must be different from the exact solutions. Such difference are called the round-off error. Although the round-off error can be very small with advanced computers, if inappropriate numerical marching methods are chosen, the round-off error will be amplified at each numerical time step and finally lead to the divergence.

Therefore, the maximum ∆t can be determined by considering the temporal evolution of the round-off errors.

The two most widely used methods for the stability analysis are the Neumann method and the matrix method.

The general idea of both methods are deriving the evolution equations of the round-off error ( n ) based on the discrete governing equations and then expressing the increasing rate of the amplitude of the round-off error

| n+1 | | n | as function of ∆x and ∆t. Then under the constraint of | n+1 | | n |
1, the maximum ∆t as function of ∆x can be obtained. More details can be referred to [START_REF] Zikanov | Essential Computational Fluid Dynamics[END_REF]. It is noteworthy that the stability analysis provides the necessary but not sufficient condition for the numerical stability.

The Courant-Friedrichs-Lewy stability condition

When the time marching is implemented with explicit scheme, the value of time step is always crucial for the numerical stability. space, its physical domain of dependence must be included in the numerical domain of dependence. A good example to explain the CFL condition is the one-dimensional linear convection equation:

∂s ∂t + Λ ∂s ∂x = 0, (3.31)
where s = s(x, t) is any quantity as function of x and t; Λ is a positive constant, as the speed of propagation.

The PDE of Eq. ( 3.31) has a general solution as:

s = F (x -Λt), (3.32)
where F is any function of x-Λt. The solution of (3.32) indicates that s(x, t) is in the form of wave propagation at the constant speed of Λ. Thus for the point of (x i , t n ), the physical domain of dependence is the straight line of the slope of 1/Λ and passing it, as shown in figure 3.2. Numerically, if Eq. (3.31) is discretized with the approximation of ∂s ∂t

n i ≈ s n+1 i -s n i ∆t , ∂s ∂x n i ≈ s n i -s n i-1 ∆x , (3.33)
the time marching is implemented as

s n+1 i = -Λ ∆t ∆x (s n i -s n i-1 ) + s n i .
(3.34)

-Numerical forcing 53

The numerical result of s n+1 i is dependent on s n i and s n i-1 . Thus the numerical domain of dependence is the grey area shown in figure 3.2. In figure 3.2, the numerical domain of dependence includes the physical domain of dependence. Thus the scheme of (3.34) can be stable, although not surely. Otherwise, if the slope of ∆t/∆x is larger than 1/Λ, the numerical scheme is unstable because it fails to capture the physical information in need.

The Courant number (also called CFL number) (Le [START_REF] Veque | Numerical Methods for Conservation Laws[END_REF][START_REF] Blazek | Computational Fluid Dynamics: Principles and Applications[END_REF] is defined as

Λ ∆t ∆x . (3.35)
Quantitatively, CFL stability condition states that the CFL number should not be larger than unity.

The PDE of Eq. ( 3.31) is meaningful to the turbulence simulation because its spatial derivative term is in a similar form with the convection term of the N-S equation and the scalar equation (Eq. ( 1.9) and (1.28)), which is generally the dominant term. When the N-S equation or the scalar equation is solved with an explicit method, the CFL number can be defined as

|u| max ∆t ∆x , (3.36)
where |u| max is the global maximum of the absolute value of single-component velocity.

Numerical forcing

The energy cascade of the kinematic energy and scalar fluctuation in the inertial length scale are of high interest in turbulence studies. The source of the fluctuation comes generally from the mean part of the quantity of interest in physical space or external input in spectral space. Thus in the numerical simulations about a turbulence in statistically stationary state, lots of techniques for numerical forcing were developed (Eswaran and Pope, 1988a,b;[START_REF] Moin | Direct numerical simulation: A tool in turbulence research[END_REF][START_REF] Overholt | A deterministic forcing scheme for direct numerical simulations of turbulence[END_REF][START_REF] Alvelius | Random forcing of three-dimensional homogeneous turbulence[END_REF][START_REF] Rao | A mathematical framework for forcing turbulence applied to horizontally homogeneous stratified flow[END_REF]. Each forcing scheme has pros and cons.

The energy cascade concerns only the fluctuating part of velocity or scalar. Thus it is a natural choice to introduce the forcing from the mean part, for example by imposing a mean gradient [START_REF] Sreenivasan | On local isotropy of passive scalars in turbulent shear flows[END_REF]Gotoh and Watanabe, 2015;[START_REF] Gauding | High-order structure functions for passive scalar fed by a mean gradient[END_REF][START_REF] Iyer | Steep cliffs and saturated exponents in three-dimensional scalar turbulence[END_REF]. With the existence of a mean gradient, due to the vortex stretching as discussed in the sections of 1.1.1, the energy is extracted from the mean motions and transferred into the fluctuations of velocity by shear stress, i.e. turbulent diffusion. For the scalar, in the case of no external forcing and when the scalar is of constant mean gradient such as θ = θ + Gz, the governing Numerical implementation of direct numerical simulations equation of the variance of scalar (Eq. (1.35)) becomes

1 2 D θ 2 Dt = -D (∇θ ) 2 -G u z θ . (3.37)
The variance of scalar fluctuations can be sustained by a source proportional to the turbulent flux of scalar (covariance of scalar and velocity) [START_REF] Tennekes | A First Course in Turbulence[END_REF]. However, the forcing by a mean gradient makes one particular direction different from the others, leading to the breakdown of the local isotropy [START_REF] Sreenivasan | On local isotropy of passive scalars in turbulent shear flows[END_REF][START_REF] Biferale | Anisotropy in turbulent flows and in turbulent transport[END_REF]Gotoh and Watanabe, 2015;[START_REF] Iyer | Steep cliffs and saturated exponents in three-dimensional scalar turbulence[END_REF].

Another widely used numerical forcing is the large scales forcing, which is generally isotropic. The idea of this type of forcing is to limit the artificial information at several largest scales, without significant interference to the universal properties in the inertial range of turbulent motions. Typically these numerical forcings are implemented in Fourier space, in which the length scales can be well-identified after the discrete Fourier transform.

To have a steady supply of turbulent energy, many works simply prescribed the energy spectrum at large scale modes, for example by freezing the amplitude of the velocity in a given range of wave numbers [START_REF] Siggia | Intermittency effects in a numerical simulation of stationary threedimensional turbulence[END_REF]. The large scale motions were totally artificial. The desired spectrum in the prescribed range should be known in advance, otherwise there might be a big distortion to the problem being studied especially when the maximum wave number was not large enough. A modification can be allowing the free evolutions of individual modes in a range of scales while maintaining the average energy in this range as constant [START_REF] Chasnov | Simulation of the kolmogorov inertial subrange using an improved subgrid model[END_REF][START_REF] Machiels | Predictability of small-scale motion in isotropic fluid turbulence[END_REF]. It allowed a more physically meaningful energy spectrum. However there was a discontinuity between forced and unforced modes, which might develop to be a large numerical error [START_REF] Kerr | Higher-order derivative correlations and the alignment of small scale structures in isotropic numerical turbulence[END_REF][START_REF] Sullivan | Deterministic forcing of homogeneous, isotropic turbulence[END_REF]. There are also some other works adopting the idea of prescribing the large scale quantities. For example [START_REF] Siggia | Numerical study of small-scale intermittency in three-dimensional turbulence[END_REF] tried to model the stain rate of large scales.

In order to avoid introducing too much artificial information, instead of prescribing the large scale quantities, some works adopted a deterministic forcing, which is a linear amplification of the existing large scale fields [START_REF] Kerr | Theoretical investigation of a passive scalar such as temperature in isotropic turbulence[END_REF][START_REF] Vincent | The spatial structure and statistical properties of homogeneous turbulence[END_REF]. The factor of amplification was generally dynamically determined for controlling the numerical error and accelerating the convergence. In the work by [START_REF] Overholt | A deterministic forcing scheme for direct numerical simulations of turbulence[END_REF], the factor of amplification was determined from a damper-like system of numerical equations, in which the large scale spectrum obtained was approaching a target spectrum. Such a scheme of forcing can be considered as a modification of the method prescribing the large scale energy spectrum, achieving a relatively better convergence efficiency. There were also works proposed an amplification factor being adjusted each numerical step to ensure a constant input power [START_REF] Schumacher | Asymptotic exponents from low-reynolds-number flows[END_REF][START_REF] Ghosal | A dynamic localization model for large-eddy simulation of turbulent flows[END_REF]. For the existing velocity field in Fourier 3.4 -Direct numerical simulations of reactive turbulent flows 55 space û(k, t), the forcing f (k, t) is

f (k, t) = in û(k, t) k∈K f |û(k, t)| 2 δ k,k f , (3.38)
where K f is a set of small amplitude wave number vectors, i.e. the forced modes, for example K f = {k : |k| 2 8}; k f is any wave number vector in K f ; in is the prescribed input power. The forcing of constant power is more physically meaningful. The simulation can generally evolve stably because the dissipation rate must be equal to the prescribed input power after the simulation is converged.

However, the forcing adding in a deterministic way means that the information of past steps always remains in the forcing. It is then hard to fully lost the effect of the initial condition. Thus there were studies that considered the stochastic forcing, which is weakly dependent on the existing fields because the random process is included. The commonly used stochastic forcings are of constant amplitude [START_REF] Alvelius | Random forcing of three-dimensional homogeneous turbulence[END_REF]Gotoh and Watanabe, 2015) and random phase. For example

f (k, t) = k∈K f A k e iφ(k,t) , (3.39)
where A k is the prescribed amplitude at the mode of k and φ(k, t) is the random phase, independent for each wave number and time step. In the work of [START_REF] Watanabe | Statistics of a passive scalar in homogeneous turbulence[END_REF] and Gotoh and Watanabe (2015), the forcings for both velocity and passive scalar were added in Fourier space. The real and complex parts of the forcings were first obtained as wave number independent and time-delta Gaussian random numbers, and then normalized to ensure the constant amplitude. By adopting forcing of random phase, the convergence of simulation is accelerated since the effect of initial condition is limited. The drawback is the undeserved discontinuity in temporal evolution. Thus a compromise can be made by adopting the forcing as linear amplification of the existing fields while adjusting the phase with a random process to a certain degree [START_REF] Perlekar | Droplet size distribution in homogeneous isotropic turbulence[END_REF]. The random number for adjusting the phase can be Gaussian distributed, equally distributed, or generated from the Uhlenbeck-Ornstein process [START_REF] Wax | Selected Papers on Noise and Stochastic Processes[END_REF]Eswaran and Pope, 1988b).

Direct numerical simulations of reactive turbulent flows

DNS has been widely used for the numerical simulation of chemical reactions in turbulent flow, but most of the studies focused on the compressible case, such as turbulent combustion [START_REF] Vervisch | Direct numerical simulation of non-premixed turbulent flames[END_REF]. Relatively fewer works focused on the DNS of the reactive scalar in incompressible turbulence [START_REF] De Bortoli | Modeling and Simulation of Reactive Flows[END_REF].

kL N E(k) Q Q )LOWHULQJ)XQFWLRQ Figure 3.3:
The DNS results of the energy spectra of the reaction rates (γ 1 R 1 R n 2 -γ 2 P ) of the first order (solid line) and second order (dashed line) reactions in the form of (3.40). k = |k| is the amplitude of wave number vector. The dealiasing is implemented as smooth dealiasing, with a filtering function (Eq. (3.6)) shown as the dotted curve.

Since no change in density or heat release is considered in incompressible turbulence, the most noticeable difference between the reactive scalar and the passive scalar is the reaction term in the governing equation.

The reaction terms are in the form of product of reactants and thus nonlinear for high order reaction [START_REF] Hill | Homogeneous turbulent mixing with chemical reaction[END_REF][START_REF] Heeb | Turbulent mixing with multiple second-order chemical reactions[END_REF]. These nonlinear reaction terms can introduce extra convolutions in the high wave number modes of reactive scalars, which can be physically meaningful especially with large reaction rates.

When the DNS is implemented with a pseudo-spectral method, the maximum resolved wave number is needed to be large enough to resolve the non-linear reaction terms. For example, consider the reversible reaction of

R 1 + nR 2 γ 1 G GGGGG B F GGGGG G γ 2 P (3.40)
in homogeneous isotropic turbulence (Wu et al., 2020) with n = 0 or 1, corresponding to first and second order reaction respectively. The net reaction rate is γ 1 R 1 R n 2 -γ 2 P , whose energy spectra obtained from DNS are shown in figure 3.3. It is clearly seen that for the second order reaction, there is a peak of the spectrum at the range of high wave number modes. However, the further development of this peak to smaller scales can not be reflected in DNS due to the dealiasing, i.e. insufficient resolved wave number.

In the DNS of turbulent velocity and passive scalar fields, the spatial resolution is required to be fine enough The reaction zone is like a flame front [START_REF] Domingo | Triple flames and partially premixed combustion in autoignition of nonpremixed turbulent mixtures[END_REF][START_REF] Wang | Flame edge statistics in turbulent combustion[END_REF], through which the concentration of reactants and the reaction rate can sharply change. To capture the information in the reaction zone, finite difference method with schemes of high accuracy order [START_REF] Craske | Energy dispersion in turbulent jets. part 1. direct simulation of steady and unsteady jets[END_REF]Zhao et al., 2018a) can be used. It generally also requires significantly more than 1 grids in the length of δ. For a canonical reaction R 1 + nR 2 γ -→ ∅ in turbulence, an estimation about the thickness of the reaction zone δ can be done by considering the balance between reaction and diffusion (convection is less important in the small scale).

Imagine the reaction zone as a flat cuboid of the lateral area of A and the thickness of δ (figure 3.4). In an infinitesimally small time period of τ , the R 1 consumed in the cuboid by reaction is equal to γR 1 R n 2 τ δA. In addition, according to the Fick's law of diffusion [START_REF] Bergman | Fundamentals of Heat and Mass Transfer[END_REF], there is a flux of R 1 entering the cuboid due to the gradient: -D ∂R1 ∂x Aτ . A balance is then built as Such a relation is similar to that derived for the thickness of flame front [START_REF] Law | Combustion Physics[END_REF]. Eq. ( 3.42) suggests a finer spatial resolution for faster reaction, with the grid size inversely proportional to the square root of the reaction rate.

-D ∂R 1 ∂x Aτ -γR 1 R n 2 δAτ = 0. ( 3
Chapter 4

Reactions in homogeneous isotropic turbulence

In this chapter, we focus on the fundamental properties of reactive scalar mixing in homogeneous isotropic incompressible turbulence. The flow statistics, from global to scale-dependent features are studied in details and theoretically modelled. This chapter is organized as follows: In section 4.1, we introduce the model system, its governing equations together with the set of dimensionless control parameters. Section 4.2 briefly specifies the numerical methods adopted for the numerical simulation. Then the DNS results and their analyses are presented in sections 4. 3-4.7. Finally, a brief summary about this chapter is given in section 4.8. *

Problem definition

In this chapter we consider reactions of the form:

R 1 + nR 2 γ 1 G GGGGG B F GGGGG G γ 2 P, (4.1)
where R 1 , R 2 and P denote three generic reactive scalars and n is an integer coefficient. The process is reversible with independent non-zero forward/backward reaction rates γ 1 and γ 2 . The order of the chemical reaction, which is defined as the sum of the powers of the reactants' concentration in the rate equation is n + 1 for the forward reaction, because the rate equation reads γ 1 R 1 R n 2 , while is of the first order for the backward reaction with 60

Reactions in homogeneous isotropic turbulence reaction rate γ 2 P . The reactants are assumed to be subject to molecular diffusion and to fluid advection. The evolution equations for the velocity and the concentration fields R 1 (x, t), R 2 (x, t) and P (x, t) read:

∂u ∂t + (u • ∇)u = ν u -∇p/ρ + f , (4.2) ∇ • u = 0, (4.3)
and

∂R 1 ∂t + (u • ∇)R 1 = D R 1 -γ 1 R 1 R n 2 + γ 2 P + qR1 , (4.4a) ∂R 2 ∂t + (u • ∇)R 2 = D R 2 -n(γ 1 R 1 R n 2 -γ 2 P ) + qR2 , (4.4b) ∂P ∂t + (u • ∇)P = D P + γ 1 R 1 R n 2 -γ 2 P + qP . (4.4c)
Here u(x, t) is the three-dimensional flow velocity, p is the pressure, ρ is the fluid density set as constant, ν is the kinematic viscosity and D is the species diffusivity (assumed as being the same for all species). To sustain the turbulent fluctuations, large-scale forcing terms f and q are introduced for the velocity and scalars, respectively.

More details about the expression of these forcing terms will be provided in section 4.2.

For comparison, a non-reactive species T undergoing both advection and diffusion is also considered. Its local concentration evolves according to the following equation,

∂T ∂t + (u • ∇)T = D T + qT . (4.5)
The equations for the above model system can be made dimensionless by choosing reference scales appropriate for the present system. Since the turbulent flow is unbounded, we take the Taylor microscale (λ) and the single component velocity fluctuation (u ) as the reference scales for space and velocity, respectively, which are defined

as: λ = 15ν ε u , u = 1 3 i u 2 i , ε = ν 2 i j (∂ i u j + ∂ j u i ) 2 .
Here ε is the mean dissipation rate; • is the ensemble average, but numerically represented by space and time average in this chapter without special notation. The scalar quantities can be non-dimensionalized by means of their equilibrium values in no-flow conditions R 1,eq , R 2,eq , P eq , while for the passive scalar the global mean T is used as the reference value. Note that at the equilibrium, the algebraic relation γ 1 R 1,eq R n 2,eq = γ 2 P eq holds. Furthermore, in this chapter for simplicity we assume that R 1,eq = R 2,eq . This leads to the following dimensionless equations:

∂u ∂t + (u • ∇)u = Re -1 λ u -∇p + f , (4.6) 4.2 -Numerical methods 61 ∇ • u = 0, (4.7) ∂R 1 ∂t + (u • ∇)R 1 = (Sc Re λ ) -1 R 1 -Da(R 1 R n 2 -P ) + qR1 , (4.8a) ∂R 2 ∂t + (u • ∇)R 2 = (Sc Re λ ) -1 R 2 -nDa(R 1 R n 2 -P ) + qR2 , (4.8b) ∂P ∂t + (u • ∇)P = (Sc Re λ ) -1 P + Da(R 1 R n 2 -P ) + qP , (4.8c) ∂T ∂t + (u • ∇)T = (Sc Re λ ) -1 T + qT , (4.8d)
where Re λ = λ • u /ν is the Taylor based Reynolds number, the Schmidt number Sc = ν/D is the ratio of viscous diffusion to molecular diffusion, the Damköhler number Da = λγ 1 R n 2,eq /u = λγ 2 /u represents the ratio of flow timescale to the chemical timescale of forward or backward reaction. Note that the particular choice R 1,eq = R 2,eq is crucial in obtaining a single Damköhler number, instead of two distinct ones that would be present in general cases.

In conclusion the control parameters of the model system are Re λ , Sc, n and Da.

Numerical methods

The model system presented in Sec. 4.1 is numerically simulated in a cubic tri-periodic domain. The flow is sustained by a large-scale forcing capable to generate a statistically steady homogeneous and isotropic turbulent flow. The expression of the forcing field in Fourier space, f (k, t) reads,

f (k, t) = 1 τ f 1 |k| 2 √ 2 û(k, t), (4.9)
with τ f a time-scale being adjusted at each time step in order to provide a constant power input, i.e. V f •udx 3 = const.. This type of forcing, called linear, has been adopted e.g. in [START_REF] Schumacher | Asymptotic exponents from low-reynolds-number flows[END_REF]. Note also that the zero mode |k| = 0 is not forced in order to prevent the development of a global mean flow, i.e., in our simulations u = 0. Similarly, the external source term on scalars ( qθ with θ = R 1 , R 2 , P or T in Eq. ( 4.8)) is also isotropic and acting at the largest scales; however it is constant in amplitude [START_REF] Alvelius | Random forcing of three-dimensional homogeneous turbulence[END_REF]Gotoh et al., Gotoh and Watanabe, 2015). In Fourier space this reads,

qθ ( k, t) = 1 | k| 2 √ 2 Q | k| e iφ θ ( k,t) , (4.10)
where Q is the constant prescribing the overall source amplitude, |k| -1 is a normalization factor to guarantee that the forcing amplitude is larger at small wave numbers. In particular, the random phase function φ θ (k, t)

is generated independently for each scalar field and delta-correlated both in time and in wave-vector [START_REF] Gotoh | Universality and anisotropy in passive scalar fluctuations in turbulence with uniform mean gradient[END_REF]Gotoh and Watanabe, 2015). As a result, qR1 , qR2 , qP and qT have amplitudes of the same order, but they are statistically independent from each other both in time and in space.

The set of dynamical equations (4.6) and (4.8) are solved numerically by means of a pseudo-spectral code [START_REF] Gauding | High-order structure functions for passive scalar fed by a mean gradient[END_REF][START_REF] Gauding | One-point and two-point statistics of homogeneous isotropic decaying turbulence with variable viscosity[END_REF], using a smooth dealiasing technique [START_REF] Hou | Computing nearly singular solutions using pseudo-spectral methods[END_REF] for the treatment of non-linear terms in the equations (Section 3.1.1).

We explore the parameter space of the problem by means of a series of simulations: the Reynolds number

t/(k/

UPV R 1 2 1/2 R 2 2 1/2 P 2 1/2 T 2 1/2 R1 R2 P T PHDQ

Spectra and coherency spectra of scalars

In this section we focus on the scale-dependent behavior of reactive scalar fluctuations and their mutual correlations. The energy spectra of the velocity and scalars are defined as for the velocity and for a passive scalar field, i.e., the KOC scaling Kolmogorov (1941b,a). The spectra are also normalized by the total energy for each scalar. The figure shows that the spectra are indistinguishable from the ones of a passive scalar, and display the same scaling in the inertial range. This proves that in the present condition the reaction terms have a negligible effect on the scalar energy transfer. Remarkably, this behaviour is also Damköhler number independent. The latter observation is qualitatively confirmed also by visualisations of the instantaneous scalar fields for different Da values. Despite perceptible larger fluctuations for the small

Energy spectra

k E (k) 2 (k ) 5/3 7 R 1 R 2 P
E(k) = 4πk 2 1 2 ûi (k)û * i (k) k , (4.11a) E θ (k) = 4πk 2 θ(k) θ * (k) k , θ = R 1 ,
Da case (the one where the chemistry is slower) it appears that the spatial structure of the fields is not affected by the magnitude of Da.

Coherency spectra

The coherency spectrum between two scalar fields θ 1 and θ 2 is defined as

Co θ1,θ2 (k) = | θ 1 (k) θ 2 * (k)| k θ 1 (k) θ 1 * (k) k θ 2 (k) θ 2 * (k) k . ( 4

Global correlation coefficients of reactive scalars

In this section the global correlation coefficients for the scalar fields is investigated. We begin with a theoretical argument for the prediction of its functional dependence on varying the dimensionless a priori control parameters Re λ , Sc, Da and n. We will later compare the prediction with the numerical results.

The global correlation coefficients between reactive scalars are defined as

r(θ 1 , θ 2 ) = θ 1 θ 2 θ 2 1 1/2 θ 2 2 1/2 . (4.13)
Here θ 1 and θ 2 are the scalars under consideration.

The global value of the cross product of scalar gradients ∇θ 1 ∇θ 2 can be estimated in terms of the global value of cross product of scalars θ 1 θ 2 normalized by the square of a characteristic length-scale λ θ . For a scalar quantity θ, define λ θ as in [START_REF] Ristorcelli | Passive scalar mixing: Analytic study of time scale ratio, variance, and mix rate[END_REF]:

λ 2 θ = θ 2 (∂ x θ) 2 . (4.14)
Such a length scalar can be interpreted as the Taylor microscale of θ. Consequently, this allows to introduce an a posteriori control parameter, the Damköhler number based on scalar Taylor micro-scale and diffusivity, denoted here Da θ , which is defined as

Da θ = Re λ ScDaλ 2 θ = λ 2 θ D γ 1 R n 2,eq . (4.15)
It has to be noted that such a number includes a combination of the three a priori control parameters (Re λ ,Sc,Da) for the model system, with the addition of the λ θ scale, which therefore plays a key role in the analysis.

The theoretical prediction is based on the following two hypotheses. First, given the fact that the reactive 4.5 -Global correlation coefficients of reactive scalars 67 scalar fluctuations are small with respect to the equilibrium global value, the chemical sources can then be linearized in the following way.

R 1 R n 2 -P ≈ ( R 1 + R 1 )( R 2 + n R 2 n-1 R 2 ) -( P + P ) ≈ R 1 + nR 2 -P , (4.16)
where we have used also the fact that R 1 ≈ R 2 ≈ P ≈ 1 (see again figure 4.1) ). Second, we assume that there are no correlations between the source term for a given scalar and other reactive scalars, which is here reasonably guaranteed from the fact that the source terms are delta-correlated in time with a fixed amplitude. Such an assumption, however, is not a general feature of reactive turbulence, and needs to be considered specifically.

In the following sections, theoretical prediction is presentation step be step in section 4.5.1, 4.5.2, and 4.5.3.

And the comparison with number results is in section 4.5.4.

Relation between correlation coefficients and coherency spectra

First, we prove that in current configuration, the global correlation coefficient between two scalars is the same as their correlation coefficient at each length scale.

In a two reactive scalar system where fluctuations are introduced by external perturbations, e.g., randomly introduced in the phase space, in the simulations for this chapter. It is reasonable to assume that:

1. the energy distribution of R 1 is the same as that of R 2 on different length scales;

2. the correlation coefficient between R 1 and R 2 conditional on different length scales is invariant.

Consider the one-dimensional case and analysis in three-dimensional space can be implemented similarly.

The fluctuating parts of R 1 and R 2 are expressed in the form of Fourier modes as

R 1 = k A k (t) sin(kx) + a k (t) cos(kx), R 2 = k B k (t) sin(kx) + b k (t) cos(kx). (4.17)
The global correlation coefficient is

r(R 1 , R 2 ) = [ k A k (t) sin(kx) + a k (t) cos(kx)][ k B k (t) sin(kx) + b k (t) cos(kx)] [ k A k (t) sin(kx) + a k (t) cos(kx)] 2 1/2 [ k B k (t) sin(kx) + b k (t) cos(kx)] 2 1/2 = k A k B k + a k b k [ k (A 2 k + a 2 k )] 1/2 [ k (B 2 k + b 2 k )] 1/2 , (4.18)
where • denotes time average.

Reactions in homogeneous isotropic turbulence

Since the energy distributions of R 1 and R 2 on different length scales are the same (assumption 1), we define

α as A 2 1 + a 2 1 B 2 1 + b 2 1 = A 2 2 + a 2 2 B 2 2 + b 2 2 = • • • = k (A 2 k + a 2 k ) k (B 2 k + b 2 k ) = α 2 .
Then the denominator of Eq. ( 4.18) can be further written as

[ k (A 2 k + a 2 k )] 1/2 [ k (B 2 k + b 2 k )] 1/2 = ( k B 2 k + b 2 k ) • α = k α(B 2 k + b 2 k ) = k (A 2 k + a 2 k ) 1/2 (B 2 k + b 2 k ) 1/2 .
Thus it yields

r(R 1 , R 2 ) = k A k B k + a k b k k (A 2 k + a 2 k ) 1/2 (B 2 k + b 2 k ) 1/2 . (4.19)
The coherency spectrum between R 1 and R 2 , Co R1,R2 , describes the correlation coefficients between two scalars corresponding to each length scale. At the mode k,

Co R1,R2 (k) = [A k (t) sin(kx) + a k (t) cos(kx)][B k (t) sin(kx) + b k (t) cos(kx)] [A k (t) sin(kx) + a k (t) cos(kx)] 2 1/2 [B k (t) sin(kx) + b k (t) cos(kx)] 2 1/2 = A k B k + a k b k (A 2 k + a 2 k ) 1/2 (B 2 k + b 2 k ) 1/2 . (4.20)
From the second assumption that the correlation coefficients of R 1 and R 2 are the same at each length scales, for any k, we obtain

Co R1,R2 (k) = A 1 B 1 + a 1 b 1 (A 2 1 + a 2 1 ) 1/2 (B 2 1 + b 2 1 ) 1/2 = A 2 B 2 + a 2 b 2 (A 2 2 + a 2 2 ) 1/2 (B 2 2 + b 2 2 ) 1/2 = . . . = k A k B k + a k b k k (A 2 k + a 2 k ) 1/2 (B 2 k + b 2 k ) 1/2 = r(R 1 , R 2 ). (4.21)
Therefore, the global correlation coefficient between R 1 and R 2 is the same as the correlation coefficient at each wave number or length scale.

Correlation coefficients of reactive scalar gradients

As shown in figure 4.2 and figure 4.3, the energy distribution and the correlation coefficients of the reactive scalars remain almost independent of the length scale. Under these two conditions, it is ready to derive that the global correlation coefficient between two reactive scalars is the same as their coherency spectrum at each length scale (section 4.5.1).

-Global correlation coefficients of reactive scalars 69

Another quantity of primary importance is the global correlation coefficients of the gradients of scalars, which can be defined along one direction (e.g. x) only because of isotropy:

r(θ 1,x , θ 2,x ) = ∂θ1 ∂x ∂θ2 ∂x ( ∂θ1 ∂x ) 2 1/2 ( ∂θ2 ∂x ) 2 1/2 . (4.22)
For various scalars, the (almost) identical spectra of the scalar energy (figure 4.2) implies the (almost) identical spectra of the energy of scalar gradient quantities. In addition, since Eq. ( 4.12) is the definition of the coherency spectrum between not only θ 1 and θ 2 but also their gradients, the coherency spectra between the gradients of two reactive scalars are also almost k independent. Therefore, the correlation coefficient of the gradients of two reactive scalars is also identical at each length scale, and supposed to be the same as the correlation coefficient of these two reactive scalars. between the reactive scalars are the same as that of their gradients is well testified, except for the cases of Da = 0.01 and 0.1, in which the condition that the coherency spectra is k independent is not satisfied at the largest scales (figure 4.3). 'D

r(R 1 , R 2 ) r(R 1 , P) r(R 2 , P) r(R 1, x , R 2, x ) r(R 1, x , P x ) r(R 2, x , P x )

Analytical prediction for reactant correlations

Under the conditions of

R 1 R 1 , R 2 R 2 , P P and R 1 ≈ R 2 ≈ P ≈ 1 (figure 4.1)
, the net reaction rate can be estimated as

R 1 R n 2 -P ≈ ( R 1 + R 1 )( R 2 + n R 2 n-1 R 2 ) -( P + P ) ≈ R 1 + nR 2 -P . (4.23)
Correspondingly, the equations for fluctuating scalars become

DR 1 Dt ≈ 1 Re λ Sc R 1 -Da(R 1 + nR 2 + M ) + qR1 , (4.24a) DR 2 Dt ≈ 1 Re λ Sc R 2 -nDa(R 1 + nR 2 + M ) + qR2 , (4.24b) DM Dt ≈ 1 Re λ Sc M -Da(R 1 + nR 2 + M ) + qM , (4.24c)
where M = -P and qM = -qP . The same form of Eq. ( 4.24c) as Eq. ( 4.24c) implies that R 1 and M behave statistically the same.

By multiplying Eq. (4.24a), (4.24b) and (4.24c) with R 2 , M and R 1 respectively and averaging ( • ) on time and space, it yields:

DR 1 Dt R 2 + 1 Re λ Sc ∇R 2 ∇R 1 ≈ -Da( R 2 R 1 + n R 2 2 + R 2 M ) + R 2 qR1 , (4.25a) DR 2 Dt M + 1 Re λ Sc ∇M ∇R 2 ≈ -nDa( M R 1 + n M R 2 + M 2 ) + M qR2 , (4.25b) DM Dt R 1 + 1 Re λ Sc ∇R 1 ∇M ≈ -Da( R 2 1 + n R 1 R 2 + R 1 M ) + R 1 qM . (4.25c)
In Eq. (4.25a), the term R 2 qR1 is estimated as 0, because the time-delta forcing to one scalar (R 1 ) can not be strongly correlated with another scalar (R 2 ). Moreover, the small net reaction rate implies that the instantaneous reactive scalar is weakly influenced by other scalar(s). Thus at the statistical stationary state,

DR 1
Dt R 2 can be assumed negligibly small, i.e.

DR 1 Dt R 2 ∼ 0.
As discussed in section 4.5.2, the correlation coefficients of the reactive scalars is roughly the same as the 4.5 -Global correlation coefficients of reactive scalars 71 correlation coefficients of their gradients, which in isotropic turbulence can be estimated as

∂ x R 2 ∂ x R 1 ≈ ∂ y R 2 ∂ y R 1 ≈ ∂ z R 2 ∂ z R 1 ≈ (∂ x R 2 ) 2 1/2 (∂ x R 1 ) 2 1/2 • R 2 R 1 R 2 2 1/2 R 2 1 1/2 ≈ R 2 2 1/2 R 2 1 1/2 λ 2 θ • R 2 R 1 R 2 2 1/2 R 2 1 1/2 ≈ R 2 R 1 λ 2 θ , (4.26)
where λ θ is the Taylor microscale for scalars (Eq. ( 4.14)).

Consequently,

1 Re λ Sc ∇R 2 ∇R 1 = 1 Re λ Sc ( ∂ x R 2 ∂ x R 1 + ∂ y R 2 ∂ y R 1 + ∂ z R 2 ∂ z R 1 ) ≈ 3 Re λ Scλ 2 θ R 2 R 1 .
(4.27)

From Eq. ( 4.25a) it yields

3 R 2 R 1 ≈ -Da θ ( R 2 R 1 + n R 2 2 + R 2 M ), (4.28)
where

Da θ = Re λ Scλ 2 θ Da = λ 2 θ γ1R n 2,eq

D

is the Damköhler number based on scalar Taylor micro-scale and diffusivity (Eq. ( 4.15)).

Similarly, from Eq. (4.25b) and (4.25c),

3 M R 2 ≈ -nDa θ ( M R 1 + n M R 2 + M 2 ), (4.29) 3 R 1 M ≈ -Da θ ( R 2 1 + n R 1 R 2 + R 1 M ). (4.30)
Because R 1 and M are statistically the same (see Eq. (4.24a) and (4.24c)), we can define

C = R 1 M , c = R 2 R 1 = R 2 M , V = R 2 1 = M 2 , v = R 2 2 , (4.31)
i.e.

C V = r(R 1 , M ) = -r(R 1 , P ), c √ V v = r(R 1 , R 2 ) = r(R 2 , M ) = -r(R 2 , P ).
Then Eq. (4.28), (4.29) and (4.30) can be rewritten as

3c ≈ -Da θ (2c + nv), (4.32a) 3c ≈ -nDa θ (C + nc + V ), (4.32b)
which then leads to the solutions as

r(R 1 , P ) ≈ Da θ 3 + n 2 Da θ + Da θ , (4.33a
) We provide here the central result of the derivation based on the above steps. The above expressions of Eq.

r(R 1 , R 2 ) = -r(R 2 , P ) ≈ -nDa θ √ 3 + n 2 Da θ + Da θ √ 3 + 2Da θ . (4.33b)

Comparison with numerical results

Da

D r(R1, R2)7KHRU\ r(R1, P) r(R2, P)7KHRU\ r(R 1 , R 2 )Re r(R 1 , P)Re r(R 2 , P)Re r(R 1 , R 2 )Re r(R 1 , P)Re r(R 2 , P)Re r(R 1 , R 2 )Re r(R 1 , P)Re r(R 2 , P)Re r(R 1 , R 2 )Re r(R
(4.33a) and (4.33b) show that the concentration field for R 1 is positively correlated to P , while is negatively correlated to R 2 . Furthermore, the correlations r(R 1 , R 2 ) and r(R 2 , P ) are opposite in sign. They also show that for large Da θ the correlations reach a saturation plateau, whose value depends on the reaction order n. For n = 1 the correlations coefficients have all the same intensity and only differs in sign. The asymptotically large Da θ limit in this case leads to the values r = ±1/2. At asymptotically large n and Da θ , it yields r(R 1 , P ) 0 and r(R 1 , R 2 ) = -r(R 2 , P ) 1. On the opposite, in the condition of vanishing values of Da θ , corresponding to a negligible role of chemical processes and predominance of mixing, all the correlations coefficients tend to zero.

In figure 4.5 we report the numerical measurements of the correlation coefficients between the reactive scalars 4.6 -Reactant variances 73 as functions of Da θ , for a set of simulations characterized by different reaction order n and different Re λ , ranging over more than a decade. It can be seen that for low Da θ values, corresponding to slow reaction rates, R 1 , R 2 and P behave as almost independent passive scalars, as expected. Thus the correlation coefficients are about zero when Da θ is small. As Da θ increases, the scalars become more and more correlated and the correlation coefficients gradually approach plateaus, which are n dependent. All these tendencies are in excellent agreement with the theoretical predictions, which are also reported on the same figure.

Da D r(R1, R2)7KHRU\ r(R1, P) r(R2, P)7KHRU\ r(R 1 , R 2 )Sc = 0.1 r(R 1 , P)Sc = 0.1 r(R 2 , P)Sc = 0.1 r(R 1 , R 2 )Sc = 0.25 r(R 1 , P)Sc = 0.25 r(R 2 , P)Sc = 0.25 r(R 1 , R 2 )Sc = 1 r(R 1 , P)Sc = 1 r(R 2 , P)Sc = 1 r(R 1 , R 2 )Sc = 4 r(R 1 , P)Sc = 4 r(R 2 , P)Sc = 4 Da E r(R1, R2)7KHRU\ r(R1, P) r(R2, P)7KHRU\ r(R1, R2)7KHRU\ r(R1, P) r(R2, P)7KHRU\ Da F r(R1, R2)7KHRU\ r(R1, P) r(R2, P)7KHRU\ r(R1, R2)7KHRU\ r(R1, P) r(R2, P)7KHRU\

Reactant variances

Chemical equilibrium and the effect of reaction on it

The probability density functions (PDF) of the scalar fields are reported in figure 4.7. It can be seen that despite the different amplitudes of the standard deviations, their normalised shapes do not deviate significantly from Gaussian. No noticeable difference is observed in the comparison of reactive scalars with the passive one.

Furthermore, side-by-side visualisations of instantaneous snapshots of reactive and passive scalars do not allow to perceive clear difference in their spatial structure (figure 4.8). In our system the global chemical equilibrium condition R 1 R n 2 = P is toughly satisfied by construction due to the irreversible reactive scalar source term. We also note that in the above described statistically steady condition the advecting flow is the main responsible for local departures in time and space from the chemical equilibrium condition. The magnitude of such deviations depends on the turbulent intensity and it grows with Re λ . We observe that the root-mean-square intensity of the local reaction rate Da(R 1 R n 2 -P ) scales approximately as ∼ Re

3')

T R 1 R 2 P / 2 1/2
3/2 λ
in the range of Reynolds number explored in this chapter (figure 4.9). Moreover, as shown in the bottom panels of figure 4.9, it is interesting to see that the reaction rate normalized by its standard deviation is of a PDF independent of Re λ .

Analytical prediction for reactant variances and comparison with DNS

As we have already mentioned the turbulent advection and the scalar forcing are the sources of scalar spatialtemporal fluctuations. In the case of passive scalars such fluctuations are smoothed out by diffusion. For the reactive case, the chemical sources function as an additional dumping mechanism. In other words, it is expected that the chemical reaction term acts as a global sink to suppress the scalar energy in addition to dissipation via molecular diffusion.

In the present model system this scenario can be understood by means of the following argument. We multiply the linearized transport equations for the reactive scalars (Eq. ( 4 obtained equations of scalar dissipation rates (ε R1 , ε R2 and ε P ) for the considered reactive fields reads:

ε R1 + ε R2 + ε P = (∇R 1 ) 2 Re λ Sc + (∇R 2 ) 2 Re λ Sc + (∇P ) 2 Re λ Sc ≈ -Da R 1 + nR 2 -P 2 + R 1 qR1 + R 2 qR2 + P qP . (4.34)
The above equation shows that the reaction is always responsible of removing the scalar energy. A consequence of this is that one expects smaller scalar fluctuations for the reactive fields as compared to a passive scalar. In particular, we expect that the scalar variance will be a monotonically decreasing function in Da θ .

To have a quantitative understanding for such scenario we compare the fluctuations of the reactive scalars with the ones of a passive scalar in the same dynamical conditions, i.e., subject to the same advective flow, and having the same diffusion and under the effect of an independent statistical realization of the source term qθ .

In order to develop also in this case a quantitative prediction for the phenomenon we need to introduce the key assumption that the scalar energy input due to the source term on the field R 1 is approximately same as the one provided on a passive scalar field T in the same conditions, i.e.

R 1 qR1 T qT . (4.35)
The soundness of this hypothesis lies on the fact that in the present conditions the reactant R 1 has fluctuation of similar intensity as the passive scalar case.

Now we derive the theoretical expressions of the variance of reactant scalars normalized by that of the passive scalar (T ).

Reactions in homogeneous isotropic turbulence By multiplying Eq. (4.24a), (4.24b) and (4.24c) with R 1 , R 2 and M , respectively and averaging ( • ) in space and time, we obtain

Da(P R 1 R n 2 ) PDF Q Re Re Re Re Re [Da(P R1R n 2 )] 2 Da(P R 1 R n 2 ) PDF Q Re Re Re Re Re [Da(P R1R n 2 )] 2 Da(P R 1 R n 2 )/ [Da(P R 1 R n 2 )] 2 1/2 PDF Q Re Re Re Re Da(P R 1 R n 2 )/ [Da(P R 1 R n 2 )] 2 1/2 PDF Q Re Re Re Re
1 2 D R 2 1 Dt + 1 Re λ Sc |∇R 1 | 2 ≈ -Da( R 2 1 + n R 1 R 2 + R 1 M ) + R 1 qR1 , (4.36a) 1 2 D R 2 2 Dt + 1 Re λ Sc |∇R 2 | 2 ≈ -nDa( R 1 R 2 + n R 2 2 + R 2 M ) + R 2 qR2 , (4.36b) 1 2 D M 2 Dt + 1 Re λ Sc |∇M | 2 ≈ -Da( R 1 M + n R 2 M + M 2 ) + M qM . (4.36c)
The dissipation terms above, e.g.

3 R 2 1 ≈ -Da θ ( R 2 1 + n R 1 R 2 + R 1 M ) + Re λ Scλ 2 θ R 1 qR1 , (4.37a) 3 R 2 2 ≈ -nDa θ ( R 1 R 2 + n R 2 2 + R 2 M ) + Re λ Scλ 2 θ R 2 qR2 , (4.37b) 3 M 2 ≈ -Da θ ( R 1 M + n R 2 M + M 2 ) + Re λ Scλ 2 θ M qM . (4.37c)
Similarly, based on Eq. ( 4.8d), the turbulent energy of T is

3 T 2 ≈ Re λ Scλ 2 θ T qT . (4.38) Define w = Re λ Scλ 2 θ R 2 qR2 , W = Re λ Scλ 2 θ R 1 qR1 = Re λ Scλ 2 θ M qM , V T = T 2 , W T = Re λ Scλ 2 θ T qT . (4.39)
Together with Eq. (4.31), (4.37a), (4.37b) and ( 4.38) we obtain

3V ≈ -Da θ (V + nc + C) + W, (4.40a) 3v ≈ -nDa θ (2c + v) + w, (4.40b)
3V T ≈ W T . (4.40c)
It is worthy noting that for all the scalar quantities the delta-correlated external forcing is exerted in the same way with constant amplitude. When R 2 1 is close to T 2 ( V V T close to 1), it is reasonable to assume W ≈ W T . Together with Eq. (4.33a), (4.33b), (4.40a) and (4.40c) we obtain

R 2 1 T 2 = P 2 T 2 = V V T ≈ 3 + n 2 Da θ + Da θ 3 + n 2 Da θ + 2Da θ . (4.41)
From Eq. (4.32), the ratio between the fluctuation magnitudes of R 2 and R 1 is determined as

R 2 2 R 2 1 = v V ≈ 3 + 2Da θ 3 + n 2 Da θ + Da θ , (4.42)
Reactions in homogeneous isotropic turbulence

which leads to R 2 2 T 2 = R 2 1 T 2 R 2 2 R 2 1 = V V T v V ≈ 3 + 2Da θ 3 + n 2 Da θ + 2Da θ . (4.43)
It is found that the fluctuations of the reactive scalars (R 1 , R 2 and P ) are close to that of passive scalar (T )

when Da θ is small, but gradually decrease as Da θ increases. The figure 4.10 shows the normalized fluctuations of the reactive scalars (R 1 , R 2 and P ) as measured from the DNS, in agreement with the above theoretical prediction.

Da D R 1 2 T 2 R 2 2 T 2 P 2 T 2 7KHRU\ R 1 2 T 2 Re R 2 2 T 2 Re P 2 T 2 Re R 1 2 T 2 Re R 2 2 T 2 Re P 2 T 2 Re R 1 2 T 2 Re R 2 2 T 2 Re P 2 T 2 Re R 1 2 T 2 Re R 2 2 T 2 Re P 2 T 2 Re Da E R 2 2 T 2 7KHRU\ R 1 2 T 2 P 2 T 2 7KHRU\ R 2 2 T 2 7KHRU\ R 1 2 T 2 P 2 T 2 7KHRU\ Da F R 2 2 T 2 7KHRU\ R 1 2 T 2 P 2 T 2 7KHRU\ R 2 2 T 2 7KHRU\ R 1 2 T 2 P 2 T 2 7KHRU\

Taylor micro scale of scalar concentration fields

As we have discussed in the above sections the correlations and fluctuations of the concentration of reactive scalars are well described by means of the control parameter Da θ , which contains the scalar Taylor micro scale of λ θ . In this section we aim at gaining more insight into this key spatial scale.

The scalar Taylor microscale was notably first studied by Stanly Corrsin [START_REF] Corrsin | Simple theory of an idealized turbulent mixer[END_REF], in the context of scalar mixing in turbulent flows. He hypothesized that such a scale is proportional to the intensity of turbulence and inversely to the Schmidt number of the problem, i.e.

λ 2 θ λ 2 ∝ 1 Sc .
However, further experimental and numerical studies [START_REF] Bahri | Self-similarity of passive scalar flow in grid turbulence with a mean cross-stream gradient[END_REF][START_REF] Corrsin | The isotropic turbulent mixer: Part ii. arbitrary schmidt number[END_REF] have reported that such a dependence is not straightforward, as it shows finite Re λ effects and different trends for asymptotically small and large Sc values (see [START_REF] Ristorcelli | Passive scalar mixing: Analytic study of time scale ratio, variance, and mix rate[END_REF] for a recent discussion).

We show the results of our simulations in figure 4.11(a). The figure reports the dimensionless λ 2 θ (actually λ 2 θ /λ 2 , because the length is non-dimensionalized by λ) as functions of Sc under the conditions of Re λ from 20 to 150. Figure 4.11(a) indicates that a clear -1 scaling law exists between λ 2 θ and Sc for the largest Re λ case, with a prefactor 0.3, implying that in the limit of intense turbulence the Da θ number can be approximated as

Da θ 0.3Re λ Da = 3T I /τ r . (4.44)
Here τ r = (γ 1 R 2,eq ) -1 is the typical time of the reaction. Remarkably, this result reveals that the unique a posteriori control parameter that we have identified with Da θ can be considered as the ratio of the largest time scale of the turbulent flow to the typical time scale associated to the chemical process.

Finally, we remark that λ θ does not vary significantly over the different scalar fields R 1 , R 2 , P and the reference passive scalar field T . This is exemplified in Figure 4.11(b) for all the simulations at Re λ = 150 and Sc = 1. The figure shows that any λ θ evaluated on a reactive field is at best 15% different from the reference λ T case. Such difference vanishes for very small (the mixing dominated limit) or very large Da and shows a weak increase trend with the order of the reaction. We can conclude that λ T can here be taken as a convenient approximation of λ θ . The estimations of Da θ in the this chapter are based on such an assumption.

Taylor micro scale of scalars advected by a coarse-grained turbulent flow field

In order to understand better the role of λ θ , we perform a series of simulations where the scalar fields are advected by a coarse-grained, i.e. spatially filtered turbulent flow, denoted as ũ. The filter functions as a spectral low pass, defined as

û = |k| K û(k), (4.45)
where K specifies the maximum wave number kept in the modified field. Such a filter retains only the large eddies of the turbulent flow, down to a wavelength 2π/K. The Taylor scale for scalars convected by the filtered flow is denoted as λθ . It is noteworthy that the length quantities are always non-dimensionalized by the Taylor scale of unfiltered flow λ, instead of the Taylor scale of the filtered flow. Interestingly, it is worth exploring what is the impact of the hierarchy of flow scales, extending from the domain size down to the dissipative Reactions in homogeneous isotropic turbulence 6F scales, on the reactive scalar dynamics. In particular, we aim to understand the dependence of the scalar correlation coefficients and the scalar Taylor microscale on the maximum wavenumber K. According to the above discussion, it is reasonably expected that the small scales of the fluid have a negligible influence, because the scaling mixing process and relevant quantities are controlled by the large eddy turnover time T I .

2 2 0.3 Sc Re Re Re Re (a) Da T R1 T Q R2 T Q P T Q R1 T Q R2 T Q P T Q R1 T Q R2 T Q P T Q (b)
Results are presented for a typical case, i.e. Sc = 1, Re λ = 80 and Da = 0.05 corresponding to Da θ = 1.42. 4.12(a) shows that the correlation coefficients vary quite weakly with the filter parameter K. The deviation becomes noticeable only for K ≤ 3, which corresponds to scales larger than the large eddy turnover scale of the flow, in the sense that the forcing is active up to |k| = 2 √ 2. Such behavior is also well captured by the theoretical predictions of (4.33a) and (4.33b) if the Damköhler number Da θ adopted is built on the measured scalar Taylor scale λθ , instead of the original Da θ . This confirms again the relevance of the scalar Taylor microscale in characterizing the present reactive scalar system. It is worth noting that, differently from the correlation coefficients, the scale λθ varies sensibly with the filter wavenumber K, as reported in figure 4.12(b).

-Summary

D r(R 1 , R 2 ) r(R 1 , P) r(R 2 , P) . E Da 3 + 2Da Da 3 + 2Da 2

Summary

In summary, the statistical properties of species undergoing reversible chemical reactions in a turbulent environment have been studied. We have addressed this by means of a model system in which the flow is statistically steady turbulence and the chemical species are retained in a dynamical equilibrium state due to the action of random large-scale source terms. It is observed that the reactive scalar fluctuations have a Gaussian distribution and energy spectra are essentially identical to the one of a passive scalar field transported by the same flow.

This can be explained by the overall small amplitude of the reaction terms in the present close-to-equilibrium conditions. However, in such a state a competition still exists between the chemical processes, which tend to dump reactant concentration fluctuations and enhance their correlation intensity, and the turbulent mixing, which on the contrary increases fluctuations and remove relative correlations.

We quantitatively describe this phenomenon by considering the linearised equations for the reactive scalar fluctuations. A unique control parameter, the Damköhler number (Da θ ), can be constructed as the ratio between the time scale of scalar diffusion across a domain of the size of the scalar Taylor micro-scale (λ 2 θ /D) and the chemical reaction time scale τ r . Importantly, Da θ characterises the functional dependence of fluctuations and correlations of the scalar quantities in the full range of explored conditions with variable reaction order, the

Reactions in homogeneous isotropic turbulence

Reynolds number and the Schmidt number. The larger is such a Damköhler number the more depleted are the scalar fluctuations as compared to the fluctuation of a passive scalar field in the same conditions, and the more intense are the correlations. A saturation in this behaviour is observed beyond Da θ O(10). These results reveal the significance of the scalar Taylor micro-scale for problems involving the mixing of chemical species.

We have shown that in the limit of intense turbulence the relation proposed by [START_REF] Corrsin | Simple theory of an idealized turbulent mixer[END_REF]

λ θ ∼ λSc -1/2
holds approximately, meaning that Da θ can also be viewed as the ratio of the large-eddy-turnover time of the flow over the typical chemical reaction time.

Chapter 5

Turbulent Kolmogorov flows and chemical reactions

In Chapter 4, the statistical behaviors of the reactive scalars in homogeneous isotropic turbulence were investigated. However, realistic turbulent flows are rarely homogeneous or isotropic, because of the presence of boundaries or spatial dependent forces, such as the channel flow and the wall bounded flow. The Kolmogorov flow [START_REF] Meshalkin | Investigation of the stability of a stationary solution of a system of equations for the plane movement of an incompressible viscous liquid[END_REF][START_REF] Green | Two-dimensional turbulence near the viscous limit[END_REF][START_REF] She | Metastability and vortex pairing in the kolmogorov flow[END_REF][START_REF] Borue | Numerical study of three-dimensional kolmogorov flow at high reynolds numbers[END_REF][START_REF] Lucas | Spatiotemporal dynamics in two-dimensional kolmogorov flow over large domains[END_REF] is an example of open turbulent flow, i.e. a flow without boundaries, which is statistically non-homogeneous along one direction and anisotropic. In the late 1950s, A. N. Kolmogorov has proposed to few of his students to study the stability properties of such flow, and an answer was proposed soon after [START_REF] Meshalkin | Investigation of the stability of a stationary solution of a system of equations for the plane movement of an incompressible viscous liquid[END_REF], with a Reynolds number threshold of √ 2 confirmed also later by [START_REF] Green | Two-dimensional turbulence near the viscous limit[END_REF]. The KF system hence corresponds to the Navier-Stokes equations studied in a periodic domain, with a constant pressure, and forced by a sinusoidal forcing. Above the critical Reynolds number, the flow becomes turbulent and we denote this as turbulent Kolmogorov flow (TKF), a flow which is statistically stationary and anisotropic in one direction.

In this chapter, the TKF and the reactive scalars convected in it are discussed based on direct numerical simulations and theoretical analyses. This chapter is divided into two parts. First the classical turbulence closure based on eddy-viscosity Boussinesq's approach and a nonlinear quadratic closure about the velocity fields in the TKF are numerical examined. Then numerical studies about reactive scalars were carried out for extending the theoretical models about the reactive scalars in chemically quasi-equilibrium state discussed in Chapter 4 to reactions in the TKF.

Numerical examinations about closure models

The TKF is now considered a classical flow model system to study non-homogenous, anisotropic, sheared turbulence [START_REF] Borue | Numerical study of three-dimensional kolmogorov flow at high reynolds numbers[END_REF][START_REF] Shebalin | Kolmogorov flow in three dimensions[END_REF][START_REF] Biferale | Anisotropic homogeneous turbulence: hierarchy and intermittency of scaling exponents in the anisotropic sectors[END_REF][START_REF] Balmforth | Stratified kolmogorov flow[END_REF][START_REF] Bibliography Boffetta | The viscoelastic kolmogorov flow: eddy viscosity and linear stability[END_REF][START_REF] Musacchio | Turbulent channel without boundaries: The periodic kolmogorov flow[END_REF]. As discussed by [START_REF] Musacchio | Turbulent channel without boundaries: The periodic kolmogorov flow[END_REF], TKF can be considered, to some respect, as a turbulent channel flow without boundaries. Many of the numerical studies devoted to TKF have focused on two-dimensional configurations [START_REF] She | Metastability and vortex pairing in the kolmogorov flow[END_REF][START_REF] Berti | Elastic waves and transition to elastic turbulence in a two-dimensional viscoelastic kolmogorov flow[END_REF]Lucas andKerswell, 2014, 2015). In this section we consider a three-dimensional Navier-Stokes incompressible turbulence forced in the x direction by a sinusoidal force depending on the z coordinate. *

The Kolmogorov flow model system

Equations of motion and numerical implementations

To generate the TKF, the forcing in the governing equations for the incompressible Navier-Stokes equations (Eq. ( 4.2)) is a sinusoidal function in the form:

f = A sin 2π z H e x , (5.1)
where A is a constant, H is the length of the side of the cubic domain chosen here as the characteristic length scale. Such force, in the x direction and depending on the z variable, makes the flow anisotropic in the z direction. It is convenient to introduce the following reference scales for velocity and time for this section:

U 0 = (AH) 1/2 , (5.2)
T 0 = H U 0 = H A 1/2 . ( 5.3)
From this one can construct the Reynolds number as

Re 0 = U 0 H ν , (5.4)
which thus becomes the only dimensionless control parameter in the system.

In the following discussion in this section, the quantities are in dimensionless units. As in section 4.2, the model system is numerically simulated in a cubic tri-periodic domain, using the means of a pseudo-spectral code Table 5.1: The dimensionless key global parameters after reaching a statistically stationary state. A = 1 is the amplitude of forcing to velocity; H = 1 is the length of the domain, also reference scales for length; U 0 = 1 and T 0 = 1 are the reference scales for velocity and time respectively as indicated in Eq. (5.2). The first column is the Run number; Re 0 = HU0 ν is the Reynolds number based on domain length; Re λ0 = λ0u 0 ν is the Reynolds number based on global Taylor microscale, where λ 0 = 15ν/ 0 u 0 , ν is the kinematic viscosity,

u 0 = 1 3 i u 2 i is the global root-mean square of single component velocity, 0 = ν 2 i j (∂ i u j + ∂ j u i ) 2
is the global energy dissipation rate, • denotes the average in time and all over the domain. N 3 is the grid size; η 0 = (ν 3 / 0 ) 1/4 is the global Kolmogorov scale; |k| max • η 0 is the resolution condition; T total is the total simulation time and T l is the large eddy turnover time, i.e. T total /T l denotes the number of large eddy turnover time spanned by the simulation.

No.

Re 0 Re λ0 ν with a smooth dealiasing technique. The time-marching scheme adopts a third order Runge-Kutta method.

0 η 0 N 3 |k| max • η 0 T total /
The global non-dimensional values of the key parameters for the simulations are reported in table 5.1.

The total integration time is long enough, as shown in Table 5.1, not only to ensure that the statistics have converged, but also to have a sample space large enough for the statistical estimations, especially for the runs with small number of grids. It is noteworthy that the adimensionalization implemented in this section is different from Chapter 4. Because of anisotropy along z direction, the quantities such as the Taylor micro scale λ and Kolmogorov length η are functions of position z in the TKF. The quantities listed in table 5.1 with the subscript of 0 means they are the global values obtained from the average also over z, i.e. • .

Reynolds decomposition and velocity moments

Let us consider a Reynolds decomposition of the velocity into mean and fluctuating quantities u = u z + u ( • z denotes the average over time and spatially in x and y directions) and note the three components of the velocity u = (u, v, w). Because of the periodicity in x and y directions, the derivatives with respect to x and y of mean quantities are 0. By taking the average in time and spatially in x and y direction of the Navier-Stokes equations, one obtains the following relations:

-∂ z τ = 1 Re 0 ∂ 2 z U (z) + sin(2πz), ∂ z w 2 z = -∂ z p z , (5.5)
where U = u z and the shear stress is τ = -u w z . Using the first line of this relation, it is found that for laminar flows when τ = 0, the mean velocity profile is also sinusoidal while the pressure field is constant [START_REF] Meshalkin | Investigation of the stability of a stationary solution of a system of equations for the plane movement of an incompressible viscous liquid[END_REF].

For turbulent flows, it is well-known that the mean velocity profile is also sinusoidal [START_REF] Borue | Numerical study of three-dimensional kolmogorov flow at high reynolds numbers[END_REF], but this is a numerical result which has, to our knowledge, no direct analytical explanation. We obtain the following z-dependence for U :

U (z) = κ sin(2πz).
(5.6)

The numerical estimations of κ are plotted in figure 5.1. The maximum value of the mean turbulent velocity is of the order of the characteristic velocity built using the forcing values, since we obtain values of κ between 1.01 and 1.12, increasing with the Reynolds number (figure 5.1 and Table 5.2). The values of κ found here are compatible with the value of κ = 1.1 found by [START_REF] Borue | Numerical study of three-dimensional kolmogorov flow at high reynolds numbers[END_REF] (the value of the Reynolds number in the latter work is not provided).

For large Reynolds number, using equations (5.5) and ( 5.6) we find that ∂ z τ is proportional to sin(2πz), obtaining finally:

τ = 1 2π 1 -(2π) 2 κ
Re 0 cos(2πz).

(5.7)

The different first and second moments of the velocity obtained after averaging Navier-Stokes equations are shown in figures 5.2 and 5.3. As proposed by [START_REF] Sarris | Box-size dependence and breaking of translational invariance in the velocity statistics computed from three-dimensional turbulent kolmogorov flows[END_REF], one criterion for the convergence of Kolmogorov flow simulations is that the left-hand-side and right-hand-side of equation ( 5.7) must be comparable. As shown in figure 5.2(b), this criterion is satisfied in our simulation (when the Re 0 is large enough for the diffusion term to be neglected). Moveover, only one component of the mean velocity is non-zero; concerning second moments, only one shear stress term is non-zero, involving the product u w . The turbulence is anisotropic since all normal stress components of the stress tensor are different:

u 2 z > w 2 z > v 2 z .
The diagonal terms have twice the spatial frequency of the forcing. Since cos(2θ) = 2 cos 2 θ -1, they can be written as: 5.8) where numerically the couple of values (α i , β i ) of each component of the velocity, for each run, are listed in table 5.2. This will be an important information for the quadratic closure done in the next section. we can write also the evolution of the mean kinetic energy:

u 2 z = α 1 + β 1 cos 2 (2πz), v 2 z = α 2 + β 2 cos 2 (2πz), w 2 z = α 3 + β 3 cos 2 (2πz), (
K(z) = α + β cos 2 (2πz), (5.9)
where the coefficients α and β are plotted in figure 5.1. The values found for the latter are in good agreement with the values reported in [START_REF] Borue | Numerical study of three-dimensional kolmogorov flow at high reynolds numbers[END_REF] (α = 0.34 and β = 0.14). We address again here that these values are Reynolds-number dependent. Since the Reynolds number of the simulation in [START_REF] Borue | Numerical study of three-dimensional kolmogorov flow at high reynolds numbers[END_REF] is not known, quantitative results cannot be precisely compared here.

Closures for the turbulent Kolmogorov flow

Introduce the Reynolds stress tensor defined as T = -u i u j z . The anisotropic stress tensor is R = -T + 2 3 KI, where K is the kinematic energy and I is the unitary tensor. The mean velocity gradient tensor A = ∂U i /∂x j , and the mean strain-rate S and rotation-rate W tensors are also introduced as: (5.11) Closure of the turbulence equations corresponds to express the Reynolds stress tensor using mean quantities, e.g. when the closure is local, using the tensors S and W. Below we first consider the simplest linear closure and estimate the eddy-viscosity, and then we address a nonlinear closure using a quadratic constitutive equation.

S = 1 2 ∂U i ∂x j + ∂U j ∂x i , (5.10)
W = A -S.

Boussinesq's eddy-viscosity hypothesis and its assessment

It is seen from equations (5.6) and (5.7) that, concerning the only non-zero non-diagonal term in the stress tensor, its z-behaviour is the same as the mean gradient term, giving an eddy-viscosity of the form: z . The z-dependence is given by the fits of equation (5.8).

ν T = τ U (z) = 1 (2π) 2 κ - 1 Re 0 . (5.12) z QRUPDOVWUHVVHV u 2 z v 2 z w 2 z (a) Re λ = 38.69 z QRUPDOVWUHVVHV u 2 z v 2 z w 2 z (b) Re λ = 49.31 z QRUPDOVWUHVVHV u 2 z v 2 z w 2 z (c) Re λ = 66.9 z QRUPDOVWUHVVHV u 2 z v 2 z w 2 z (d) Re λ = 122.85
The eddy-viscosity does not depend on z, but depend on the Reynolds number and the coefficient κ. For Run 7 we find a value of ν T = 0.0224, as an overall estimation of the eddy-viscosity with a constant value which does not depend on z. However, the estimation of an eddy-viscosity does not validate the linear closure. The Boussinesq's hypothesis, which is at the basis of all eddy-viscosity turbulence models, corresponding to a linear proportionality between tensors (Boussinesq, 1877) : R = 2ν T S.

(5.13)

For the flow considered here, there are some symmetries so that the strain as well as the stress have a simplified form: and

S = a 2       0 0 1 0 0 0 1 0 0       (5.
R =       2 3 K -σ 2 u 0 τ 0 2 3 K -σ 2 v 0 τ 0 2 3 K -σ 2 w       , (5.15)
where a = U (z), σ 2 u = u 2 z and the same for σ 2 v and σ 2 w .

It is then clear, as is also the case for turbulent channel flows [START_REF] Speziale | On nonlinear kl and k-ε models of turbulence[END_REF][START_REF] Nisizima | Turbulent channel and couette flows using an anisotropic k-epsilon model[END_REF][START_REF] Pope | Turbulent Flows[END_REF], that such linear relation between tensors can be realized only when diagonal terms are zero, i.e.

in an isotropic situation. However, the TKF is anisotropic and as seen in figure 5.3, the three normal stresses are all different, which means that a precise proportionality cannot be found.

In order to quantify this alignment, we consider the inner product between tensors:

A : B = {A t B} = A ij B ij ,
where {X} is a notation for the trace of X. The norm is then ||A|| 2 = A : A. As a direct test of Boussinesq's hypothesis, we first represent here the normalized inner product of R and S tensors (which is 5.1 -Numerical examinations about closure models 91 similar to the cosine of an "angle" between vectors, see [START_REF] Schmitt | Experimental study of the constitutive equation for an axisymmetric complex turbulent flow[END_REF]Schmitt, 2007a)):

ρ RS = R : S ||R|| ||S|| . (5.16)
The ratio ρ RS is thus a number between -1 and 1, which characterizes the validity of Boussinesq's hypothesis:

it is 1 when this hypothesis is valid, and when close to 0 it corresponds to "perpendicular" tensors. This is shown in figure 5.5. It is seen that a plateau is obtained, and that this hypothesis is approximately valid when the velocity gradient is not small, whereas it completely fails for some range of values around the positions where the mean velocity gradient vanishes.

More precisely, from Run 7, we find ρ RS = 0.93 for z = 1/2. And by putting a threshold at ρ RS = 0.9, we find that 0.9 ≤ ρ RS ≤ 1 for z ∈ [0, 0.13] ∪ [0.39, 0.59] ∪ [0.87, 1]. Hence it is larger than 0.9 for approximately 46% of the volume considered: for about half the volume, the linear relation between strain and stress tensors is approximately valid. The two remaining tensors of the tensor basis are: 5.22) and

T 2 = a 2 2       -1 0 0 0 0 0 0 0 1       (
T 3 = a 2 12       1 0 0 0 -2 0 0 0 1      
.

(5.23)

The quadratic constitutive equation finally writes, replacing invariants in equation ( 5.18):

R = 2τ a S + σ 2 u -σ 2 w 1 a 2 T 2 + 6σ 2 v -4K 1 a 2 T 3 .
(5.24) Equation ( 5.24) is a quadratic constitutive equation which is expressing a nonlinear closure of the turbulent Kolmogorov flow. The first coefficient is twice the eddy-viscosity and therefore a constant, as discussed above.

Whereas the other coefficients involve ratio of [constant + cos 2 ] divided by cos 2 terms.

When a = U (z) 0, for z 1/4 and z 3/4, cos(2πz) = 0 and all S, T 2 and T 3 vanish, but in the threeterms development of R, the second term and the third are non-zero constants, since the coefficients diverge (the a 2 terms cancel). In those positions, we see that R is a diagonal tensor which is not vanishing (see figure 5.6);

we see numerically that the second term is also very small and that the third term is dominant. This means that in those positions, Boussinesq's linear eddy-viscosity development is no more valid and the anisotropic stress tensor is a constant perpendicular to the linear term and approximately proportional to T 3 = S 2 -1 3 η 1 I.

Summary

In this section, the classical turbulence closure based on eddy-viscosity Boussinesq's approach and a nonlinear quadratic closure about the velocity fields in the turbulent Kolmogorov flow (TKF) are numerical examined. It was found that the mean velocity profile has the same form, with a damping of a factor κ, with respect to the mean velocity value calculated from the forcing term. The value of this parameter was found to weakly increase with the Reynolds number with indications of a possible asymptotic saturation at very large Re. anisotropic. It was also shown that a quadratic nonlinear constitutive equation can be proposed for this flow, involving a linear term and two nonlinear terms in the form of traceless and symmetric tensors. For about half of the flow domain, the linear term is dominating. Whereas for the vanishing mean velocity gradient regions, only one non-linear term remains non-zero and becomes constant. Hence an effective viscosity coefficient can indeed be estimated for TKF, but contrary to what has been stated previously [START_REF] Rollin | Variations on kolmogorov flow: turbulent energy dissipation and mean flow profiles[END_REF], globally all linear and nonlinear terms are needed for the complete closure. Some of the values obtained here are in agreement with a previous work [START_REF] Borue | Numerical study of three-dimensional kolmogorov flow at high reynolds numbers[END_REF], even if the Reynolds number was not provided in their work. We have checked the Reynolds-number dependence of the different parameters considered here, with 8 different runs with different grid sizes from 128 3 to 512 3 , and with Reynolds numbers from Re λ = 38.7

to 198. The parameters are found to converge for the largest Reynolds numbers considered.

Reactions in turbulent Kolmogorov flow

In Chapter 4, the reactive scalar mixing in homogeneous isotropic incompressible turbulence was studied by direct numerical simulation and theoretical analysis. The fluctuations and correlations of the reactive scalars were the main focus, and found to be uniquely dependant on the control parameter of Da θ . However, according to the definition of Da θ (Eq. ( 4.15)), this number is based on scalar Taylor micro-scale, diffusivity and reaction rate, no the velocity fields. In addition, in the modeling analysis, the fundamental condition is that the reacting system is constrained to a quasi-equilibrium state, weakly relevant to the turbulent flow environment. Therefore, this lead us to consider extending our theoretical analysis about the statistical properties of turbulent reactive scalars in Chapter 4 to non-homogeneous and anisotropic turbulence, for example the TKF discussed above.

In this section we consider the same reactions and passive scalar as that discussed in section 4.1, sustained by the same isotropic constant amplitude forcing (Eq. (4.10)). Specifically, only the second order reaction (n=1) is discussed here. The global non-dimensional values of the key parameters for the simulations are reported in table 5.3. It is noteworthy that the adimensionalization implemented in this section is in the same way as Chapter 4, instead of section 5.1. Similarly with table 5.1, the quantities of the Taylor micro scale λ, Kolmogorov length η and dissipation rate are functions of position z in TKF. The subscript of 0 means λ 0 , η 0 and 0 are the global values obtained from the average also over z, i.e. • .

Scalar Taylor micro-scale

As the Taylor microscale of θ, λ θ is a key quantity in the definition of Da θ (Eq. ( 4.15)). In this chapter, because the flow is supposed to be homogeneous in y direction, λ θ is quantified as TKF, on the Da θ can be determined as

λ 2 θ = θ 2 z (∂ y θ) 2 z . ( 5
z/H R 1 2 z R 2 2 z P 2 z T 2 z (a) z/H R 1 2 z T 2 z R 2 2 z T 2 z P 2 z T 2 z (b) Da R 1 2 z T 2 z R 2 2 z T 2 z P 2 z T 2 z 7KHRU\ R 1 2 z T 2 z R 2 2 z T 2 z P 2 z T 2 z (c)
R 2 1 z T 2 z ≈ R 2 1 z T 2 z ≈ R 2 1 z T 2 z ≈ 3 + 2Da θ 3 + 3Da θ .
(5.28)

As shown in figure 5.9(c), the normalized fluctuations of the reactive scalars measured from the DNS are in agreement with the theoretical predictions by Eq. (5.28).

Summary

In this section, the reversible reaction discussed in Chapter 4 was introduced into the TKF. The Taylor microscales of reactive scalars as functions of the position were found to fluctuate only slightly. Thus parameter of Da θ (the ratio of scalar diffusion time scale to the reaction time scale) is weakly influenced by the anisotropy of the velocity fields. The statistical properties including the fluctuations and correlations of the reactive scalars are numerical examined, whose dependence on Da θ were found to be in good agreement with the theoretical predictions proposed by the modeling analysis discussed in Chapter 4. This indicates that, in the quasi-equilibrium state, the relations between the statistical properties of reactive scalars and the Da θ is weakly influenced by the background velocity field.

Chapter 6

Chemical reactions sustained by Dirichlet boundary conditions

In Chapter 4, the theoretical model was based on the linearization of reaction term, which is valid only when the fluctuations of the scalars are not comparable to the mean quantities. The chemical source was not strong enough to change the scaling behavior of scalars. The main reason is that the reacting system is not strongly deviated from the global equilibrium state in the homogeneous isotropic configuration. In this chapter, a combustor-like configuration of reactive scalars, in which the scalar fields are of Dirichlet boundary conditions in one direction, is explored. The fluctuations of the scalars are maintained by an intrinsic mean gradient instead of isotropic forcing. Accordingly, the reacting system is supposed to be strongly deviated from the chemical equilibrium state and of rich physical interest. The statistical and scaling properties of reactive scalars are investigated, in comparison with a passive scalar. The rest of this chapter is organized as follows. In section 6.1, the definitions of the configuration and problem studied are provided. Section 6.2 elaborates the details of the numerical methods; then in section 6.3, the modeling analysis and numerical results are presented and discussed. Finally the conclusions of this chapter are summarized in section 6.4. *

Introduction of the model system

In this chapter, The velocity fields are exactly the same as that described in section 4.1. While the scalars are different. The reactions considered here are second order reversible reaction, in the following form

R 1 + R 2 γ 1 G GGGGG B F GGGGG G γ 2 P, (6.1)
where R 1 , R 2 and P denote the three involved reactive scalars. The process is reversible with the respective forward and backward reaction rate coefficients γ 1 and γ 2 . The reactants are assumed to undergo diffusion and to be transported in a passive manner by an incompressible velocity flow field upon which they do not exert any effect. The evolution equations scalar concentrations (R 1 (x, t), R 2 (x, t) and P (x, t)) read

∂R 1 ∂t + (u • ∇)R 1 = D R 1 -γ 1 R 1 R 2 + γ 2 P + ṡR1 , (6.2a) ∂R 2 ∂t + (u • ∇)R 2 = D R 2 -γ 1 R 1 R 2 + γ 2 P + ṡR2 , (6.2b)
∂P ∂t + (u • ∇)P = D P + γ 1 R 1 R 2 -γ 2 P + ṡP . (6.2c)
Here u(x, t) is the flow velocity and D is the species diffusivity (assumed identical for all species). The source terms ṡ are precised later below.

In the following analyses, to gain primary insights of the flow physics, a non-reactive species T is also considered for comparison with the following governing equation

∂T ∂t + (u • ∇)T = D T + ṡT . (6.3)
In Chapter 4, it was shown that large-scale statistically homogeneous and isotropic reactive scalar sources/sinks are not able to sustain strong deviation from the chemical equilibrium, which instead needs to be realized by imposing non-zero mean gradients profiles for the reactants. The canonical homogeneous shear flow, although with adjustable mean scalar gradients, can not be adopted because the scalar concentrations defining the reaction rates in Eqs. (6.2) are undetermined. Moreover, other flow cases such as the shear layer [START_REF] Mellado | Gradient trajectory analysis of a scalar field with external intermittency[END_REF] are not suitable either because of the unstationary evolution or inherent spatial non-homogeneity.

Thus we propose here a new flow configuration, which is schematically illustrated in figure 6.1. In a cubic domain a large-scale forcing term f (Eq. (4.9)) is exerted into the momentum equation For the scalars quantities, the periodic boundary conditions are set along x and y directions, while a Dirichlet boundary condition is used along the z direction. The shadowed layers near the boundaries are the "buffer layers" generated artificially, in which the quantities of scalars are close to the preset boundary values, as defined in Eq. (6.5). The part between buffer layers is denote as the bulk region. Such setup is statistically stationary and ensures the local positiveness of scalar concentrations. fields are periodic only in x and y directions. Along the z direction, the following Dirichlet boundary conditions are implemented:

       R 1 = R 0 , R 2 = 0, P = 0, T = R 0 when z = 0, R 1 = 0, R 2 = R 0 , P = 0, T = 0 when z = H, (6.4)
where H is the length of domain in z direction and R 0 is the constant boundary condition.

Numerically it is found that to realize a reasonably large fluctuation for the scalar fields, buffer layers in the vicinity of the Dirichlet boundaries are needed, as shown in figure 6.1 by the shadowed parts with the bulk region in between. Inside both the upper and bottom buffer layers, the artificial sources ṡθ with θ = R 1 , R 2 , P or T are added in the scalar equations Eq. (6.2) and Eq. (6.3). Specifically, ṡθ is designed here as

ṡθ =        1 τ (θ 0 -θ)
, in the buffer with z > H -δ or z < δ, 0, in the bulk.

(6.5)

Here θ 0 is the boundary value of θ, as defined in Eq. (6.4); τ is a characteristic time scale to control the strength Chemical reactions sustained by Dirichlet boundary conditions of the source terms. Obviously, small values of τ imply a large source to reduce the defect of the scalar quantities from the boundary values. In our present simulation cases, τ is set at the order of 100 times of the numerical time step. Another parameter δ is the buffer layer thickness, which can be tailored to adjust the scalar source (i.e. lager thickness corresponds to larger species source) and the scalar mean gradient.

Numerical implementation

First of all, it is convenient to non-dimensionalize the above set of equations by choosing reference scales appropriate for the present system. The domain size H, the overall scalar difference R 0 in Eq. (6.4) and the overall fluctuating velocity u are used as the reference quantities for the length scale, scalar and velocity, respectively. For simplicity, in the following, symbols by default denote the corresponding nondimensionalized ones. It then yields The reaction rates computed with the mean quantities as functions of z. The solid lines are for the forward reaction rates Da 1 R 1 R 2 z and the dashed lines are for the backward reaction rates Da 2 P z . The clear difference can be observed. Vertical dotted lines mark the interfaces between the buffer layers and the bulk region. The key involved parameters in the present simulations are listed in table 6.1.

∂R 1 ∂t + (u • ∇)R 1 = (Sc Re) -1 R 1 -Da 1 R 1 R 2 + Da 2 P + ṡR1 , (6.6a) ∂R 2 ∂t + (u • ∇)R 2 = (Sc Re) -1 R 2 -Da 1 R 1 R 2 + Da 2 P + ṡR2 , (6.6b

Re

The same as described in section 4.2, the turbulent flow is statistically stationary homogeneous and isotropic, sustained by a large-scale forcing term of (4.9). The isotropic velocity field is obtained by numerically solved by using a pseudo-spectral code [START_REF] Gauding | High-order structure functions for passive scalar fed by a mean gradient[END_REF][START_REF] Gauding | One-point and two-point statistics of homogeneous isotropic decaying turbulence with variable viscosity[END_REF] with a smooth dealiasing technique [START_REF] Hou | Computing nearly singular solutions using pseudo-spectral methods[END_REF] for the treatment of non-linear terms in the equations. Differently, the scalar equations Eq. (6.6) are solved by the finite difference method.

The velocity field is initialized by prescribing the spectrum of kinematic energy in the Fourier space û(k; 0) [START_REF] Schumacher | Asymptotic exponents from low-reynolds-number flows[END_REF], where both the modulus and phases are randomly determined, under the constraint Meanwhile, the scalars are linearly initialized as:

                       R 1 (x, y, z; 0) = 1 -z, R 2 (x, y, z; 0) = z, P (x, y, z; 0) = 0, T (x, y, z; 0) = 1 -z.
(6.8) Figure 6.3 shows the typical visualization of R 1 and the reaction rate (Da 1 R 1 R 2 -Da 2 P ) on the isosurface of R 1 = 0.5 under the condition of Γ = 10. It can be seen that with the properly defined buffer layer, the spatial fluctuation of scalar can be self sustained with the total scalar quantities confined in the prescribed [0, 1] range.

The statistical stationarity of the setup can be appreciated from the temporal evolution of the scalar statistical moments. As an example, figure 6.4 shows that the evolutions of the spatial average of scalar concentrations and root mean square of scalar fluctuation in the bulk region for the case Γ = 10. This indicates that the reacting system is strongly deviated from chemical equilibrium. The following analyses will be focused on the bulk region. Data samples are collected in a time span of about ten times of the integral time T I = k/ with k = 3u 2 /2, once the statistically stationary state is reached.

6.3 -Result analyses 105 t/(k/ ) PHDQ R 1 z R 2 z P z T z (a) mean t/(k/ ) UPV R 1 2 1/2 z R 2 2 1/2 z P 2 1/2 z T 2 1/2 z (b) fluctuation

Result analyses

Properties of the buffer layer and the bulk region

In the present configuration, the two buffer layers function as a source which sustain the mean scalar gradient with positive scalar concentrations, while the bulk domain is where the turbulent mixing occurs and where the present analysis is focused on. The means and root-mean-square (r.m.s.) of the fluctuations of the passive scalars (T ) in the configurations with different buffer layer thickness δ are shown in figure 6.5 (a) and (b), respectively. As can be expected, in the bulk region the mean scalar profile follows a linear relation with respect to z, i.e. constant gradient, under the action of isotropic turbulent velocity. From the prescribed geometrical and boundary conditions, the constant gradient is about 1/(1 -2δ). For the scalar fluctuation, figure 6.5(b) demonstrates that in the bulk region the scalar r.m.s. is almost constant, similar to the homogeneous shear turbulence [START_REF] Mellado | Gradient trajectory analysis of a scalar field with external intermittency[END_REF]. In most of the buffer layer, the scalar r.m.s. is negligibly small, because the strong modulation effect from the source term ṡT makes the scalar T roughly constant. Specifically, the strength of such modulation is determined by the control parameter τ . In this sense, the present flow can be effectively tailored by the control parameters τ and δ.

z PHDQ T z = T z = T z = (a) z UPV T 2 1 2 z = T 2 1 2 z = T 2 1 2 z = (b)
Further insights of the relation between the buffer and bulk regions can be gained by the following analytical approach. By taking the ensemble average ( • z ) of Eq. (6.6d), it yields

∂ t T z = (D T + D) T z + ṡT z , (6.9)
where the Reynolds stress term has been absorbed into the turbulent diffusivity D T according to the following relation

uT z = -D T ∂ T z ∂z . (6.10)
Because of symmetry, only a half of the domain, in the range of z ∈ [0, 1/2] needs to be studied. Under the statistical stationary condition, the temporal derivative term in Eq. (6.9) vanishes. Let us denote D T in the bulk region and buffer layer as D T,1 and D T,2 , respectively. Combining the specific form of ṡT (Eq. (6.5)) and 6.3 -Result analyses 107 neglecting the laminar diffusivity D in Eq. (6.9), an ordinary differential equation about T z is then found:

       D T,2 d 2 T z dz 2 = 1 τ ( T z -1) when 0 ≤ z ≤ δ, D T,1 d 2 T z dz 2 = 0 when δ ≤ z ≤ 1/2.
(6.11)

The numerical evolution of D T , calculated according to Eq. (6.10), is shown in figure 6.6 (a). Similar to the r.m.s. profile of T , D T is also negligibly small in most of the buffer layer, following the same modulation mechanism as found for the source ṡT . From the analytical point of view, we approximate D T,1 and D T,2 as zindependent constants. As discussed before, the strong modulation effect from ṡT , or more specifically from the controlling parameter τ , leads to a small value of D T,2 , while physically D T,1 is determined by the flow integral time T I . Thus we further assume that 1)) needs to be determined numerically.

D T ,2 D T ,1 ∼ τ T I , or quantitatively D T ,2 D T ,1 = K τ T I , where the proportionality coefficient K (K = O(
To solve this set of ODE, four boundary conditions are needed, including T z (0) = 1, T z (1/2) = 1/2, the continuity of T z at z = δ, and the continuity of the flux of T z at z = δ. Because of the different diffusivity in the buffer layer' and the bulk region, the continuity of the flux of T z at z = δ can be expressed as

D T,2 d T z dz (z → δ -) = D T,1 d T z dz (z → δ + ). (6.12)
Therefore, the analytical solution of Eq. (6.11) for T z is obtained as

T z (z) =        C 1 sinh( 1 D T ,2 τ z) + 1 when 0 ≤ z ≤ δ, C 2 z + C 3 when δ ≤ z ≤ 1/2.
(6.13)

where

C 1 = -1 2(A -B(δ-0.5) D T ,1 ) , C 2 = BC 1 D T,1 , C 3 = 1 2 - C 2 2 , with A = sinh( 1 τ D T ,2 δ) and B = cosh( 1 τ D T ,2 δ) D T ,2
τ . Numerically, it is found that the model solution with 6.10). (b) Theoretical prediction of the mean of T (dashed lines) compared with the DNS results (solid lines with same color). The vertical dotted lines mark the interfaces between the buffer layers and the bulk region. modulation effect from the source term ṡT (and other ṡ terms as well) in the buffer.

Statistical properties of the reactive scalars

In this section, we focus on the difference between the z-dependent statistical properties of the reactive and passive scalars, or the effects of the chemical reaction. It can be seen from Eq. (6.6) and the corresponding boundary conditions (6.7) that R 1 (z) = R 2 (1 -z) statistically, i.e. R 1 and R 2 are symmetric with respect to the middle of the domain with z = 1/2. Therefore, for the sake of brevity the results for R 2 will be omitted in the rest analyses. Figure 6.7 shows on the middle plane with z = 1/2 the PDFs of the scalar quantities R 1 and P , together with that of the passive scalar T . Overall, the PDF of T , denoted as p T , is symmetric and has a central maximum at T = 0.5. Moreover, since the effect of the buffer layers can be compared to the mixing process in the shear layer, two other local peaks appear at the tails of p T at T = 0 and T = 1 (same for other PDFs) [START_REF] Mellado | Gradient trajectory analysis of a scalar field with external intermittency[END_REF]. Larger Γ lead to stronger skew of p R1 toward the R 1 = 0 side. Such skewness property is the consequence of chemical reactions, because faster forward chemical reactions tend to deplete the reactants R 1 and R 2 but enrich the product P , enhancing p R1 at the R 1 = 0 end and extending the p P toward the larger P side. T . We can see clearly the mirror symmetry between p T (z) and p T (1 -z), which, however, breaks down for the R 1 case, because of the strong influence from the chemical sources. In addition, the PDF of T is of particular importance. In the modelling analysis discussed in the following section, the moments of the reactive scalars can be theoretically predicted based on p T undergoing the same turbulent environment.

-Result analyses

109 R 1 RUT 3') R 1 R 1 R 1 R 1 R 1 T (a)
The dependence of the PDF of the net reaction rate

R net = Da 1 R 1 R 2 -Da 2 P = Da 1 (R 1 R 2 -P/Γ) on Γ
is presented in figure 6.9 (at z = 1/2). Toward the fast chemical limit with large Γ, the PDF peaks higher at the R net = 0 end and meanwhile becomes more extended toward the higher R net side. When Da 1 and Da 2 are comparable, the PDF peaks at some moderate value of R net . Since all the cases are under the control of the identical turbulence velocity, such difference must be caused by the chemical mechanism, which can be more clearly viewed from the spatial distribution of the reaction rates. From the comparison between figure 6.10 (a) and (d), there is a clear difference between the distribution of R net for Γ = 100 and Γ = 1. For large Γ (and large Da 1 as well), the large R net regions are highly concentrated in thin stripes, while for small Γ, regions with high R net are much broader distributed, which explains the local bump in the PDF profile in figure 6.9. More detailed understanding of such property can be clarified from the separated results of the forward and backward reaction rates. It can be seen that a similar difference appears between figure 6.10 (b) and (e), while, figure 6.10 (c) and (f) are weakly influenced or even uninfluenced by Γ. Because the forward reaction rate is determined

Da 1 R 1 R 2 Da 2 P 3') Figure 6.9: PDF of the reaction rate (Da 1 R 1 R 2 -Da 2 P ) under the conditions of different Γ, at the position of z = 1/2.
z is slightly greater than 1/2, the gradient of R 1 z is equal to that of T z , resulting in equality of r.m.s. of R 1 and r.m.s of T . For the product P , its r.m.s reaches a maximum at the edge of the bulk region and a minimum in the middle (z = 1/2). A tentative explanation is that the fluctuation of P is determined by the fluctuation of R 1 and R 2 because of the chemical kinetics and flow dynamics, i.e. turbulent transport. At the edge of the bulk, the r.m.s. of R 1 is large, but R 2 fluctuates weakly, which can not lead to high peak of the r.m.s. of P .

Therefore, such maximum must come from the contribution from turbulent transport, or specifically, the large gradient of P z close to the bulk edge (see the P z results), due to the gradient hypothesis. In parallel, at z = 1/2 the gradient of P z and the turbulent transport part vanish, leading to the minimum of the r.m.s. of

P at z = 1/2.
To have further understanding of the effects of chemical reaction on the scalar statistics, the present reactive turbulent system needs to be analyzed theoretically. Let us define X = R 1 -R 2 . Subtracting Eq. (6.6a) by Eq. (6.6b), it yields (6.15)

∂ t X + (u • ∇)X = (Sc Re) -1 X + ṡR1 -ṡR2 , (6.14) (a) (b) (c) (d) (e) (f)
Comparing the governing equation and boundary conditions of X with that of the passive scalar T (Eq. (6.6d) and (6.7)), it yields

X = R 1 -R 2 = 2T -1. (6.16)
Let us define p X (x; z) as the pdf of X at the position of z (similar definition for other random quantities).

Relation (6.16) gives:

p X (x; z) = 1 2 p T ( x + 1 2 ; z). (6.17)
First consider the case of infinitely large Da 1 . This implies that R 1 and R 2 cannot coexist because of the finite

chemical source Da 1 R 1 (x, t)R 2 (x, t), leading to R 1 (x, t)R 2 (x, t) = 0. Therefore, a positive X(x, t) is equivalent to R 1 (x, t) = X(x, t) and R 2 (x, t) = 0, while X(x, t) < 0 implies R 1 (x, t) = 0 and R 2 (x, t) = -X(x, t).
For the quantity P , subtracting Eq. (6.6d) to Eq. (6.6a), together with the boundary conditions, we conclude that P (x, t) = T (x, t) -R 1 (x, t). (6.18) Therefore, a relation between P (x, t) and X(x, t) can be obtained directly. Here, the infinitely large Da 1 leads to the following relations:

       R 1 (x, t) = 2T (x, t) -1, R 2 (x, t) = 0, P (x, t) = 1 -T (x, t), when T (x, t) 1/2, R 1 (x, t) = 0, R 2 (x, t) = 1 -2T (x, t), P (x, t) = T (x, t), when T (x, t) < 1/2. (6.19)
Consequently, given that the passive scalar field T is known, the mean and r.m.s. of R 1 (at infinite Γ) can be respectively determined as

R 1 z (z) = 1 0 xp X (x; z)dx = 1 0 1 2 xp T ( x + 1 2 ; z)dx =2 T |T > 1 2 z (z) - 1 1 2 p T (t; z)dt, (6.20) 6.3 -Result analyses 115 and R 2 1 z (z) = R 2 1 z (z) -R 1 2 z (z) = 1 0 1 2 x 2 p T ( x + 1 2 ; z)dx -[ 1 0 1 2 xp T ( x + 1 2 ; z)dx] 2 . (6.21)
The similar derivation can be done for P . The predictions are shown in figure 6.12 and 6.13.

For the finite but large Γ, R 1 (x, t) and R 2 (x, t) can locally coexist, i.e. R 1 R 2 > 0. Since the overall forward reaction is still strong (if Γ is sufficiently larger than unity), we assume here that there exits a upper limit for

R 1 (x, t)R 2 (x, t), i.e. R 1 (x, t)R 2 (x, t) ≤ C Γ , (6.22)
where C is a constant to be determined. Moreover, for any given X

(x, t) ∈ [-1, 1], another constraint is R 1 (x, t)R 2 (x, t) ≤ 1-|X(x, t)|.
The reason can be explained by the fact that the species concentrations R 1 (x, t) and R 2 (x, t) need to be confined in the range of [0, 1]. Thus for a given value of X 6.24) in the case of X(x, t) < 0. In summary, R 1 (x, t)R 2 (x, t) ≤ 1 -|X(x, t)|. Putting these together, it gives

(x, t) = R 1 (x, t) -R 2 (x, t), if R 1 (x, t) ≥ R 2 (x, t) i.e. X(x, t) ≥ 0, we have R 1 (x, t)R 2 (x, t) ≤ R 2 (x, t) = R 1 (x, t) -X(x, t) ≤ 1 -X(x, t). (6.23) Similarly, R 1 (x, t)R 2 (x, t) ≤ R 1 (x, t) = R 2 (x, t) + X(x, t) ≤ 1 + X(x, t), (
R 1 R 2 ∈ [0, min( C Γ , 1-|X(x, t)|)] = [0, β max ]. For a given X(x, t) = α and R 1 (x, t)R 2 (x, t) = β, R 1 = α+ √ α 2 +4β 2 . If the conditional PDF of R 1 R 2 on X, i.e. p R1R2|X (β|α; z), is known, R 1 z as function of z can be determined as R 1 z (z) = 1 -1 p X (α; z) βmax 0 α + α 2 + 4β 2 p R1R2|X (β|α; z)dβdα. (6.25)
A hypothesis assumed here is that with given

X(x, t), R 1 (x, t)R 2 (x, t) is equally distributed in [0, β max (α)].
Together with the numerical results of the PDF of the passive scalar T , the mean R 1 z (z) and variance

R 2 1 z (z) = R 2 1 z (z) -R 1 2
z (z) can then be calculated (similar analyses for R 2 and P ). As shown in figure 6.12 and 6.13, when Γ 10, the modeling and numerical results can satisfactorily match if the constant C in .12: The scalar mean: (a) R 1 z and (b) P z , as a function of z obtained from theoretical analysis (dashed lines) based on Eq. ( 6.25) and DNS (solid lines with the same colors). The grey dashed lines are from the theoretical prediction at infinitely large Da 1 according to Eq. (6.20). We see that the prediction for Γ = ∞ is close to the curves for Γ = 100 and also the predictions for large Γ are close to the DNS results when Γ=10, 30, 100. The vertical dotted lines mark the interfaces between the buffer layers and the bulk region.

Eq. (6.22) is set as 0.7. When Γ < 10, these predictions do not hold.

Correlation Coefficients

For scalars θ 1 and θ 2 under consideration, the correlation coefficients are defined (based on the fluctuating parts) as

r(θ 1 , θ 2 )(z) = θ 1 θ 2 z θ 2 1 1/2 z θ 2 2 1/2 z . (6.26)
The scalar correlation is jointly determined by the chemical reaction and the turbulent mixing. In the present flow configuration, the numerical simulations of the z-dependent correlation coefficients are shown in figure 6.14.

In Chapter 4, it was found that a competition exists between the chemical reaction and turbulent mixing.

Specifically, the chemical reaction tends to dump reactant concentration fluctuations and enhance their correlation intensity, while turbulent mixing increases fluctuations and removes relative correlations. For the present non-equilibrium configuration, the reactants R 1 and R 2 consume each other, especially when the forward reaction is strong, which explains the negative r(R 1 , R 2 ) in figure 6.14(a). In the bulk region, R 1 and R 2 are less becomes abnormally more negative. This abnormal property can be explained by the fact that at sufficiently large Γ the chemical source Da 1 R 1 R 2 in the buffer layer is much weaker, because either R 1 or R 2 is negligibly small. Therefore, r(R 1 , R 2 ) is close to be that of the non-reactive case.

z r(R 1 , R 2 ) (a) z r(R 1 , P) (b)
For r(R 1 , P ), since R 1 and P function as the mutual sources rather than sink, the result is different from r(R 1 , R 2 ). Overall, r(R 1 , P ) increases from -1.0 at z = 0.0 to 1.0 at z = 1.0. At z = 0 the reaction rate of P is mainly determined by R 2 since R 1 remains close to constant as 1.0. Because of the stoichiometric relation, the defect of R 1 from 1.0 is determined by either R 2 or P . Therefore, r(R 1 , P ) ∼ -1.0. In a similar manner, at z = 1 the reaction rate or the generation rate of P is mainly determined by R 1 and thus r(R 1 , P ) ∼ 1.0. For the nonreactive case, at the middle plane with z = 0.5, R 1 z and R 2 z are exactly equal. Thus, the concentration of P is not influencing either R 1 or R 2 , yielding r(R 1 , P ) ∼ 0. With increasing Γ, R 1 will be more consumed and P will be more produced. As shown in figure 6.7, larger Γ leads to the PDF of R 1 more skewed toward the R 1 = 0 side, while the PDF of P skews differently toward the large P side. Therefore, r(R 1 , P ) will unanimously decrease and shift downwards with increasing Γ, as demonstrated in figure 6.14(b).

Scalar Energy Spectra

The z-dependent scalar energy spectra is also investigated. At a specific z, the energy spectra corresponding to a two-dimensional scalar field is defined as (6.27) where k is the 2D wave number and k = |k|, • k denotes the average in time, θ(k) is the Fourier coefficients of the mode of k, θ * (k) is the corresponding complex conjugate.

E θ (k, z) = 2πk 2 θ(k) θ * (k) k , θ = R 1 , R 2 , P or T,
In Chapter 4, it was shown that the scalar energy spectra at the quasi-equilibrium state with different chemical sources are almost identical, because of the negligibly small reaction rates. For the present nonequilibrium reactive turbulence cases, the chemical source plays important roles in determining the structure and statistics of the scalar quantities. It is interesting to consider if such chemical source will largely influence or change the energy spectra. Consider the R 1 case for instance in figure 6.15 (a)-(c), higher Γ makes more scalar energy to shift from small wave number range to the large wave number range, indicating that stronger chemical reactions tend to lump the local scalar quantity and strengthen the scalar intermittency. A similar tendency appears for the scalar P , as shown in figure 6.15 (d)-(f). Such Γ effect becomes stronger at z = 1/2

(more difference between the curves in figure 6.15 (e)), since P z reaches maximum at z = 1/2.

In addition, the coherency spectrum between two scalars θ 1 and θ 2 is defined as: especially when the reaction is strong. In figure 6.16, the spectrum peaks toward the high wave number end, indicating that the correlation between R 1 and R 2 is mainly from the small scale contribution, in consistence with the stripe like structures visible in figure 6.10.

Co θ1,θ2 (k) = | θ 1 (k) θ 2 * (k)| k θ 1 (k) θ 1 * (k) k θ 2 (k) θ 2 * (k) k , ( 6

Summary

To maintain the definite positiveness of the species concentration, an original flow configuration has been proposed, where the species are supplied from the buffer boundaries with adjustable thickness to drive chemical reactions at the strongly out of equilibrium state. This allows to go beyond our previous exploration of reactive scalar fields in homogeneous turbulence, where only a moderate out of equilibrium state could be attained, and their results could be explained in term of a linearization approach of the reaction rates (Chapter 4).

Here, a theoretical modelization that take as an input the PDF of a passive scalar field can satisfactorily predict the local mean and fluctuation of reactants for Γ 10. The correlation coefficient between the scalar quantities are determined by two counteracting effects, the turbulent mixing and the chemical kinetics. For larger Γ, the spatial distribution of the forward reaction and net reaction takes a stripe like structure, making the scalar field more intermittent.

In comparison with the existing results of the isotropic turbulence case as discussed in Chapter 4, we conclude that the scalar correlation is jointly determined by the chemical source and physical conditions, e.g., flow configuration and boundary setup. Under the non-equilibrium condition with strong chemical sources, the chemical reaction also plays important roles in determining the scalar energy spectra. Consistently, large Γ will shift both the scalar energy and coherency spectra from the small wave number range to the large wave number range, which is in a sharp contrast with the near equilibrium case.

Chapter 7

Conclusion and future works

Summary of the thesis

In this thesis, the statistical properties of reactive scalars undergoing reversible chemical reactions in incompressible turbulence were studied by the means of direct numerical simulations and theoretical analyses. The cases studied include the reactions close to and strongly deviated from the chemical equilibrium states, in homogeneous isotropic and non-homogeneous anisotropic turbulent flow. Theoretical analysis on the statistical properties of scalars at different order of moments were carried out based on appropriately proposed approximations and models. The theoretically derived results were then compared with numerical results obtained by direct numerical simulation.

In the direct numerical simulations, for solving the PDEs of the quantities with periodic boundary conditions, the spatial derivatives were approximated by using a pseudo-spectral method with smooth dealiasing. Otherwise, for the special configurations in which the boundary conditions are not periodic, a finite difference method with high-order schemes was used. The statistical analysis of the quantities of interest was implemented in a statistical stationary state, which was maintained by numerical forcing for the velocity and scalars. In the anisotropic configuration, the forcing is the inherent mean gradient. While in the isotropic configuration, the velocity and scalar fluctuations are sustained with constant input power and amplitude respectively.

In the case of reversible chemical reactions in homogeneous isotropic turbulence, the reacting system is found to be in a chemical quasi-equilibrium state. In the homogeneous isotropic case, the theoretical analysis about the statistical properties of the reactive scalars are weakly affected by the small scale properties of the turbulent flow. This inspired us to consider the non-homogeneous anisotropic turbulence, such as the turbulent Kolmogorov flow (TKF). First, we gained a deeper insight into the TKF with a forcing of the sinusoidal form. It was found that the mean velocity profile has the same form, with a damping of a factor κ, with respect to the mean velocity value calculated from the forcing term. The value of this factor was found to weakly increase with the Reynolds number with indications of a possible asymptotic saturation at very large Re. The only non-zero shear stress term is proportional to the cosine function, and the normal stress components all involve a square cosine expression. The normal stresses are never equal, showing that as expected the turbulence is anisotropic. It was also shown that a quadratic nonlinear constitutive equation can be proposed for this flow, involving a linear term and two nonlinear terms in the form of traceless and symmetric tensors. For about half of the flow domain, the linear term is dominating.

Whereas for the vanishing mean velocity gradient regions, only one non-linear term remains non-zero and becomes constant. Hence an effective viscosity coefficient can indeed be estimated for TKF, but contrary to 7.1 -Summary of the thesis 125 what has been stated previously [START_REF] Rollin | Variations on kolmogorov flow: turbulent energy dissipation and mean flow profiles[END_REF], globally all linear and nonlinear terms are needed for the complete closure. Then, the reversible reactions discussed for the homogeneous isotropic case was also studied with the background flow as the TKF. The theoretical predictions about the dependence of the fluctuations and the correlations of the reactive scalars on the Da θ in homogeneous isotropic flow are found to work also for the TKF case. It indicates that, in the quasi-equilibrium state, the dependence of the statistical properties of reactive scalars on the scalar diffusion and reaction rate is weakly influenced by the background velocity field.

For the irreversible reactions in homogeneous isotropic turbulence, the linearization of the reaction term is valid because the fluctuations of the reactive scalars are small when compared to the mean quantities in the quasi-equilibrium state. However, a reacting system strongly deviated from the global equilibrium state and with relatively larger fluctuations of reactive scalars can be more practically meaningful and physically interesting. This is because the reaction flow possesses stronger non-linearities, which is a symbolic feature of turbulence problems. Thus, a combustor-like configuration of reactive scalars, in which the scalar fields are of Dirichlet boundary conditions in one direction, was explored. In such a configuration, the entire flow consists of two buffer layers and a bulk region. The fluctuations of the scalars are maintained by intrinsic mean gradient instead of isotropic forcing. Accordingly, the reacting system is supposed to be far from the chemical equilibrium state. The statistical and scaling properties of reactive scalars are investigated, in comparison with a passive scalar. A theoretical model based on the PDF of the passive scalar can satisfactorily predict the mean and fluctuation of reactants, if the associated Damköhler number (Da 1 ) is sufficiently large. For large Da 1 , the spatial distributions of the forward reaction and net reaction assume stripe like structures, which makes the scalar field more intermittent as compared to the reference passive scalar case. The correlation coefficient between the scalar quantities are determined by two counteracting effects, the turbulent mixing and the chemical kinetics. Under the non-equilibrium condition with strong chemical sources, the chemical kinetics also plays important roles in determining the scalar energy spectra. Consistently, large Da 1 will shift the scalar energy from the small wave number range to the large wave number range, which is different from the near equilibrium case. In addition, the coherency spectra of the reactive scalars are strongly wave number dependent, in sharp distinction with the near equilibrium case.

In summary, reactions in incompressible turbulence of both homogeneous isotropic and non-homogeneous anisotropic cases were explored in this thesis. Theoretical analyses about the statistical properties of the reactive scalars, in comparison with a passive scalar, were implemented and examined with direct numerical simulations for different configurations.

Scientific contributions

The scientific contributions of this thesis can be summarized as following:

• The topic of this thesis, which is a combination of reaction and passive scalar in turbulence, is novel;

• We developed a configuration for the numerical simulation of turbulent reaction, in which the positiveness and strong fluctuation of reactive scalars are ensured;

• In the study about reactions in homogeneous isotropic turbulence, we introduced a key parameter (Da θ ) characterizing the competition between scalar diffusion and reaction, and showed its important role in determining the fluctuations and correlations of reactive scalar;

• In study about Chemical reactions sustained by Dirichlet boundary conditions, we proposed modeling analysis building a link between the moments of reactive scalar and the PDF of non-reactive scalar.

Future works

Compared with studies about turbulent combustion and passive scalars in turbulence, relatively few works have been devoted to the reaction in incompressible turbulence, especially reversible reactions. One unavoidable issue concerning reversible reactions is the chemical equilibrium. The statistical properties of reactive scalars can be hugely different in a reacting system close to and strongly deviated from the chemical equilibrium state. For the case with quasi-equilibrium state discussed in this thesis, the fluctuations of scalars have to be constrained as about 10% of the mean quantities for ensuring the positiveness of the concentrations. However, such a phenomenon does not make good practical sense for typical turbulence problems. In the real ocean, the fluctuations displayed by planktonic populations can attain values that are comparable to the ones of the mean population density, because there exist regions where nearly no individual or where huge accumulation can be observed. Simple solutions can be to enforce the numerical negative concentrations as zero or to express the concentration quantities in the form of power functions. However, these operations introduce too much arbitrary interference into the problem. It is thus important to propose configurations ensuring chemical quasiequilibrium, the positiveness of reactive scalars and not small scalar fluctuations, preferably without too much artificial interference or at least only with physically reasonable corrections to the concentration quantities.

A straightforward way to amplify the fluctuations is by deriving the system into a non-equilibrium state.

Another reason that the non-equilibrium state is of interest is that the reaction term can be more important.

One operation offered in this thesis is by introducing Dirichlet boundary conditions for the reactive scalars (as 7.3 -Future works 127 in Chapter 6). However, such a configuration breaks the isotropy of the scalars. Further consideration about a configuration retaining the isotropy of the scalar fields and meanwhile possessing a reaction source non-negligible compared with the convection or diffusion can be interesting. Moreover, it is also interesting to investigate the irreversible reactions, in which the reaction term can be non-negligible because the chemical equilibrium is not involved.

Based on the homogeneous isotropic cubic domain investigated in Chapter 4, an anisotropic configuration with Dirichlet boundary conditions and buffer layers for scalars was developed and discussed in Chapter 6. The reactive scalar fields in the latter configuration are much more complex, because the non-linearity induced by the reaction source is well reserved. Thus, the modeling analyses about the statistical properties of reactive scalars is much more challenging. Such a configuration is of great potential for deeper and broader exploration beyond the theoretical analysis carried out so far in this thesis. In the present works, the theoretical predictions about the mean profiles and fluctuations of the reactive scalars could match the numerical tests only for the cases with forward reaction infinitely large or much larger than the backward reaction. In addition, the current theoretical predictions about the mean profiles and fluctuations of the reactive scalars requires the prerequisite that the PDF of the passive scalar undergoing the same convection and diffusion is known. In future works on such a configuration, an important focus can be the improvement of the modeling analysis for a wider application range and proposing more compact expressions about the moments of the reactive scalars. Another direction can be the modification of the current hypotheses, which work poorly for the weak forward reaction case. Or any new theoretical frameworks based directly on the governing equations of the reactive scalars can also be considered.

Furthermore, more in-depth and comprehensive analysis about other statistical properties of reactive scalars, for example the correlation coefficients, can also be a interesting subject for future works.

Concerning the Kolmogorov flow case studied in this thesis, in future works we suggest also to explore different shapes of forcing, to better understand the expressions of the different moments of the velocity field. It remains also to be understood from analytical arguments why the eddy-viscosity does not depend on z for such flow, contrary to what is found in similar but different flow such as channel flow [START_REF] Schmitt | Direct test of a nonlinear constitutive equation for simple turbulent shear flows using dns data[END_REF] or boundarylayer flows. Moreover, the closure models for passive scalar in the Kolmogorov flow can also be examined, and the extension to the reactive scalars will be more interesting.

The scaling behavior of turbulent scalars has always been of high interest since the Kolmogorov-Obukhov-Corrsin theoretical framework. The reaction, especially fast reaction, is supposed to introduce cliff structure into the scalar fields, and thus to promote the intermittency. For examining this idea and further exploring, comparison of the relevant statistics, such as the structure functions and the energy spectra, between passive 128 Conclusion and future works and reactive scalars can be useful. Since the reaction term are commonly non-linear, the convolution of the existing scalar fields can be reflected in high wave number motions. Thus, for a convincing direct numerical simulation about the turbulent problems concerning reaction, the demand in the computational resource is higher, since the numerical spatial resolution needs be finer than that required for passive scalar.

In one word, this thesis carried out basic and preliminary studies about the reactive scalars in incompressible turbulence by direct numerical simulations and theoretical analyses. The new discoveries always accompany with new problems. As the complement or development of the works done in this thesis, a lot of interesting questions are still waiting for further investigations in future works.

  une combinaison des paramètres importants caractérisant le mouvement turbulent et la diffusion scalaire, y compris le nombre de Reynold basé sur Taylor (Re λ ), le nombre de Schmidt (Sc), le nombre de Damköhler traditionnel (Da) et l'échelle de Taylor pour les scalaires. (λ θ ). Les prédictions théoriques sur les coefficients de corrélation montrent qu'avec un petit Da θ (réaction faible), les scalaires réactifs sont faiblement corrélés avec des coefficients de corrélation négligeables. Au fur et à mesure que Da θ augmente (réaction plus forte), les scalaires réactifs deviennent progressivement soit positifs soit négatifs en corrélation avec les coefficients de corrélation évoluant asymptotiquement vers certaines constantes. Les transitions entre zéro et les quantités saturées des coefficients de corrélation se produisent toujours lorsque Da θ est à peu près dans la plage de [0,1,10]. En ce qui concerne l'amplitude de la fluctuation scalaire, bien que la réaction montre un effet global sur le dumping des fluctuations scalaires, nous pouvons supposer que l'apport d'énergie aux scalaires réactifs pour lesquels l'ordre de la réaction est de 1 (en R 1 et P) est à peu près le même qu'avec l'apport d'énergie au scalaire non réactif. Avec cette hypothèse supplémentaire, la modélisation peut être étendue aux fluctuations des scalaires réactifs normalisés par la fluctuation du scalaire non réactif, qui se sont également avérés être des fonctions du paramètre de contrôle de Da θ . Les fluctuations normalisées des scalaires réactifs étaient théoriquement comme étant 1 lorsque Da θ est petit (réaction faible), indiquant que les scalaires réactifs sont des mêmes fluctuations avec le scalaire non réactif subissant le même environnement turbulent et le même forçage externe. Lorsque Da θ augmente (réaction plus forte) jusqu'à des valeurs de l'ordre de 0.1, les fluctuations normalisées commencent à diminuer. Les tendances à la baisse cessent lorsque Da θ atteint environ 10, les fluctuations normalisées saturant à certaines constantes. Les transitions se sont avérées se produire à peu près dans la même plage de Da θ que celle des coefficients de corrélation. Selon la définition de Da θ , la compétition entre les processus chimiques et le mélange turbulent est essentielle pour déterminer les fluctuations de concentration des réactifs et leur intensité de corrélation. Plus précisément, les processus chimiques ont tendance à réduire les fluctuations de concentration des réactifs et à améliorer leur intensité de corrélation. Alors qu'au contraire, le mélange turbulent augmente les fluctuations et supprime les corrélations relatives. Nos analyses de modélisation théorique ont été fortement validées par les simulations numériques, car les dépendances fonctionnelles des fluctuations et les corrélations des grandeurs scalaires caractérisées par Da θ sont bien confirmées par les résultats DNS dans la gamme complète des conditions explorées avec un ordre de réaction variable (2 à 4), le nombre de Reynolds basé sur Taylor (20 à 150) et le nombre de Schmidt (0.1 à 4).

  linéaire entre les tenseurs de déformation et de contrainte est approximativement valide. Il a également été montré qu'une équation constitutive non linéaire quadratique peut être proposée pour ce flux, impliquant un terme linéaire et deux termes non linéaires sous forme de tenseurs sans trace et symétriques. Pour environ la moitié du domaine d'écoulement, le terme linéaire domine. Alors que pour les régions de gradient de vitesse moyenne de fuite, un seul terme non linéaire reste non nul et devient constant. Par conséquent, un coefficient de viscosité efficace peut en effet être estimé pour TKF, mais contrairement à ce qui a été dit précédemment, globalement tous les termes linéaires et non linéaires sont nécessaires pour la fermeture complète. Après l'analyse du champ de vitesse dans le TKF, les réactions réversibles discutées pour le cas isotrope homogène ont également été introduites pour étendre notre analyse théorique sur les propriétés statistiques des scalaires réactifs turbulents à la turbulence non homogène et anisotrope, par exemple la TKF. Ici, nous n'avons considéré qu'une réaction du second ordre avec un nombre de Schmidt égal à 1. En tant que grandeur clé dans la définition de Da θ , le carré de l'échelle micrométrique scalaire de Taylor λ θ s'est avéré être une fonction cosinus carrée de la position z. Cependant, ils n'ont pas de grandes fluctuations, correspondant à moins de 10% de la valeur moyenne. Ainsi, λ θ peut grossièrement être considéré comme constant dans TKF. Pour les fluctuations normalisées et les coefficients de corrélation des scalaires réactifs, qui sont les principaux sujets des prédictions théoriques proposées dans le cas isotrope homogène précédent, nos résultats numériques ont montré qu'ils sont faiblement influencés par la position en TKF. Plus important encore, les dépendances des fluctuations normalisées et des coefficients de corrélation des scalaires réactifs sur le Da θ , théoriquement dérivés dans le flux isotrope homogène, pourraient également être validées par les résultats numériques dans le cas TKF actuel. Ceci indique que, dans l'état de quasi-équilibre, les dépendances des propriétés statistiques des scalaires réactifs sur la diffusion scalaire et la vitesse de réaction sont faiblement influencées par le champ de vitesse sous-jacent. Pour les réactions irréversibles en turbulence isotrope homogène, la linéarisation du terme de réaction est valable car les fluctuations des scalaires réactifs sont faibles par rapport aux quantités moyennes à l'état de quasi-équilibre. Cependant, un système réactif fortement dévié de l'état d'équilibre global et avec des fluctuations relativement plus importantes des scalaires réactifs peut être plus significatif sur le plan pratique. Par exemple dans l'océan, les fluctuations affichées par les populations planctoniques peuvent atteindre des valeurs comparables à celles de la densité moyenne de population. De plus, les scalaires réactifs avec des fluctuations plus importantes sont plus intéressants physiquement, car l'écoulement réactif possède une non-linéarité plus forte, ce qui est une caractéristique principale des écoulements turbulents. Pour ces raisons, il sera intéressant d'explorer différentes configurations d'écoulement où de fortes sources chimiques pourraient produire des écarts prononcés par rapport à l'état d'équilibre global. Ainsi, une configuration semblable à une chambre de combustion de scalaires réactifs, dans laquelle les champs scalaires possèdent des conditions aux limites de Dirichlet dans une direction, a été explorée. Dans une telle configuration, l'écoulement global se compose de deux couches tampons générées par une technique numérique et d'une région centrale de mélange entre les deux. Les fluctuations des scalaires sont maintenues par un gradient moyen intrinsèque imposé par les conditions aux limites, au lieu d'un forçage isotrope comme adopté auparavant. En conséquence, le système réactif est supposé être fortement dévié de l'état d'équilibre chimique. Les propriétés statistiques et d'échelle des scalaires réactifs dans un système de réaction réversible du second ordre avec un nombre de Schmidt de 1 ont été étudiées, en comparaison avec un scalaire passif non réactif. Dans ce cas, il existe deux nombres de Damköhler: Da 1 et Da 2 pour les réactions directe et inverse respectivement. Nous avons maintenu Da 2 comme une constante et avons fait varier le rapport Γ = Da 1 /Da 2 de 1 à 100.

  de réaction. Les visualisations des taux de réaction ont montré que, pour les petits Γ, les régions avec une réaction directe élevée sont distribuées sur des valeurs étendues. Alors que pour Γ grand, les distributions spatiales de la réaction directe et de la réaction total se distribuent selon des structures en forme de bande, ce qui rend le champ scalaire plus intermittent par rapport au scalaire passif de référence. L'intermittence peut également être étudiée via les spectres d'énergie scalaire bidimensionnels. À partir des spectres d'énergie des scalaires réactifs normalisés par les spectres d'énergie du scalaire passif, il a été constaté qu'un Γ plus élevé fait passer plus d'énergie scalaire de grandes échelles de longueur à de petites échelles, indiquant que les réactions chimiques plus fortes ont tendance à concentrer la quantité scalaire locale et à renforcer l'intermittence scalaire. Les coefficients de corrélation entre les scalaires réactifs ont été déterminés par deux effets antagonistes, le mélange turbulent et la cinétique chimique. Fait intéressant, les réactifs R 1 et R 2 se sont révélés plus positivement corrélés avec un Γ plus grand, bien qu'ils soient censés se consommer plus rapidement. Ce résultat peut s'expliquer par le fait que les PDF de R 1 et R 2 étaient de plus en plus positivement biaisés avec l'augmentation de Γ dans la configuration actuelle. De plus, les spectres de cohérence des scalaires réactifs étaient fortement dépendants du nombre d'onde, en nette distinction avec le cas de quasi-équilibre. En particulier, la corrélation entre les scalaires réactifs construits avec une forte réaction provient principalement de la contribution à petite échelle. A l'instar des analyses faites dans le cas homogène et isotrope, des prédictions théoriques sur les moments des scalaires réactifs (r.m.s des fluctuations et les grandeurs moyennes) ont également été réalisées ici. Une condition préalable nécessaire est que le PDF du scalaire passif (T) soit connu. Essentiellement, les analyses de modélisation actuelles se concentrent sur le lien entre le scalaire réactif et le scalaire passif, par lequel les moments des scalaires réactifs peuvent être calculés à partir du PDF connu du scalaire passif. Pour le cas de Γ infiniment grand (réaction directe extrêmement forte), les réactifs R 1 et R 2 ne peuvent pas coexister. Dans cette condition, des relations définies entre les scalaires réactifs et le scalaire passif peuvent être obtenues. Alors que pour le cas de Γ fini, nous avons dû proposer une hypothèse que le produit de R 1 et R 2 (R 1 R 2 ) a localement une limite supérieure, qui est une constante divisée par Γ. De plus, la PDF conditionnelle de R 1 R 2 sur le scalaire passif doit être modélisée. Ici, celle-ci a été initialement considérée comme étant uniformément distribuée. Ainsi, les prédictions théoriques sur les fluctuations moyennes et les quantités moyennes des scalaires réactifs pourraient être obtenues en faisant l'intégrale du scalaire réactif dans l'espace d'échantillonnage statistique du scalaire passif. Lorsque nous avons validé les résultats de la modélisation avec des simulations numériques, il a été constaté que les prédictions théoriques pouvaient correspondre aux résultats numériques lorsque Γ est significativement grand. Les modèles devraient être encore améliorés pour comprendre les situations avec un Γ plus petit.

1 Introduction 1 . 1

 111 ) nous avons développé une configuration pour la simulation numérique de la réaction turbulente, dans laquelle la positivité et la forte fluctuation des scalaires réactifs sont assurées; (3) dans l'étude sur les réactions en turbulence isotrope homogène, nous avons introduit un paramètre clé (Da θ ) caractérisant la compétition entre diffusion scalaire et réaction, et montré son rôle important dans la détermination des fluctuations et corrélations du scalaire réactif; (4) dans l'étude sur les réactions chimiques soutenues par les conditions aux limites de Dirichlet, nous avons proposé une modélisation établissant un lien entre les moments du scalaire réactif et le PDF du scalaire non réactif. Travaux futurs Par rapport aux études sur la combustion turbulente et les scalaires passifs en turbulence, relativement peu de travaux ont été consacrés à la réaction en turbulence incompressible, en particulier les réactions réversibles. Un problème inévitable concernant les réactions réversibles est l'équilibre chimique. Les propriétés statistiques des scalaires réactifs peuvent être extrêmement différentes dans un système réactif proche et fortement dévié de l'état d'équilibre chimique. Pour le cas avec état de quasi-équilibre discuté dans cette thèse, les fluctuations des scalaires doivent être contraintes à environ 10% des quantités moyennes pour assurer la positivité des concentrations. Cependant, un tel phénomène n'a pas de bon sens pratique pour des problèmes de turbulence typiques. Dans l'océan réel, les fluctuations affichées par les populations phytoplanctoniques peuvent atteindre des valeurs comparables à celles de la densité moyenne de population, car il existe des régions où la densité est très faible et des situations où une énorme accumulation peut être observée. Des solutions simples peuvent être d'imposer des concentrations négatives numériques à zéro ou d'exprimer les quantités de concentration sous la forme de fonctions de puissance. Cependant, ces opérations introduisent trop d'interférences arbitraires dans le problème. Il est donc important de proposer des configurations assurant le quasi-équilibre chimique, la positivité des scalaires réactifs et non de petites fluctuations scalaires, de préférence sans trop d'interférences artificielles ou du moins seulement avec des corrections physiquement raisonnables des quantités de concentration. Une façon simple d'amplifier les fluctuations consiste à faire passer le système dans un état de non-équilibre. Une autre raison pour laquelle l'état de non-équilibre est intéressant est que le terme de réaction peut être plus important. Une opération proposée dans cette thèse consiste à introduire des conditions aux limites de Dirichlet pour les scalaires réactifs. Cependant, une telle configuration brise l'isotropie des scalaires. Une réflexion plus approfondie sur une configuration conservant l'isotropie des champs scalaires et possédant en même temps une source de réaction non négligeable par rapport à la convection ou à la diffusion peut être intéressante. De plus, il est également intéressant d'étudier les réactions irréversibles, dans lesquelles le terme de réaction peut être non négligeable car l'équilibre chimique n'est pas impliqué. Sur la base du domaine cubique isotrope homogène étudié dans un premier temps dans cette thèse, une configuration anisotrope avec des conditions aux limites de Dirichlet et des couches tampons pour les scalaires a été développée et discutée. Les champs scalaires réactifs dans cette dernière configuration sont beaucoup plus complexes, car la nonlinéarité induite par la source de réaction est bien présente et active. Ainsi, les analyses de modélisation sur les propriétés statistiques des scalaires réactifs sont beaucoup plus difficiles. Une telle configuration est d'un grand potentiel pour une exploration plus profonde et plus large au-delà de l'analyse théorique menée jusqu'à présent dans cette thèse. Dans les présents travaux, les prédictions théoriques sur les profils moyens et les fluctuations des scalaires réactifs pourraient correspondre aux tests numériques uniquement pour les cas avec une réaction directe infiniment grande ou beaucoup plus grande que la réaction inverse. De plus, les prédictions théoriques actuelles sur les profils moyens et les fluctuations des scalaires réactifs nécessitent la condition préalable que le PDF du scalaire passif subissant la même convection et diffusion soit connu. Dans les travaux futurs sur une telle configuration, un objectif important peut être l'amélioration de la modélisation pour une gamme d'applications plus large et la proposition d'expressions plus compactes sur les moments des scalaires réactifs. Une autre direction peut être la modification des hypothèses actuelles, qui fonctionnent mal pour le cas de réaction directe faible. Ou tout nouveau cadre théorique basé directement sur les équations gouvernantes des scalaires réactifs peut également être envisagé. En outre, une analyse plus approfondie et plus complète des autres propriétés statistiques des scalaires réactifs, par exemple les coefficients de corrélation, peut également être les sujets intéressants pour les travaux futurs sur cette configuration. Concernant le cas d'écoulement de Kolmogorov étudié dans cette thèse, dans les travaux futurs, nous proposons également d'explorer différentes formes de forçage pour mieux comprendre les expressions des différents moments du champ de vitesse. Il reste également à comprendre à partir d'arguments analytiques pourquoi la viscosité turbulente est une constante indépendante de la position dans un tel écoulement, contrairement à ce que l'on trouve dans un écoulement similaire mais différent tel qu'un écoulement en canal ou des écoulements en couche limite. De plus, les modèles de fermeture du scalaire passif dans le « Kolmogorov flow » peuvent également être examinés, et l'extension aux scalaires réactifs sera plus intéressante. Les lois d'échelle des scalaires turbulents a toujours été d'un grand intérêt depuis le cadre théorique de Kolmogorov-Obukhov-Corrsin. La réaction, en particulier la réaction rapide, est censée introduire une structure « de falaise » (« ramp-cliff ») dans les champs scalaires, et ainsi favoriser l'intermittence. Pour examiner cette idée et l'explorer plus avant, la comparaison des statistiques pertinentes, telles que les fonctions de structure et les spectres d'énergie, entre les scalaires passifs et réactifs peut être utile. Etant donné que le terme de réaction est généralement non linéaire, la convolution des champs scalaires existants peut se refléter dans les nombres d'onde élevés. Ainsi, pour une simulation numérique directe convaincante sur les problèmes turbulents concernant la réaction, la demande en ressource de calcul est plus élevée, car la résolution spatiale numérique doit être plus fine que celle requise par le scalaire passif. En un mot, cette thèse a mené des études de base et préliminaires sur les scalaires réactifs en turbulence incompressible, par des simulations numériques directes et des analyses théoriques. Les nouvelles découvertes accompagnent toujours de nouveaux problèmes. En complément ou en développement des travaux réalisés dans cette thèse, Vitae . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0.1.1 PhD Thesis: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0.1.2 Education: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0.1.3 Participation to summer schools and conferences: . . . . . . . . . . . . . . . . . . . . . . . 0.1.4 Papers published and in preparation: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0.2 Objectives of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0.3 Organization of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Turbulent flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Navier-Stokes equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.2 Energy cascade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.3 Kolmogorov's 1941 theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Passive Scalars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  because the intermittency of the scalar field is relatively stronger. . . . . . . . . . . . . . . . . . . 1.4 A visualization of a passive scalar field, by Chen and Kraichnan (1998). . . . . . . . . . . . . . . 1.5 The scaling exponents ζ and ζ θ from experiments by[START_REF] Schmitt | Relating lagrangian passive scalar scaling exponents to eulerian scaling exponents in turbulence[END_REF][START_REF] Schmitt | Linking eulerian and lagrangian structure functions' scaling exponents in turbulence[END_REF], compared with the theoretical predictions of K41 and KOC theories[START_REF] Huang | Arbitrary Order Hilbert Spectral Analysis Definition and Application to Fully Developed Turbulence and Environmental Time Series[END_REF]. . . . . . . . . . . . . . . . . . . 2.1 Scatter plots about the example θ 1 -θ 2 joint distribution in the case of (a) r(θ 1 , θ 2 ) = -0.9; (b) r(θ 1 , θ 2 ) = 0; and (c) r(θ 1 , θ 2 ) = 0.9. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 (a) Figure 1 in Corrsin (1961): A quantitative sketch of the variance spectrum of reactive scalar undergoing a first order reaction. (b) Figure 1 in Corrsin (1964a): A quantitative sketch of the variance spectrum of reactive scalar undergoing a second order reaction. Here C is the coefficient of reaction rate; k or m is the wave number; G is the wave number dependant scalar variance. . . 2.3 Figure 2(b) in Molemaker and de Arellano (1998): one-dimensional spectrum of the covariance between two chemical species (averaged in space and time and normalized). . . . . . . . . . . . . 3.1 The comparison of the energy spectra of velocity numerically obtained by pseudo-spectral methods with and without dealiasing. k = |k| is the amplitude of wave number vector. The dotted curve is the filtering function of smooth dealiasing. . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 The illustration of the physical domain of dependence of Eq. (3.31) (thick dashed line) and the numerical domain of dependence (gray area). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii LIST OF FIGURES

  function (Eq. (3.6)) shown as the dotted curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.4 An illustration for the balance between reaction and diffusion in the reaction zone. . . . . . . . . 57 4.1 Evolution of the root mean square of scalar fluctuations and mean values for the case of n = 2, Da = 0.1, Sc = 1 and Re λ = 150. Time is normalized by the integral time k/ε with k = 3u 2 /2. The mean quantities were represented with three-dimensional volumetric average here. The dash vertical line marks the initial time for the computation of statistical quantities. . . . . . . . . . 63 4.2 Energy spectra of reactive and passive scalar fields, i.e. E θ (k) with θ = R 1 , R 2 , P, T , and velocity field E u (k) (solid light-blue line) in the condition Re λ = 150, Sc = 1, n = 1 and Da = 10, 1, 0.1, 0.01, 0.001 (from high to bottom). Each spectra is compensated with the KOC scaling, (kη) 5/3 and normalized by the global energy θ 2 . For clarity, the energy spectra of scalars are shifted vertically by a multiplicative factor 0.1 . . . . . . . . . . . . . . . . . . . . . . 64 4.3 Coherency spectra of the reactive scalars, under the condition of Re λ = 80, Sc = 1, n = 1. The horizontal axis is normalized in terms of the Kolmogorov scale η. The dash vertical line marks the maximum wave number at which the scalar source terms acts. . . . . . . . . . . . . . . . . . 65 4.4 The global correlation coefficients of the reactive scalars and their gradients along x direction, under the condition of Re λ = 150, Sc = 1, n = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.5 Correlation coefficients between R 1 and R 2 (r(R 1 , R 2 )), R 1 and P (r(R 1 , P )), R 2 and P (r(R 2 , P )) as functions of Da θ , under the condition of (a): n = 1 (b): n = 2 and (c): n = 3 and the Schmidt number Sc = 1. Theoretical predictions are shown in black lines. . . . . . . . . . . . . . . . . . . 72 4.6 Correlation coefficients between R 1 and R 2 (r(R 1 , R 2 )), R 1 and P (r(R 1 , P )), R 2 and P (r(R 2 , P )) as functions of Da θ , under the condition of (a): Re λ = 20 (b): Re λ = 40 and (c) Re λ = 80. The order of R 2 (n) is 1. Theoretical predictions are shown in black lines. . . . . . . . . . . . . . . . . 73 4.7 PDF of the scalar quantities (main panel) and normalized with respect to their standard deviations (inset, the black lines are Gaussian curves) under the condition of Re λ = 150, Sc = 1, n = 2 and Da = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.8 Visualisation of two-dimensional instantaneous sections of the R 1 field at Re λ = 150, Sc = 1 and for three different Da values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.9 PDF of reaction rate Da(R 1 R n 2 -P ) at Da = 0.1 and Sc = 1 for n = 1 (left) and n = 2 (right) for different Reynolds numbers Re λ = 20, . . . , 150. The top panels shows the raw PDF while the bottom ones show the PDF of the normalised variable with respect to its standard deviation. The two insets display the trend of the variance of the reaction rate as a function of the Reynolds number in log-log scale, the red line has a slope 3. . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.10 The fluctuations of the reactive scalars normalized by the fluctuation of passive scalar (T ) as functions of Da θ , under the condition of (a): n = 1 (b): n = 2 and (c): n = 3 and the Schmidt number Sc = 1. Theoretical predictions are shown in black lines. . . . . . . . . . . . . . . . . . . 78 LIST OF FIGURES ix 4.11 a) Taylor micro scales of scalars (computed on T ) as functions of Sc. The red line draw is 0.3/St. A power law fit of the form aSc b on the Re λ = 150 data set gives a = 0.29 ± 0.1 and b = -0.93 ± 0.02. b) Taylor micro scales of the reactive scalars with respect to the passive scalar one with different order of reaction, as functions of Da for all the simulations at Re λ = 150 and Sc = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.12 (a) Correlation coefficients of scalar evolving in a coarse-grained turbulent flow filed, compared with theoretical predictions and (b) Taylor micro scales of the scalars convected by filtered velocity ( λθ ) as functions of the maximum wave number of filtered velocity (K), at Sc = 1, Re λ = 80 and Da = 0.05 (Da θ = 1.42). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 The coefficients obtained by fitting the profiles of mean velocity and kinetic energy, κ in Eq. 5.6 and α and β in Eq. 5.9, as function of Re λ0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 (a) The adimensional mean quantities of each component of the velocity of Run 8. The only nonzero term is u z , having a maximum value of κ, where κ = 2.82. (b) The different adimensional shear stress terms of Run 8. The only non-zero term is u w z . Its z dependence is given by relation (5.7). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 The different normal stresses: one finds u 2 z > w 2 z > v 2 z . The z-dependence is given by the fits of equation (5.8). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 The mean kinetic energy K(z) = 1 2 u i u i z . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Simulation results for the test of the validity of Boussinesq's hypothesis, representing the alignment ρ RS between R and S. The mean velocity profile is superposed in dotted line for reference. 5.6 The amplitudes of the terms at the right hand side of Eq. (5.24) as function of z. The mean velocity profile is also represented as a dotted line, for reference. The horizontal red dotted lines mark the 0 value for the amplitudes. The simulation results of Run 1, 3, 5, 7 are shown here. . . 5.7 The square of the Taylor microscale of θ under the condition of Da=8.5 compared with that of the velocity, as functions of z. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Correlation coefficients between reactive scalars (a) as function of Da θ at the position of z = H/2; (b) as function of z under the condition of Da=8.5. . . . . . . . . . . . . . . . . . . . . . . . . . 5.9 (a) Variances of reactive scalars, under the condition of Da=8.5, and passive scalar as function of z. (b) Variances of reactive scalars normalized by the variance of passive scalar as function of Da θ at the position of z = H/2. (b) Variances of reactive scalars normalized by the variance of passive scalar as function of z under the condition of Da=8.5. . . . . . . . . . . . . . . . . . . . . 6.1 Schematic diagram of the flow configuration and computational domain. For the scalars quantities, the periodic boundary conditions are set along x and y directions, while a Dirichlet boundary condition is used along the z direction. The shadowed layers near the boundaries are the "buffer layers" generated artificially, in which the quantities of scalars are close to the preset boundary values, as defined in Eq. (6.5). The part between buffer layers is denote as the bulk region. Such setup is statistically stationary and ensures the local positiveness of scalar concentrations. . . . . x LIST OF FIGURES 6.2 The reaction rates computed with the mean quantities as functions of z. The solid lines are for the forward reaction rates Da 1 R 1 R 2 z and the dashed lines are for the backward reaction rates Da 2 P z . The clear difference can be observed. Vertical dotted lines mark the interfaces between the buffer layers and the bulk region. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 6.3 The three-dimensional instantaneous snapshot of (a) R 1 , and (b) reaction rate (Da 1 R 1 R 2 -Da 2 P ) on the isosurface of R 1 = 0.5, under the condition of Γ = 10. . . . . . . . . . . . . . . . . . . . . 104 6.4 Evolution of the spatial average of the mean and root mean square of scalar fluctuation in the bulk region for the case of Γ = 10. Time is normalized by the integral time k/ with k = 3u 2 /2. The dashed vertical lines marks the initial time for the computation of statistical quantities. . . . 105 6.5 Numerical profiles of: (a) the mean, and (b) the r.m.s. of T with different buffer layers thickness δ. The vertical dotted lines mark the interfaces between the buffer layers and the bulk region. . . 106 6.6 (a) Turbulent diffusivity calculated from Eq. (6.10). (b) Theoretical prediction of the mean of T (dashed lines) compared with the DNS results (solid lines with same color). The vertical dotted lines mark the interfaces between the buffer layers and the bulk region. . . . . . . . . . . . . . . . 108 6.7 Dependence of PDFs at z = 1/2 on Γ for: (a) R 1 and T , and (b) P . The insert panel in (b) plots the peak of p P as function of Γ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 6.8 PDF of (a) R 1 under the condition of Γ = 10 and (b) T at different positions. . . . . . . . . . . . 110 6.9 PDF of the reaction rate (Da 1 R 1 R 2 -Da 2 P ) under the conditions of different Γ, at the position of z = 1/2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 6.10 The instantaneous two-dimensional snapshot of reaction rate at the position of z = 1/2. The upper row (a,b,c) correspond to Γ = 100; the lower row (d,e,f) correspond to Γ = 1. The first column (a,d) show the net reaction rate (Da 1 R 1 R 2 -Da 2 P ); the second column (b,e) show the forward reaction rate (Da 1 R 1 R 2 ); the third column (c,f) show the backward reaction rate (Da 2 P ).112

  . . . . . . . . . . . . . . . . . . . . . . . . . . . z obtained from theoretical analysis (dashed lines) and DNS (solid lines with the same colors). The grey dashed lines are from the theoretical prediction at infinitely large Da 1 according to Eq. (6.21). We see that the prediction for Γ = ∞ is close to the curves for Γ = 100 and also the predictions for large Γ are close to the DNS results when Γ=10, 30, 100. The vertical dotted lines mark the interfaces between the buffer layers and the bulk region. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 LIST OF FIGURES xi 6.14 DNS simulations of the correlation coefficients between (a) R 1 and R 2 ; (b) R 1 and P as function of z, for different Γ cases. The vertical dotted lines mark the interfaces between the buffer layers and the bulk region. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.15 The ratio of the energy spectra of (a) R 1 at z = 1 4 ; (b) R 1 at z = 1 2 ; (c) R 1 at z = 3 4 ; (d) P at z = 1 4 ; (e) P at z = 1 2 ; (f) P at z = 3 4 to the energy spectrum of T at the same z. . . . . . . . . . 6.16 2D coherency spectra between R 1 and R 2 at the positions of (a) z = 1/4, (b) z = 1/2 and (c) z = 3/4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4. 1

 1 Parameters for the simulations: Re λ is the Taylor scale based Reynolds number; Sc is the Schmidt number; Da = λγ 1 R n 2,eq /u = λγ 2 /u is the Damköhler number based on Taylor scale; n is the order of R 2 in the reaction; N 3 is the total number of grid points; η is the Kolmogorov length or dissipative length; |k| max is the maximum wave number amplitude kept by the dealiasing procedure; |k| max • η is the spatial resolution condition; dt/τ η is the time step normalized by the Kolmogorov time scale τ η . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 The dimensionless key global parameters after reaching a statistically stationary state. A = 1 is the amplitude of forcing to velocity; H = 1 is the length of the domain, also reference scales for length; U 0 = 1 and T 0 = 1 are the reference scales for velocity and time respectively as indicated in Eq. (5.2). The first column is the Run number; Re 0 = HU0 ν is the Reynolds number based on domain length; Re λ0 = λ0u 0 ν is the Reynolds number based on global Taylor microscale, where

  eddy turnover time spanned by the simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 The numerical values of the coefficients in Eq. 5.8 for each run. . . . . . . . . . . . . . . . . . . . 5.3 Key global parameters after reaching statistically stationary state. Re λ0 = u λ 0 /ν is the global Taylor microscale, where u is the global root mean square of single component velocity, ν is kinematic viscosity, λ 0 is the global Taylor micro scale; Sc is the Schmidt number; n is the reaction order of R 2 in the reaction; N 3 is total number of grids; η 0 = (ν 3 / 0 ) 1/4 is the global Kolmogorov length scale; |k| max is the maximum wave number amplitude kept by the dealiasing procedure, |k| max • η is the resolution condition; Da = λ 0 γ 1 R n 2,eq /u = λ 0 γ 2 /u is the Damköhler number based on Taylor scale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Wu, W., L. Wang, E. Calzavarini, F. G. Schmitt. Reactive scalars in incompressible turbulence with strongly non-equilibrium chemistry. To be sumbitted.

4

  focuses on the reaction processes in homogeneous isotropic turbulence. In Chapter 5, an anisotropic and non-homogeneous turbulent flow: the turbulent Kolmogorov flow is discussed by testing an original turbulence closure model and by investigating chemical reactions in such flow. As a development of the homogeneous isotropic case, a configuration with Dirichlet boundary conditions for scalars is introduced and investigated in Chapter 1

  g. the drawing by Leonardo da Vinci of the flowing water out of a ditch, as shown in figure 1.1.

Figure 1 . 1 :

 11 Figure 1.1: Drawing of turbulence by Leonardo da Vinci (15th century).

  's turbulent cascade. Energycontaining eddies break up into smaller eddies and energy is transferred until it is dissipated by viscosity.

  these substances are important in many circumstances, such as in the atmospheric studies, chemical and thermal engineering, and biological research. They are all one-dimensional quantities and share the same form of governing equation. Thus they are described as scalars (in comparison with vector such as velocity). Scalars in the turbulent environment are transported and dispersed by the turbulent eddies from largest to smallest scales, under the action of thermal or molecular diffusion as well. In many cases, the advected scalars have strong effects on the turbulent flow, i.e. the velocity field. For example, some considerable change in temperature can significantly change the fluid density. These scalars are called active scalars, which are beyond the consideration of this thesis. The passive scalars are undergoing dispersion and mixing processes mainly driven by the turbulent flow, but does not couple back on the flow dynamics, leading to negligible change in the property of background flow

Figure 1 . 3 :

 13 Figure 1.3: The probability density function of normalized velocity (upper curve) and temperature (lower curve) in a grid-generated turbulence (shown in the work of[START_REF] Jayesh | Probability distribution of a passive scalar in grid-generated turbulence[END_REF]). For legibility purpose, the upper curve has been shifted by 2 decades with respect to the lower one. Note that in the labels, the θ and u are the fluctuations and θ and u are the corresponding rms. The deviation of PDF from Gaussian curve for the scalar is larger than that of the velocity, because the intermittency of the scalar field is relatively stronger.

Figure 1 .

 1 4 is a numerical result showing the intermittency in a passive scalar field, in which the large scale of plateau structures are separated by small scale cliff structures.

Figure 1 . 4 :

 14 Figure 1.4: A visualization of a passive scalar field, by Chen and Kraichnan (1998).

  3 grid points was implemented for an isotropic steady turbulence with Re λ = 805. Two passive scalars with Sc = 0.72 were studied: one (θ) is forced by a random source that is Gaussian and white in time, and the other (q) is forced by the mean uniform scalar gradient. The local scaling exponents of the two passive scalars: ζ θ n,0 (r) and ζ q n,0 (r) as functions of separation distance r and order n were introduced. They were computed based on the zeroth (isotropic) terms Introduction of the Legendre expansions of the structure function. Then the scaling exponents ζ θ n and ζ q n were determined as the slopes of the ζ θ n,0 (r) -r and ζ q n,0 (r) -r curves (obtained by fitting with least-square method). Mostly importantly, ζ θ n and ζ q n were found to be very close to each other. Such a result revealed the universality in the scaling behavior of turbulent passive scalar. In figure1.3, non-Gaussian tails are more clearly found in the PDF profiles of passive scalars compared with the PDF of velocity. Moreover, figure1.5 shows the ζ θ -q curve is more bending than the ζ-q curve, i.e. the structure function of the passive scalar is of more abnormal scaling behavior. Both these results suggest that the passive scalar field is more intermittent than the velocity.

Figure 1 . 5 :

 15 Figure 1.5: The scaling exponents ζ and ζ θ from experiments by[START_REF] Schmitt | Relating lagrangian passive scalar scaling exponents to eulerian scaling exponents in turbulence[END_REF][START_REF] Schmitt | Linking eulerian and lagrangian structure functions' scaling exponents in turbulence[END_REF], compared with the theoretical predictions of K41 and KOC theories[START_REF] Huang | Arbitrary Order Hilbert Spectral Analysis Definition and Application to Fully Developed Turbulence and Environmental Time Series[END_REF].

  death-saturation model. For an individual species A, the model considering only birth and death states that the temporal evolution of A follows dA dt = bA -dA. (2.16) Here the b and d are the growth rate and death rate respectively. The terms of bA and -dA indicate that this process consists of the combination of the reaction of A b -→ 2A and A d -→ ∅. In such case, it is easy to find that A evolves exponentially as e (b-d)t . Thus if b is smaller than d, A will decay to infinitesimally small quickly. The case with b larger than d is of more interest, in which A can grow explosively. However this is unrealistic in biology. The problem lies in the modelling about the growth. The growth of a species replies on the supply from the environment, such as food and oxygen. Thus the growth subjects to the finiteness of the resources. The reaction A b -→ 2A has to be replaced with A + B b -→ 2A, where B represents the food. The modified coefficient of reaction rate b is function of A, as b = b(1 -A/A max ), where A max is the carrying capacity of the environment, also the saturation value of A. It means that A cannot grow without superior limit, but it will saturate when approaching A max .

Figure 2 . 1 :

 21 Figure2.1: Scatter plots about the example θ 1 -θ 2 joint distribution in the case of (a) r(θ 1 , θ 2 ) = -0.9; (b) r(θ 1 , θ 2 ) = 0; and (c) r(θ 1 , θ 2 ) = 0.9.

Figure 2 . 2 :

 22 Figure 2.2: (a) Figure 1 in Corrsin (1961): A quantitative sketch of the variance spectrum of reactive scalar undergoing a first order reaction. (b) Figure 1 in Corrsin (1964a): A quantitative sketch of the variance spectrum of reactive scalar undergoing a second order reaction. Here C is the coefficient of reaction rate; k or m is the wave number; G is the wave number dependant scalar variance.

  that the covariance is almost invariant for very slow and very rapid secondorder reactions. In a direct numerical simulation about the influence of convective turbulence on chemical reactions in the atmospheric boundary layer by Molemaker and de Arellano (1998), two species A and B representing the reactants in second order irreversible reaction were introduced. The one-dimensional spectrum of covariance between the reactants (averaged in space and time and normalized) were plotted (figure 2.3), and a -2/3 scaling law was found. Chemical reaction introduces various complexities to the passive scalar

Figure 2 . 3 :

 23 Figure 2.3: Figure 2(b) in Molemaker and de Arellano (1998): one-dimensional spectrum of the covariance between two chemical species (averaged in space and time and normalized).

  carried out a two-dimensional numerical simulation about the temporal evolution of the spatial distribution of reactive scalar in the advection of fluid field. The reactants were in elementary reactions in open or closed flows. The evolution of the reactive scalar field was found to be dependent on the Damköhler number. When the Damköhler number is small, the variance decays and the reactive scalar field is gradually of spatially homogeneous distribution. But when Damköhler number is large, the structure of filaments with strong perturbations soon grows in the reactive scalar field.

Figure 3 . 1 :

 31 Figure 3.1: The comparison of the energy spectra of velocity numerically obtained by pseudo-spectral methods with and without dealiasing. k = |k| is the amplitude of wave number vector. The dotted curve is the filtering function of smooth dealiasing.

3. 1 -

 1 Spatial discretization and approximation of derivatives 47 high orders of accuracy at other positions. For a function F (x) with x ∈ [0, L], the discrete function and coordinate can be F = [F 0 , F 1 , . . . , F N ] and x = [x 0 , x 1 , . . . , x N ] with x 0 = 0 and x N = L. A Dirichlet boundary condition such as F (0) = α can numerically be employed as F 0 = α. For a Neumann boundary condition, which prescribes the gradient of the quantity at the boundary, an asymmetric scheme can be used. If it is a second order center scheme used in the inner domain, the numerical implementation of the Neumann boundary condition ∂F ∂x 0 or L = β at the grid of x 1 and x N can be

Figure 3 . 2 :

 32 Figure 3.2: The illustration of the physical domain of dependence of Eq. (3.31) (thick dashed line) and the numerical domain of dependence (gray area).

Figure 3 . 4 :

 34 Figure 3.4: An illustration for the balance between reaction and diffusion in the reaction zone.

4 . 3 - 4 . 3

 4343 Re λ varies in the range Re λ ∈[20, 150], the Schmidt number spans the interval [0.[START_REF] Wu | Active transport in the Ocean: Turbulence, Chemistry & Biology[END_REF] 4] and the Damköhler number changes from O(10 -4 ) to O(10), while the reaction order n is increased from one up to n = 3 (i.e. from second to fourth order reaction). The values of the key parameters for the simulations are reported in table 4.1. Temporal evolution the mean and fluctuation component of scalar fields 63 Temporal evolution the mean and fluctuation component of scalar fields We begin by looking at the temporal evolution of the two first statistical moments of reactive fields, i.e., their mean values and root-mean-square fluctuations.

Figure 4 .

 4 Figure 4.1 shows the typical temporal evolution of the fluctuating and the mean (volumetric average) parts of scalars for a simulation with n = 2, Da = 0.1, Sc = 1 and Re λ = 150. After a sufficiently long simulation time, a statistically steady state is established where the global mean value for the reactive scalar fields is close to the respective equilibrium quantities, i.e., R 1 ≈ R 2 ≈ P ≈ 1. Furthermore, in spite of the presence of a vigorous external mechanical forcing and random scalar source terms, the reactive scalar dynamics is characterised by relatively small global fluctuations from the equilibrium state. We observe that the scalar fluctuations are proportional to the amplitude of the mechanical forcing, which poses a limitation for the numerical implementation of the model system, i.e. the positiveness of the scalar concentration fields (θ ≥ 0). In order to fulfill this constraint, in the simulations for this chapter the r.m.s. of scalars reaches at maximum 10% of the mean value. The statistical convergence is reached by means of simulations extending in time ∼ 45T I , where T I = k/ε with k = 3u 2 /2 is the integral time scale. The temporal averages are performed after at least 8T I from the beginning of the simulation (see figure 4.1).

Figure 4 . 1 :

 41 Figure 4.1: Evolution of the root mean square of scalar fluctuations and mean values for the case of n = 2, Da = 0.1, Sc = 1 and Re λ = 150. Time is normalized by the integral time k/ε with k = 3u 2 /2. The mean quantities were represented with three-dimensional volumetric average here. The dash vertical line marks the initial time for the computation of statistical quantities.

Figure 4 . 2 :

 42 Figure 4.2: Energy spectra of reactive and passive scalar fields, i.e. E θ (k) with θ = R 1 , R 2 , P, T , and velocity field E u (k) (solid light-blue line) in the condition Re λ = 150, Sc = 1, n = 1 and Da = 10, 1, 0.1, 0.01, 0.001 (from high to bottom). Each spectra is compensated with the KOC scaling, (kη) 5/3 and normalized by the global energy θ 2 . For clarity, the energy spectra of scalars are shifted vertically by a multiplicative factor 0.1

  R 2 , P and T, (4.11b) where • k denotes the average in time and over all the modes in the shell of thickness ∆k centred at k = |k|, ûi (k) and θ(k) are the Fourier coefficients of the mode of k, û * i (k) and θ * (k) are the corresponding complex conjugates.

Figure 4 .

 4 Figure 4.2 depicts the log-log plots of the three-dimensional energy spectra of scalars of a typical case at Re λ = 150, Sc = 1, n = 1, compensated with k -5/3 , which is the scaling expected in the inertial regime both

Figure 4 . 3 :

 43 Figure 4.3: Coherency spectra of the reactive scalars, under the condition of Re λ = 80, Sc = 1, n = 1. The horizontal axis is normalized in terms of the Kolmogorov scale η. The dash vertical line marks the maximum wave number at which the scalar source terms acts.

Figure 4 .

 4 4 presents the global correlation coefficients of the reactive scalars and their gradients against Da with Re λ = 150, Sc = 1, n = 1. The speculation that the correlation coefficients

Figure 4 . 4 :

 44 Figure 4.4: The global correlation coefficients of the reactive scalars and their gradients along x direction, under the condition of Re λ = 150, Sc = 1, n = 1.

Figure 4 . 5 :

 45 Figure 4.5: Correlation coefficients between R 1 and R 2 (r(R 1 , R 2 )), R 1 and P (r(R 1 , P )), R 2 and P (r(R 2 , P )) as functions of Da θ , under the condition of (a): n = 1 (b): n = 2 and (c): n = 3 and the Schmidt number Sc = 1. Theoretical predictions are shown in black lines.

Figure 4 . 6 :

 46 Figure 4.6: Correlation coefficients between R 1 and R 2 (r(R 1 , R 2 )), R 1 and P (r(R 1 , P )), R 2 and P (r(R 2 , P )) as functions of Da θ , under the condition of (a): Re λ = 20 (b): Re λ = 40 and (c) Re λ = 80. The order of R 2 (n) is 1. Theoretical predictions are shown in black lines.

Figure 4 .

 4 Figure 4.6 further confirms the range of validity of the prediction, by displaying the same correlations now for the case of different Schmidt numbers in the range from 0.1 to 4 and for n = 1. Again the trends are well captured by the theoretical predictions.

Figure 4 . 7 :

 47 Figure 4.7: PDF of the scalar quantities (main panel) and normalized with respect to their standard deviations (inset, the black lines are Gaussian curves) under the condition of Re λ = 150, Sc = 1, n = 2 and Da = 0.1.

Figure 4 . 8 :

 48 Figure 4.8: Visualisation of two-dimensional instantaneous sections of the R 1 field at Re λ = 150, Sc = 1 and for three different Da values.

Figure 4 . 9 :

 49 Figure 4.9: PDF of reaction rate Da(R 1 R n 2 -P ) at Da = 0.1 and Sc = 1 for n = 1 (left) and n = 2 (right) for different Reynolds numbers Re λ = 20, . . . , 150. The top panels shows the raw PDF while the bottom ones show the PDF of the normalised variable with respect to its standard deviation. The two insets display the trend of the variance of the reaction rate as a function of the Reynolds number in log-log scale, the red line has a slope 3.

Figure 4 . 10 :

 410 Figure 4.10: The fluctuations of the reactive scalars normalized by the fluctuation of passive scalar (T ) as functions of Da θ , under the condition of (a): n = 1 (b): n = 2 and (c): n = 3 and the Schmidt number Sc = 1. Theoretical predictions are shown in black lines.

Figure 4 .

 4 Figure 4.11: a) Taylor micro scales of scalars (computed on T ) as functions of Sc. The red line draw is 0.3/St. A power law fit of the form aSc b on the Re λ = 150 data set gives a = 0.29 ± 0.1 and b = -0.93 ± 0.02. b) Taylor micro scales of the reactive scalars with respect to the passive scalar one with different order of reaction, as functions of Da for all the simulations at Re λ = 150 and Sc = 1.

Figure

  Figure 4.12(a) shows that the correlation coefficients vary quite weakly with the filter parameter K. The

Figure 4 .

 4 Figure 4.12: (a) Correlation coefficients of scalar evolving in a coarse-grained turbulent flow filed, compared with theoretical predictions and (b) Taylor micro scales of the scalars convected by filtered velocity ( λθ ) as functions of the maximum wave number of filtered velocity (K), at Sc = 1, Re λ = 80 and Da = 0.05 (Da θ = 1.42).

Re 0 FRHIILFLHQWV 0 Figure 5 . 1 :Figure 5 . 2 :

 005152 Figure 5.1: The coefficients obtained by fitting the profiles of mean velocity and kinetic energy, κ in Eq. 5.6 and α and β in Eq. 5.9, as function of Re λ0 .

Figure 5 . 3 :

 53 Figure 5.3: The different normal stresses: one finds u 2z > w 2 z > v 2 z . The z-dependence is given by the fits of equation(5.8).

  Re λ = 122.85

Figure 5 . 4 :

 54 Figure 5.4: The mean kinetic energy K(z) = 1 2 u i u i z .

Figure 5 . 5 :

 55 Figure 5.5: Simulation results for the test of the validity of Boussinesq's hypothesis, representing the alignment ρ RS between R and S. The mean velocity profile is superposed in dotted line for reference.

Figure 5 . 6 :

 56 Figure 5.6: The amplitudes of the terms at the right hand side of Eq.(5.24) as function of z. The mean velocity profile is also represented as a dotted line, for reference. The horizontal red dotted lines mark the 0 value for the amplitudes. The simulation results of Run 1, 3, 5, 7 are shown here.

Figure 5 .

 5 Figure 5.7 shows the numerical estimations of the square of Taylor microscale of the reactive scalars under the condition of Da=8.5, in comparison with the profile of λ 2 . Although it can be found that λ 2 θ is of the same

Figure 5 .

 5 Figure 5.9: (a) Variances of reactive scalars, under the condition of Da=8.5, and passive scalar as function of z. (b) Variances of reactive scalars normalized by the variance of passive scalar as function of Da θ at the position of z = H/2. (b) Variances of reactive scalars normalized by the variance of passive scalar as function of z under the condition of Da=8.5.

Figure 6 . 1 :

 61 Figure6.1: Schematic diagram of the flow configuration and computational domain. For the scalars quantities, the periodic boundary conditions are set along x and y directions, while a Dirichlet boundary condition is used along the z direction. The shadowed layers near the boundaries are the "buffer layers" generated artificially, in which the quantities of scalars are close to the preset boundary values, as defined in Eq. (6.5). The part between buffer layers is denote as the bulk region. Such setup is statistically stationary and ensures the local positiveness of scalar concentrations.

  Figure 6.2:The reaction rates computed with the mean quantities as functions of z. The solid lines are for the forward reaction rates Da 1 R 1 R 2 z and the dashed lines are for the backward reaction rates Da 2 P z . The clear difference can be observed. Vertical dotted lines mark the interfaces between the buffer layers and the bulk region.

Figure 6 . 3 :

 63 Figure 6.3: The three-dimensional instantaneous snapshot of (a) R 1 , and (b) reaction rate (Da 1 R 1 R 2 -Da 2 P ) on the isosurface of R 1 = 0.5, under the condition of Γ = 10.

Figure 6 . 4 :

 64 Figure 6.4: Evolution of the spatial average of the mean and root mean square of scalar fluctuation in the bulk region for the case of Γ = 10. Time is normalized by the integral time k/ with k = 3u 2 /2. The dashed vertical lines marks the initial time for the computation of statistical quantities.

Figure 6 . 5 :

 65 Figure 6.5: Numerical profiles of: (a) the mean, and (b) the r.m.s. of T with different buffer layers thickness δ. The vertical dotted lines mark the interfaces between the buffer layers and the bulk region.

D T , 2 DFigure 6

 26 Figure6.6: (a) Turbulent diffusivity calculated from Eq. (6.10). (b) Theoretical prediction of the mean of T (dashed lines) compared with the DNS results (solid lines with same color). The vertical dotted lines mark the interfaces between the buffer layers and the bulk region.

Figure 6 . 7 :

 67 Figure 6.7: Dependence of PDFs at z = 1/2 on Γ for: (a) R 1 and T , and (b) P . The insert panel in (b) plots the peak of p P as function of Γ.

Figure 6 . 10 :

 610 Figure 6.10: The instantaneous two-dimensional snapshot of reaction rate at the position of z = 1/2. The upper row (a,b,c) correspond to Γ = 100; the lower row (d,e,f) correspond to Γ = 1. The first column (a,d) show the net reaction rate (Da 1 R 1 R 2 -Da 2 P ); the second column (b,e) show the forward reaction rate (Da 1 R 1 R 2 ); the third column (c,f) show the backward reaction rate (Da 2 P ).

  Figure6.12: The scalar mean: (a) R 1 z and (b) P z , as a function of z obtained from theoretical analysis (dashed lines) based on Eq. (6.25) and DNS (solid lines with the same colors). The grey dashed lines are from the theoretical prediction at infinitely large Da 1 according to Eq. (6.20). We see that the prediction for Γ = ∞ is close to the curves for Γ = 100 and also the predictions for large Γ are close to the DNS results when Γ=10, 30, 100. The vertical dotted lines mark the interfaces between the buffer layers and the bulk region.

Figure 6 .

 6 Figure6.14: DNS simulations of the correlation coefficients between (a) R 1 and R 2 ; (b) R 1 and P as function of z, for different Γ cases. The vertical dotted lines mark the interfaces between the buffer layers and the bulk region. correlated (with less negative r(R 1 , R 2 )) at larger Γ, but more correlated (with more negative r(R 1 , R 2 )) at smaller Γ. The reason can be given as follows. For the non-reactive case with Da 1,2 = 0, the influence from the chemical reaction on the correlation coefficient vanishes. In the present flow configuration, the fluctuations of R 1 and R 2 are forced exclusively by their respective mean gradients, which are opposite with each other. Thus R 1 and R 2 are perfectly negatively correlated. In the bulk region with a stronger chemical source Da 1 R 1 R 2 , as Γ increases, the chemical reaction depletes more the reactants, resulting in the skewness of the PDFs of of R 1 and R 2 toward the R 1,2 = 0 side, as shown in figure 6.7. Therefore, the product R 1 R 2 tends to be more positive. As a result, according to the definition in Eq. (6.26), r(R 1 , R 2 ) increases, i.e. R 1 and R 2 become less negatively correlated. Interestingly, in the buffer layer when Γ is sufficiently large (e.g. Γ = 100) r(R 1 , R 2 )

Figure 6 . 15 :Figure 6

 6156 Figure6.15: The ratio of the energy spectra of (a)R 1 at z = 1 4 ; (b) R 1 at z = 1 2 ; (c) R 1 at z = 3 4 ; (d) P at z = 1 4; (e) P at z = 1 2 ; (f) P at z = 3 4 to the energy spectrum of T at the same z.

  Such a state is maintained by the competition between the reversibility of reaction, which enforces a global chemical equilibrium in the statistically stationary state, and the variability of the scalar forcing, which introduces stochasticity into the scalar fields. The overall reaction Conclusion and future works rate is small in such a dynamical equilibrium state, and thus too weak to change the scaling behavior and statistical distribution of the scalars. Moreover, the fluctuations of the scalars are found to be small compared with the mean quantities. Thus a linearization of the reaction terms in the scalar equations is a viable analysis approach. Based on the linearization of the reaction terms and other appropriately proposed hypotheses, modeling analysis about the correlations and variances of fluctuations of reactive scalars were carried out. A theoretical analysis showed that the correlation coefficients and the fluctuations of reactive scalars are functions of a unique control parameter: the Damköhler number (Da θ ), which is constructed as the ratio between the time scale of scalar diffusion across a domain of the size of the scalar Taylor micro-scale (λ 2 θ /D) and the chemical reaction time scale τ r . The larger is such a Damköhler number, the more depleted are the scalar fluctuations as compared to the fluctuation of a passive scalar field in the same conditions, and the more intense are the correlations. A saturation in this behaviour is observed beyond Da θ O(10). It indicates that the chemical processes tend to reduce reactant concentration fluctuations and enhance their correlation intensity. While on the contrary, the turbulent mixing increases fluctuations and removes relative correlations. Importantly, the functional dependence of the fluctuations and the correlations of the scalar quantities characterized by Da θ are well confirmed by the DNS results in the full range of explored conditions with variable reaction order, the Reynolds number and the Schmidt number. The key role of Da θ also implies the significance of the scalar Taylor micro-scale λ θ in the mixing of chemical species. A deeper insight into λ θ reveals that Da θ can also be viewed as the ratio of the large-eddy-turnover time of the flow over the chemical reaction time.

  Moments of the reactive scalars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3.4 Correlation Coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3.5 Scalar Energy Spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Summary of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Scientific contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Future works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	CONTENTS	v
	6.3.3 7 Conclusion and future works
	7.1 List of Figures

1.2.1 Advection-diffusion equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Obukhov-Corrsin theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.3 Probability density function, structure function and intermittency . . . . . . . . . . . . . 2 Reactive scalars in incompressible turbulence 2.1 Chemical kinetics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 General descriptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.2 Models in chemical and biological systems . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Advection-diffusion-reaction equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Interactions between turbulence and reactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Analytical studies on reactive scalars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Numerical studies on reactive scalars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.3 Turbulence and reactive scalars in biological oceanography . . . . . . . . . . . . . . . . . 3 Numerical implementation of direct numerical simulations 3.1 Spatial discretization and approximation of derivatives . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Pseudo-spectral methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Introduction of the model system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Numerical implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Result analyses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3.1 Properties of the buffer layer and the bulk region . . . . . . . . . . . . . . . . . . . . . . . 6.3.2 Statistical properties of the reactive scalars . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Drawing of turbulence by Leonardo da Vinci (15th century). . . . . . . . . . . . . . . . . . . . . 1.2 Illustration of turbulent cascade (Ardeshiri, 2016). L is the length scale of domain and η is the smallest length scale, kolmogorov microscale. Energy is injected from outside of the flow at the largest scales. Then the largest eddies break up into eddies of smaller and smaller size and energy is meanwhile transferred. This process continues until the smallest eddies and the energy is dissipated by viscosity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 The probability density function of normalized velocity (upper curve) and temperature (lower curve) in a grid-generated turbulence (shown in the work of

  Thus δ can be estimated as D γ R

		-n/2 2	. The factor of R	-n/2 2	is for the special case, but
	generally, δ roughly satisfies
	δ ∝	D γ	.		(3.42)

.41) Since the reaction zone is supposed to contain the local maximums of the reactant concentrations, which are close to zero outside, the R 1 and R 2 in Eq. (3.41) can be considered as local maximums and the ∂R1 ∂x can be Numerical implementation of direct numerical simulations estimated as -R 1 /δ.

Table 4 . 1 :

 41 Parameters for the simulations: Re λ is the Taylor scale based Reynolds number; Sc is the Schmidt number; Da = λγ 1 R n 2,eq /u = λγ 2 /u is the Damköhler number based on Taylor scale; n is the order of R 2 in the reaction; N 3 is the total number of grid points; η is the Kolmogorov length or dissipative length; |k| max is the maximum wave number amplitude kept by the dealiasing procedure; |k| max • η is the spatial resolution condition; dt/τ η is the time step normalized by the Kolmogorov time scale τ η .

	No.	1	2	3	4
	Re λ	20	40	80	150
	Sc	0.1-4	0.1-4	0.1-1/2-4	1
	Da	0.0005 -50	0.0003 -30	0.0005 -20	0.001 -10
	n	1	1 -3	1 -3	1 -3
	N 3	64 3	64 3 * /128 3 †	128 3 * /256 3 †	256 3
	|k| max η	3.12	1.48 * /2.95 †	1.26 * /2.52 †	1.05
	dt/τ η	0.044	0.06 * -0.03 † 0.044 * -0.022 †	0.034
	2011;

1

  Re λ Sc |∇R 1 | 2 , can be estimated similarly as Eq. (4.26). Under the isotropic 4.6 -Reactant variances 77 and statistical stationary conditions, the turbulent energy of R 1 , R 2 and P are approximately determined as

Table 5 . 2 :

 52 The numerical values of the coefficients in Eq. 5.8 for each run.

	Consequently,

  The only non-zero shear stress term is proportional to the cosine function, and the normal stress components all involve a square cosine expression. The normal stresses are never equal, showing that as expected the turbulence is

	DPSOLWXGHVRIWKHWHUPV	2 a ||S|| 2 u 2 w a 2 ||T2|| 6 2 v 4K a 2 ||T3||	u z	DPSOLWXGHVRIWKHWHUPV	2 a ||S|| 2 u 2 w 6 2 a 2 ||T3|| v 4K a 2 ||T2||	u z
	z			z
	(a) Re λ = 38.69		(b) Re λ = 49.31
	DPSOLWXGHVRIWKHWHUPV	2 a ||S|| 2 u 2 w a 2 ||T2|| 6 2 v 4K a 2 ||T3||	u z	DPSOLWXGHVRIWKHWHUPV	2 a ||S|| 2 u 2 w 6 2 a 2 ||T3|| v 4K a 2 ||T2||	u z
	z			z
	(c) Re λ = 66.9		(d) Re λ = 122.85

Table 5 . 3 :

 53 Key global parameters after reaching statistically stationary state. Re λ0 = u λ 0 /ν is the global Taylor microscale, where u is the global root mean square of single component velocity, ν is kinematic viscosity, λ 0 is the global Taylor micro scale; Sc is the Schmidt number; n is the reaction order of R 2 in the reaction; N 3 is total number of grids; η 0 = (ν 3 / 0 ) 1/4 is the global Kolmogorov length scale; |k| max is the maximum wave number amplitude kept by the dealiasing procedure, |k| max • η is the resolution condition; Da = λ 0 γ 1 R n 2,eq /u = λ 0 γ 2 /u is the Damköhler number based on Taylor scale.

	Re λ0 Sc n λ 0 u	N 3	η 0	|k| max • η 0	Da
	54.7	1	1	1	1 128 3 0.052	1.55	0.0085-8.5

Table 6 .

 6 Sc is the Schmidt number (ν/D); N 3 is the number of total grid points; |k| max • η is the resolution condition, where |k| max is the maximum amplitude of wave number kept by the dealiasing procedure, η if the Kolmogorov lenght scale; τ η is the Kolmogorov time scale; Γ = Da 1 /Da 2 , with Da 1 and Da 2 as the Damkholer numbers for forward and backward reactions respectively; L I is the integral length scale; T I is the integral time scale; ∆t is the numerical time step.

1: Non-dimensionalized parameters for the simulations: Re = u H/ν is the Reynolds number based on large scale, where u is the single-component root-mean-square velocity, H is the length of the domain, ν is the viscosity; Re λ = u λ/ν is the Taylor scale λ based Reynolds number;

Reactive scalars in incompressible turbulence

Acknowledgments

This thesis is under the joint support by Shanghai Jiao Tong University and French Region "Hauts-de-France" in the framework of a cotutella PhD programme. I will be eternally grateful of such a wonderful opportunity for improving my academic ability and broadening my horizons.

It can never be enough for me to express my gratitude to my advisors: Mr. Lipo Wang, Mr. François G. Schmitt and Mr. Enrico Calzavarini. Thank you for every comment you have given on my works. Thank you for every knowledge you have taught me. Thank you for every major or minor correction you have made to my writing. Thank you for every help you have given to my study and life. Thank you for every criticism on me. In my mind, you are most respectful because of your profound knowledge and meticulous attitude in science. It is absolutely my great honor to study under your supervision. This experience is the unforgettable memory and priceless fortune in my whole life. I would like to thank Dr. Michael Gauding (CORIA (CNRS UMR 6614), Rouen, France) for providing the original code and patience in the further instruction; thank Dr. Yongxiang Huang (Xiamen University) for the kind academic suggestions or friendly help in life; thank all the staffs in UM-JI, especially Fiona Li, for your sweet help in my daily affairs. I acknowledge the computing resources for the numerical simulations in this thesis including the cluster in UM-SJTU Joint Institute, the CALCULCO computing platform supported by SCoSI/ULCO (Service COmmun du Système d' Information de l' Université du Littoral Côte d' Opale), the National Supercomputer Center in Guangzhou, the High Performance Computing Center (HPCC) at Université de Lille and the Center for High Performance Computing at Shanghai Jiao Tong University.

I would like to also thank all my friends who made me feel not lonely during my PhD: Himani Garg, Ali Amar Katia, Hoan, Xiaodan Cao, Yang Gao, Peipei Zhao, Shibo Gu, Qianhui Li, Xiang Li, Hui Cao, Yuqing Guo, Kun Wang, Yu Cang, Guiwen Tan, Qun Hu, Chenlin Guo, Kunlin Li, Chun Gong . . . Thank goodness for letting me acquaint all of you.

Chapter 6. This configuration has the advantage of producing strongly out of equilibrium chemical reactions.

Finally, in Chapter 7, a summary of the thesis and some ideas about the future works are identified.

Test of a quadratic constitutive equation

We have seen above that the linear model cannot produce an anisotropic Reynolds stress tensor for anisotropic flow such as the TKF. [START_REF] Pope | A more general effective-viscosity hypothesis[END_REF] has proposed to use invariant theory in turbulence modeling, representing the stress tensor as a development into a tensor basis. Originally it was on the form R = 10 i=1 a i T i with 10 basis tensors. By considering a quadratic development, only three tensors are used. This is complete for two dimensional flows [START_REF] Pope | A more general effective-viscosity hypothesis[END_REF], and is also a good approximation for fully 3-dimensional flows [START_REF] Jongen | General explicit algebraic stress relations and best approximation for three-dimensional flows[END_REF]. It can be used for channel flows [START_REF] Schmitt | Direct test of a nonlinear constitutive equation for simple turbulent shear flows using dns data[END_REF][START_REF] Modesti | A priori tests of eddy viscosity models in square duct flow[END_REF] and we propose to use it also here for the TKF.

In this framework, the anisotropic stress tensor writes as a three-terms development

where the three tensors of the basis are all symmetric and traceless [START_REF] Pope | A more general effective-viscosity hypothesis[END_REF]:

(5.17)

The coefficients a i can be written using scalar invariants of the flow, which correspond to scalar fields whose values are independent of the system of reference. Invariants can be defined as the traces of different tensor products [START_REF] Spencer | Part iii. theory of invariants[END_REF]. Some of the first invariants are the following:

All these invariants can be here estimated numerically. Furthermore, the coefficients a 1 , a 2 and a 3 can be expressed using the invariants. This is done by projecting the constitutive equation onto the tensor basis: successive inner product of this equation with tensors T i provides a system of scalar equations involving the invariants [START_REF] Jongen | General explicit algebraic stress relations and best approximation for three-dimensional flows[END_REF]. For 2D flows such as TKF, we have η 3 = 0 and η 5 = η 1 η 2 /2, and the system of scalar equations is inverted to provide finally the quadratic constitutive equation using invariants:

(5.18)

For the TKF the invariants write: (5.19) (5.20) shape with λ 2 , their fluctuations are hugely different. The value of λ is proportional to the mean gradient of velocity ∂ u z ∂z , and therefore of great variety. However, λ 2 θ is not of large fluctuation, within 10% of the mean value. Thus Da θ can be reasonably considered as independent on z and scalar. According to Eq. ( 4.15), the Da θ can be roughly considered as independent of z.

Correlation coefficient

The correlation coefficients between reactive scalars are defined based on the fluctuating parts of scalars as

Here θ 1 and θ 2 are the scalars under consideration. The statistics of the correlation coefficients are implemented for each specific z position. Thus they are functions of z.

According to the theoretical prediction about the correlation coefficients between scalars based on Da θ given in Chapter 4, under the condition of n = 1, we have

(5.27)

In figure 5.8(a), the correlation coefficients between reactive scalars at an example position are shown, and found to be in excellent agreement with the theoretical predictions by Eq. (5.27). Moreover, as shown in figure 5.8(b), the correlation coefficients between scalars are independent of z. This result further confirms that the theoretical model about the correlations between reactive scalars proposed in Chapter 4 relies weakly on the specific flow environment.

Variances of scalars

Similar with the correlation coefficients, the scalar fluctuations are also functions of z in the TKF. Figure 5.9(a)

illustrates such a dependence. The fluctuations of reactive scalars are evidently smaller than that of the passive scalar. It indicates that, like in the homogenous isotropic turbulence, the reaction here also tends to dump the fluctuations of scalars. In addition, the profile of the scalar fluctuations are toughly of sinusoidal shapes, although the amplitude is small compared with the mean. A theoretical prediction about the fluctuations of reactive scalars normalized by that of the passive scalar was proposed in the modeling analysis in Chapter 4. These quantities in TKF are shown in figure 5.9(b) as functions of z. Roughly, they can be considered as independent of z. Following the procedure described in section 4.6.2, under the condition of n = 1, the analytical predictions about the dependence of the normalized fluctuations of reactive scalars, at any z in the by the correlation between R 1 and R 2 , larger Γ will reduce in most of the flow field their coexistence, which explains the stripe like distribution in figure 6.10 (a). Since for the present chemical kinetics the backward reaction rate is solely determined by P , the effect of Γ on the scalar correlation is not relevant in determining the backward reaction rate. Therefore figure 6.10 (c) and (f) are almost identical. In summary, under different Γ the PDF and spatial distribution of the net reaction rate will be mainly determined by the forward part.

Moments of the reactive scalars

The reactive scalar θ (e.g. R 1 , R 2 or P ) can be decomposed into the z-dependent mean part and the fluctuating part as θ(x, t) = θ z (z, t) + θ (x, t), whose numerical results are shown in figure 6.11.

Different from the linear profile of the passive scalar T , the profiles of R 1 z are concave in the bulk region, because of chemical consumption of R 1 with R 2 . With increasing Γ, R 1 z decreases while P z increases, because stronger forward reaction depletes more R 1 and produce more P . The scalar means tend to saturate at the infinite large Γ. The main panel of (c) shows the mean profile of P normalized by its maximum, whose function as Γ is plotted in the inset plot. There is a prefect superposition for all Γ values. In all the plots, the vertical dotted lines mark the interfaces between the buffer layers and the bulk region. as function of z obtained from theoretical analysis (dashed lines) and DNS (solid lines with the same colors). The grey dashed lines are from the theoretical prediction at infinitely large Da 1 according to Eq. (6.21). We see that the prediction for Γ = ∞ is close to the curves for Γ = 100 and also the predictions for large Γ are close to the DNS results when Γ=10, 30, 100. The vertical dotted lines mark the interfaces between the buffer layers and the bulk region.

Last, I would like to give my most sincere gratitude to my family. Thank you for raising me up. Thank you for your unconditional support. Thank you for being an inseparable part of my life. Thank you for your love.