Langevin approach to generate synthetic turbulence

We present an analytical scheme, easily implemented numerically, to generate synthetic Gaussian turbulent flows by using a linear Langevin equation, where the noise term acts as a stochastic stirring force. The characteristic parameters of the velocity field are well introduced, in particular the kinematic viscosity and the spectrum of energy. As an application, the diffusion of a passive scalar is studied for two different energy spectra. Numerical results are compared favorably with analytical calculations.


I. INTRODUCTION
The statistical approach to turbulence has a long history on its own [1,2]. Actually, it is the firm recognition that both fundamental and applied aspects of turbulence can be consistently addressed from the knowledge of the statistical moments and correlations of the velocity flow that has made this approach particularly appealing. Since in any case such an statistical point of view has adopted multiple perspectives during these last decades, let us first put our formulation into an appropriate context.
Our primary aim here is to describe an analytical and numerical methodology to generate a statistically stationary, homogeneous, and isotropic two-dimensional turbulent flow with zero mean velocity and well-defined energy spectrum. Needless to say, such a regime of steady turbulence can only be maintained by means of an external input of energy to compensate the dissipative nature of the viscous forces. This external input could be conceptually associated with stirring forces which are stochastic in nature so as to produce a random velocity field which is to represent the turbulent flow [3]. It is quite obvious that the statistical properties of the steady stirring forces will ultimately determine those of the turbulent flow.
Ideally, the statistical properties of any turbulent flow should come as the output of a first principles, Navier-Stokes based, formulation of the problem. However, we will adopt here a somewhat reversed perspective aimed at developing a methodology to construct what would be a sort of synthetic turbulence [4][5][6][7][8]. Rather than to retain the nonlinear coupling which makes possible the redistribution of energy from the largest length scales down into the smaller ones, we assume that the energy is incorporated into the system in an individual wave number basis. In other words, our model would represent a collection of uncoupled stirrers each one acting on a single-length scale and introducing its own wave number dependent energy contribution, but chosen to reproduce, and this is the main point in our approach, the desired spectral distribution in steady state.
This practical approach to produce a turbulent flow is justified by the fact that what we want is to study a physical process within the turbulent medium described only by its statistical properties. Our strategy is accomplished by using a generalized Ornstein Uhlenbeck-like Langevin equation for the stream function η(r, t) [4,6]. Being more specific, our proposal is based on the use of the following Langevin equation where ν is the kinematic viscosity and ζ(r, t) is a Gaussian white noise with zero mean value and correlation In this last expression, ǫ 0 , and λ are parameters which control, respectively, the intensity and correlation length of the random flow. Q[λ 2 ∇ 2 ] will play the most relevant role in our scheme, representing the random stirring forces. The random flow generated in this way has Gaussian properties due to the linear nature of such an equation. Apart from that, other limitations concerning some invariances of the turbulent flow appear also as a consequence of such a linearity [6]. However, the big advantage of our model is that it facilitates the control over the characteristic parameters of the turbulence, i.e. its integral time and length scales and its spectrum, just by appropriately prescribing the input parameters of the noise entering into the Langevin equation. A final remark concerning Gaussianity is worth mentioning at this point. Certainly, we recognize that such an statistical property is under scrutiny [9]. However particular experimental scenarios supporting it [10], together with the lacking of conclusive results on intermittency [3], leave this question rather open and allows us to use such a Gaussian property, at least as a working simplifying hypothesis, mainly when the focus is in the study of physical process inside this medium. This paper is organized as follows. In the next section we summarize the Langevin approach to generate the random flow (synthetic turbulence). In Sec. III the technical details of the numerical simulations are presented. As a sort of application we consider the classical problem of the diffusion of a passive scalar in Sec. IV. We devote Sec. V to draw some conclusions and perspectives.

II. LANGEVIN APPROACH
As our approach starts with the generation of a scalar stream function [11], we are going to review first the relationships between the stream function properties and those of the velocity field. Let η(r, t) be the stream function from which we define the two-dimensional incom- Its mean value is to be taken zero and the two-point velocity autocorrelation is defined as usual through where homogeneity in space and time is explicitly denoted [2,3]. Our scheme is entirely implemented in two dimensions (2D). However, it could be generalized to 3D just by taking a vector stream function η and obtaining the velocity field as v = ∇ × η. In this case each component of the stream function must satisfy a Langevin equation similar to (1) with independent noises. Making use of the spatial isotropy we define the radial correlation function as where r = |r 1 − r 2 | and s = |t 1 − t 2 |. The autocorrelation of η(r, t) is correspondingly assumed to have the properties of homogeneity, isotropy, and stationarity. From the definition (3) we can express the velocity correlation (5) in terms of (6) where J 0 (kr) is the Bessel function of zeroth order and C(k, s) is the Fourier transform of C(r, s).
The physical parameters of steady turbulence, i.e., its intensity u 2 0 , and characteristic (integral) time and length scales follow from their standard definitions [2,3], With these definitions in mind let us move to the discussion of the analytical scheme. The Fourier transform of (1) reads Now we introduce the structure function S(k, t) related with the equal-time correlation function (6), but in Fourier space By using standard stochastic calculus [12], S(k, t) is shown to verify On the other hand,C(k, s) obeys, in the steady state, the equation with the initial condition From (11) and (12) we get The energy spectrum E(k, t) is defined in terms of S(k, t) as Using (11) we can also obtain the equation of evolution where This quantity can be regarded as the input of energy due to the stirring forces. The stationary state is achieved when the input term is exactly balanced by the dissipation term, The stirring intensity u 2 0 can also be related with the spectrum and the characteristic time t 0 (9), and the characteristic length (9) can be related as well with the viscosity ν and λ. Explicit relations will be obtained later on for particular spectra.
As it is quite obvious, our scheme does not incorporate the non-linear term of the Navier-Stokes equation. Notice, nevertheless, that this scheme is versatile enough to reproduce a large variety of energy spectra, just by appropriately prescribing the differential operator Q[λ 2 ∇ 2 ] in (1). In our approach the noisy term represents the stirring mechanism which feeds energy continuously into the system according to the chosen spectrum. At the same time this energy is dissipated at time and spatial scales controlled by ν. As there are only linear terms in the Langevin equation, no kind of energy cascade is possible. The noise forces introduce the whole hierarchy of turbulent structures evolving according to their own time and length scales, but without any interaction between them.

III. SIMULATION ALGORITHM
We have chosen for the discretization of the real space a standard two-dimensional square lattice N × N with elementary unit spacing ∆ in both directions. In most of our simulations N = 128 and ∆ = 0.5, unless other values are specified. The discrete Fourier space is discretized accordingly k = (k x , k y ) = 2π/N ∆(µ, υ) [14] (Greek indices are used in Fourier space on all that follows). We can take advantage of the fact that (9), when written in the discrete Fourier space, transforms into a set of linear and not coupled ordinary differential equations. In these circumstances the exact integration between t and t+∆t, (∆t = 0.1t 0 ) gives in terms of new random variables β µυ (t) and γ µυ (t) are defined according to where Q µυ in the last two expressions denote the discrete Fourier transform of the operator Q[λ 2 ∇ 2 ]. In the Fourier space the derivative operators have been translated into The correlation of the random variables can be expressed as and similar equations for the γ µυ (t). Using (2), the symmetry properties of the Q operator Q µυ = Q −µ−υ = Q * µυ , and further integrating and using the expressions defined above for the discrete operators in Fourier space we finally have We can now construct an explicit expression for β µυ (t) adapted to the result just obtained where α µυ are Gaussian random numbers, which satisfy α * µυ (t)α ρσ (t) = δ µρ δ υσ . We would proceed analogously for γ µυ (t). The correlation of the stream function in the steady state is When dealing with the diffusion of the passive scalars in Sec. IV we will always start our simulations in a steady configuration of the random flow. According to the result (28) this is easily accomplished by taking as initial condition At each time step and after having generated the stream function in Fourier space we proceed to antitransform and using it in relation with the appropriate discretized forms for the velocity field. We skip some details to refer directly to the the discrete version of the energy spectrum which finally reads [4]

IV. DIFFUSION OF PASSIVE SCALARS
Before we start the study of the scalar diffusion, we will discuss the selection of spectra. There are several possible choices and between all of them we have selected two spectra which are well behaved in all the range of wave numbers.

A. Kraichnan's and Kárman-Obukhov spectra
The spectrum introduced by Kraichnan (K) [15] describes turbulent velocity fields with a widely distributed band of excitations and a peak centered at well-defined wave number ∼ k 0 In this case, the operator Q has to be chosen according to where λ = k −1 0 . Using (14) and substituting in (7), the velocity correlation function in steady state reads R(r, s) = ǫ 0 8π(λ 2 + νs) 2 From it, we can obtain u 2 0 , and the integral time and length scales according to The dimensionless Reynolds number, defined according to Re = l 0 u 0 /ν, is thus expressed in terms of the noise parameters as Our second choice is the Kárman-Obukhov's (KO) spectrum which was introduced [16] to study Kolgomorov turbulence with a long " − 5/3 ′′ tail in the spectrum for large k. To this end we select the family of Kárman-Obukhov spectra, whose general behavior is For convenience we have taken n = 3. The choice of the Q operator is then In this case we obtain the following results for the three basic parameters The corresponding Reynolds number is expressed as A closed expression for the correlation function cannot be obtained in this case. These two spectra have been simulated according to the recipes of the previous section. In Fig. 1 we have plotted the evolution of an initially flat spectrum towards its final steady state. By looking at the patterns corresponding to different time intervals we can judge the role of the kinematic viscosity. In particular the modes with smaller k relax slower than modes with larger k. Continuous and discrete analytical results are also compared in the steady state to see the differences introduced by the discretization procedure specially for large k. The clear differences between both spectra are going to influence the diffusion of passive scalar as it will be discussed later on.

B. Scalar diffusion in a continuous scheme
Scalar diffusion in turbulent flows is a classical problem [17] which still deserves a lot of attention [18]. Our starting point is the Eulerian equation of motion for a scalar distribution ψ(r, t) which is assumed to be passively advected by the previously prescribed isotropic homogeneous and stationary random flow v(r, t), In the standard notation used in (44), D stands for the "bare" molecular diffusion coefficient. By averaging over realizations of v(r, t) we simply obtain from the above equation the temporal evolution of < ψ(r, t) >. Actually taking a δ-like initial condition, this quantity is nothing but the probability density for the spatiotemporal dispersion of a unit amount of the randomly advected scalar [19,20]. Thus the first nonzero moment, is all that we need to compute an effective diffusion coefficient. The procedure outlined above, although simply enunciated, is quite involved in its analytical resolution.
The best way to proceed is thus to follow the standard analytical strategies furnished by the theory of Gaussian stochastic processes. The central difficulty arises from the non-Markovian nature of the process at hand. Controlled perturbative schemes are thus necessary. In particular a consistent expansion, [20] based on the smallness of the correlation time, t 0 , leads to a closed equation for < ψ(r, t) >, linear in the autocorrelation tensor, from which a diffusive regime is identified through the common linear law for the scalar dispersion < r 2 > − < r > 2 =< (∆r) 2 >= 4D eff t = 4(D + ∆D)t. (46) Moreover from such an expansion, the explicit expressions for the leading contribution to ∆D can be simply evaluated as where R ′′ (0, s) = ∂ 2 R(r, s)/∂r 2 | r=0 . Computing these left integrals for both spectra we end up with a common expression which reads In particular for the Kraichnan's flow An analogous expression would be obtained for the Kárman-Obukhov's spectrum.
Both expressions identify the zeroth-order contribution, u 2 0 t 0 , which can be also viewed as an exactly correct limiting case of the classical Roberts analysis [21] for the diffusion of a scalar field advected by a rapidly varying random velocity field.

C. Discrete scheme
Needless to say that given the discrete nature of our simulations, the numerical results for the scalar dispersion will be more favorably compared when referring to the discrete version of the analytical results given above. In particular (47) transforms into For the spectra here analyzed the explicit expression reads

D. Results
Numerical simulation proceeds first by constructing the random velocity field with the desired statistical properties and then by seeding the system with a δ-like initial condition for the dispersed scalar at the center of the lattice. Equation (44) is numerically simulated by using a first-order Euler algorithm for the time variable and symmetric forms for the derivative operators. The values of ∆ and ∆t where those of Sec. III which satisfy numeric stability criteria.
As anticipated in Sec. III, the initial condition η(r, 0) is chosen to correspond to the steady state of η(r, t) (understood here in a statistical sense) (see Eq. (29)). In this way we are sure to be in an isotropic and homogeneous turbulent environment from the beginning of the simulation. Randomly advected by the turbulent velocity field the scalar spreads over the lattice, Fig. 2. At each time step we measure the variance < (∆r) 2 > and we fit, after transients, its temporal evolution to a linear law, to obtain D eff . Summarized in Fig. 3 we present two series of results for ∆D corresponding to two different choices for the correlation time of the random flow t 0 . Simulation results are compared with the appropriate, perturbative obtained theoretical predictions both in their continuous and discrete versions. As evidenced in that figure the agreement is entirely satisfactory. The general conclusion is that a large enough relative correction to molecular diffusion, up to 25% in Fig. 3, can be accurately predicted as long as the correlation time of the random flow, t 0 , is small enough relative to the typical time scale for the scalar advection l 0 /u 0 . In addition we note that the differences between the discrete and continuous theoretical results are easily interpreted when comparing the corresponding expressions for ∆D. Actually the correlation R(r, s) evaluated at the origin r = 0 or to first neighbors may be significantly different given the small value of λ here employed (λ = 1) relative to the elementary spacing ∆ (∆ = 0.5) chosen in all our simulations.
The role of the parameter t 0 is better analyzed in relation with the results of Fig. 4. Plotted against the inverse of the value of the viscosity, a direct measure of t 0 for l 0 fixed as prescribed here, the effective scalar dispersion increases with ν −1 . As expected from the intrinsic limitation of the perturbative approach here developed, numerical simulations and theoretical predictions, although showing similar trends, differ progressively as t 0 increases. Complementary to the results depicted in Fig. 4, chosen to consider the role of correlation time of the random flow, we propose Fig. 5 to examine the influence of the correlation length l 0 . Our theoretical result Eq. (4.19) predicts a very small dependence on l 0 , which is the behavior observed in Fig. 5. Note in this respect that in going from l 0 = 1.0 to l 0 = 3.0, ǫ 0 had to be increased nearly two orders of magnitude to keep u 2 0 and t 0 fixed, but ∆D hardly changes.
The last worthy remark refers to the comparison of K and KO spectra. To this end we proposed in Fig. 6 a standard representation of ∆D vs. u 2 0 t 0 for both flows with identical integral values l 0 and t 0 . What we see is that the KO spectrum has more statistical dispersion and lower values of the effective diffusion than K spectrum. These facts can be understood by comparing both spectra in Fig. 1. The most remarkable difference between both spectra concerns the inertial subrange they mean to represent: KO spectrum is much broader than K wavenumber dispersion. This in turn is directly reflected in the different behavior of both spectra for large k. Thus KO spectrum shows a richer variety of turbulent structures at short distances. For large k, KO spectrum still allows structures of the same order than the lattice size ∆. This is not the case of K spectrum where very small structures have a very low weight as can be seen in the value of the intensity of both spectra for k = 2. So in KO spectrum there is less energy concentrated in the interval of those k ′ s where the maximum of the energy takes place. This would finally result in a reduction of the effectiveness of stirring in the case of KO spectrum, since small k-eddies can be thought to be more effective in dispersing the scalar.

V. SUMMARY AND PERSPECTIVES
An stochastic method, based in a Langevin equation, to generate Gaussian synthetic turbulent flows has been presented. Characteristic parameters of the flow such as its intensity and integral time and space scales are well controlled. A relevant aspect of our method is that multiple choices of flow spectra can be generated. Two of those spectra are explicitly presented here. As a practical application of this methodology we have considered the study of the diffusion of a passive scalar under the influence of two different flow spectra. The role of the different parameters, in particular the kinematic viscosity, which control the characteristic of the flow is discussed.
The generation of well-controlled flow spectra can be a very useful tool in other problems of practical application. For example, this approach was already used in the study of phase separation dynamics under stirring [23]. Reactive fronts, i.e., flames propagating under turbulent convection is presently under study following also this technique [24].