Dynamical Simulation and Statistical Analysis of Velocity Fluctuations of a Turbulent Flow behind a Cube

A statistical approach for the treatment of turbulence data generated by computer simulations is presented. A model for compressible flows at large Reynolds numbers and low Mach numbers is used for simulating a backward-facing step airflow. A scaling analysis has justified the commonly used assumption that the internal energy transport due to turbulent velocity fluctuations and the work done by the pressure field are the only relevant mechanisms needed to model subgrid-scale flows. From the numerical simulations, the temporal series of velocities are collected for ten different positions in the flow domain, and are statistically treated. The statistical approach is based on probability averages of the flow quantities evaluated over several realizations of the simulated flow. We look at how long of a time average is necessary to obtain well-converged statistical results. For this end, we evaluate the mean-square difference between the time average and an ensemble average as the measure of convergence. This is an interesting question since the validity of the ergodic hypothesis is implicitly assumed in every turbulent flow simulation and its analysis. The ergodicity deviations from the numerical simulations are compared with theoretical predictions given by scaling arguments. A very good agreement is observed. Results for velocity fluctuations, normalized autocorrelation functions, power spectra, probability density distributions, as well as skewness and flatness coefficients are also presented.


Introduction
In spite of considerable progress in computer technology, numerical methods, and turbulence modeling during the last several decades, reliable prediction of complex turbulent flows at high Reynolds number remains an elusive target. Turbulent flows and, in particular, wall-bounded flows exhibit wide ranges of spatial and temporal scales that increase with Reynolds number. Tsherefore, direct numerical simulation (DNS) is limited to relatively low Reynolds numbers, as the number of grid points required for DNS increases proportionally to Re 9/4 [1]. Smaller timesteps result in an extra 3/4 power for the total cost scaling as Re 3 . Due to the limitations of DNS, great expectations have been placed on large eddy simulation (LES) [2]. Large eddy simulation is an important technique in the study of turbulent flows. In LES, the governing equations are spatially averaged allowing the large-scale motion to be solved. On the other hand, from this averaging process, the so-called subgrid terms arise, requiring constitutive models to their evaluation [3]. LES requires less computational effort than direct numerical simulations, which attempts to solve all scales present in the turbulent flow [4]. Other important characteristic is the unsteady feature of LES. This implies that a statistical treatment is needed in order to permit an accurate characterization of the simulated turbulent flow. In addition, several theoretical studies on small-scale two-dimensional nonlocal turbulence, where the interactions of small scales with the large vortices dominate in the small-scale dynamics, have been developed in the current literature [5]. Bouris and Bergeles [6] have found that the two-dimensional large eddy simulation using a fine grid resolution, especially in the near wall region, gives a good representation of the quasi-two-dimensional mechanisms of the flow since they are directly simulated instead of being modeled as with statistical turbulence models. In addition, the two-dimensional LES performed by them has proven to be much better than any of the Reynolds-averaged Navier-Stokes (RANS) models when the major two-dimensional mechanisms of the flow and the statistical turbulence quantities are examined.
In a general case, a formal statistical treatment is based on probability averages evaluated over an ensemble of several realizations of the same process, which defines a stochastic set. For an ergodic process, the probability average can be replaced by a temporal average, and the statistical analysis is more feasible. Nevertheless, when the turbulence is dominated by large and coherent structures, typically strongly correlated, the ergodic hypothesis cannot be assumed and only a probability or statistical average (i.e., ensemble averages) should be used to describe correctly the statistical quantities of the flow [7,8].
In an LES context, the total time of simulation needs to be long enough to ensure the ergodicity of the process and to get converging statistics.
The main goal of this paper is to perform a statistical treatment of turbulent velocity signals resulting from numerical simulations, in particular, large eddy numerical simulations (LESs). The large eddy simulations are performed for the limit of high Reynolds number (Re) and low Mach number (Ma) compressible flows. The scalings show the relative importance of the subgrid terms when the flow obeys the Re 1 and Ma 1 limits. A scaling analysis is also developed in order to estimate the deviation ε between the time average and probability average associated with the ergodicity hypothesis. In addition, from the turbulent flow over a backward-facing step simulated a long time, behavior analysis is carried out in order to quantify the integral scales for ten different positions on the flow domain. Based on this correlation time, a stochastic set is built and a statistical analysis is performed. The velocity time series of the flow are analyzed statistically using T. F. Oliveira et al. 3 the formal probability average approach and confronted with the statistics given by the conventional time average analysis. The deviation between the two approaches is characterized by the ergodic parameter ε. We look at how long of a time average is necessary to obtain well-converged statistical results, and evaluate the mean-square difference between the time average and an ensemble average as the measure of convergence. Certainly this is an interesting question since the validity of the ergodic hypothesis is implicitly assumed in every turbulent flow simulation and its analysis. Turbulent intensities, skewness and flatness factors are also examined. All statistical quantities investigated are calculated using the probability average approach and the associated error bars are always evaluated. More recently, an extension of the ideas and of the method explored in this paper has been investigated experimentally for a three-dimensional flow [9].

Average equations.
Let us consider a generic flow property that can be a function of space and time φ(x,t). The spatial average φ(x,t) is defined as where x is the position vector, r is the displacement vector regarding x, and Ω denotes the volume of r-space over which the integral is taken (i.e., space average). The function G(x − r) is a filter function G : R 3 → [0,1] and satisfies This averaging process still regards the linearity and the commutability with the spatial and temporal derivatives where s = x,t. The properties (2.3) are derived from the continuity of φ and the properties of the filter function presented in (2.2) [10]. Now, a density-weighted average process is more appropriat for compressible models. This process corresponds to the well-known Favre filtered [11], defined as where ρ is the density of the fluid. Note that according to (2.4), ρφ = ρ φ. This identity is largely applied for averaging the governing equations shown below. The averaged mass and momentum balance equations written using index notation in the three-dimensional 4 Mathematical Problems in Engineering Euclidian space are given, respectively, by Here, u i are the components of the velocity vector u, p is the pressure, μ is the dynamical viscosity coefficient, D i j are the components of the strain-rate tensor D, and the symbol δ i j (the Kroenecker delta) denotes the components of the identity tensor. Now, since the velocity covariance is defined as σ u = u i u j − u i u j , the averaged product of velocities that appears in the second term on the left-hand side of (2.6) can be written as u i u j = u i u j + σ u . So, the averaged momentum equation written in the indicial notation becomes Defining the tensor Σ i j ≡ −ρσ u , a modified Cauchy equation can be written as follows: Note that we have used the Favre filtered velocity for the material derivative, D/Dt = ∂/∂t + u · ∇, namely ρDu/Dt = ρD u/Dt, with the constitutive equation for the stress tensor given by where I is the identity tensor. In this formulation, Σ represents the momentum transport by velocities fluctuations in the subgrid scales. The same averaging process is applied to the energy equation leading to where q i are the components of the heat flux vector q, given by the Fourier constitutive equation q = −k D ∇ . Here, denotes temperature, k D is the thermal conductivity, and e T is the total energy, namely, e T = e + u · u/2, with e being the internal energy. The terms T. F. Oliveira et al. 5 on the right-hand side of (2.11) identified by II and III represent the work done by the shear stress in the subgrid scale, whereas term I is the convective transport of total energy. Term I can be decomposed into two contributions as expressed below: (2.12) Using the equation of state for perfect gases, terms II and IV are directly related by the following expression: In (2.13), γ = c p /c v , where c p and c v are the specific heat at constant pressure and volume, respectively. Adding terms II and IV, the vector Q j that represents the transport of internal energy in the subgrid scales is defined, namely, (2.14) For instance, the terms III, V and the vector Q j , resulting from the averaging process of energy equation, need to be modeled in a general case.

Scaling analysis.
Before discussing a model for the subgrid terms, some scaling analysis will be performed in order to evaluate the relative importance of each subgrid contribution. The scales of the large eddies are set by the geometry and the speed of the stirring mechanism, while cut-off scales of the small eddies are determined by the action of viscosity. Here, one concentrates on the small scales for a flow with large eddies of given velocity, length, and time scales U , , T . The important Kolmogorov microscale for the smallest eddies is based on a further assumption that the smallest eddies depend only on the rate at which energy is put into the large eddies, that is, on one particular combination of U , . The friction only acts on the smallest scale and the energy is supplied only at the large scale. The rate of dissipation = 2νD : D is measured per unit of mass, and can be related to the macroscales by assuming that a significant fraction of the kinetic energy per unit of mass k = (1/2)u · u in the large eddies is dissipated in the turnover time of the large eddies, that is, per unit time, Now, the dimensions of this dissipations per unit mass are L 2 T −3 , while the dimensions of the kinematic viscosity ν are L 2 T −1 . Hence, by a simple dimensional analysis, we obtain the velocity, length, time, and strain-rate scalings of the Kolmogorov microscale: U k = (ν ) 1/4 , k = (ν 3 / ) 1/4 , T k = (ν/ ) 1/2 , and S k = ( /ν) 1/2 [12]. Introducing the Reynolds number of the large-scale eddies Re = U /ν, one obtains Mathematical Problems in Engineering The Kolmogorov microscale for the smallest eddies depends on the velocity and length scales of the large eddies in the combination = U / . Note that the turnover time of the smallest eddies T k is shorter than the turnover time of the large eddies T by the factor Re −1/2 . Hence, mixing takes place faster and more efficiently on small scales than on large scales. Large-scale mixing, however, is described by the Taylor diffusivity D T = U . So, for a container of height H, the time for eddy diffusion is then H 2 /D T = T H 2 / 2 [12].
A second microscale, the Taylor microscale, uses a different combination to yield a slightly large scale. The Taylor microscale λ T can be thought of as the boundary layer thickness on the edge of a large eddy, that is, with t being the turnover time of the large eddies T = /U . Hence, using the Reynolds number of the large eddies we can show that λ T = (ν /U ) 1/2 = Re −1/2 . Now, we turn to the scaling of the turbulence energy equation (2.11). First, note that the vector Q j can be expressed in terms of the temperature in the subgrid scale, namely, This diffusion mechanism is promoted by the velocity fluctuation transport in this scale. In such case, a typical scale of these velocity fluctuations is given by U = √ u · u , where u = u − u. The temperature fluctuations scale like U 2 . Then, a typical scale for the vector Q j is given by where R is the gas constant given by Carnot's relation R = c p − c v . The important velocity gradient occurs at the smallest scales as mentioned above. Therefore, S i j ∼ ( /ν) 1/2 ∼ (U k / k ) ∼ Re 1/2 U / , and consequently It should be important to note that using the Taylor microscale to evaluate S i j ∼ U /λ T leads to the same result given in (2.20) since λ T = Re −1/2 as mentioned before, and again S i j ∼ Re 1/2 U / . The scaling of the term V in terms of the large eddies is given by T. F. Oliveira et al. 7 Now, it is possible to determine the scaling for the ratio of the terms III/Q j and V/Q j , respectively, where Re = U /ν is the Reynolds number, Ma = U /c is the Mach number, and c is the speed of the sound. The resulting scaling indicates that for high Reynolds and low Mach numbers typically for the values of these parameters investigated in the present paper, the work done by shear stress and the kinetic energy transport done by subgrid eddies are very small in comparison to the transport of internal energy Q j . The scalings are supported by Knight et al. [13], who have evaluated subgrid terms directly by their numerical simulation. After this dimensional analysis, we can say that the only two terms to be modeled for the limit Re 1 and Ma 1 are the subgrid stress tensor Σ i j and the subgrid internal energy transport vector Q j .

Constitutive relations for the remaining subgrid terms.
The constitutive equation used to describe the subgrid stress tensor is the well-known Smagorinsky model [14], where the nonlocal turbulent viscosity μ t is calculated under conditions of inertial equilibrium subrange of turbulence [12], namely, Here γ is the average shear rate defined as γ = (2 D : D) 1/2 ∼ /k. The filter width Δλ ∼ k 3/2 / is set equal to 2δ g , where δ g is the grid spacing. It is an indication that the smallest eddies are represented by two grid points. Note that μ t ∼ k 2 / . The factor C S is known as Smagorinsky's constant. Several values have been proposed for this constant ranging from 0.1 to 0.2 [15,16]. In the present work, one has used C S = 0.20, as suggested by Deardorff [16]. Thus, the model for the subgrid stress tensor takes the form The subgrid internal energy transport tensor Q j is related to the diffusion of temperature in the subgrid scales due to velocities fluctuations and may be modeled as being a diffusive heat transport given by a nonlocal Fourier law in the form  [17], For the edge of a turbulent boundary layer, Pr t = 0.6 [18] and this value has been used in the simulations. The set of governing equations is made nondimensionalized by using a characteristic length and velocity and U , respectively, and the properties of the nondisturbed flow. From this point through all over the work, we will omit any superscript notation and assume that all properties are dimensionless averaged quantities. The set of dimensionless governing equations simulated consists of the continuity, written such as in (2.5), and the momentum and energy-averaged equations given in dimensionless terms by where Pr is the Prandtl number, Pr = c p μ/k D .

Statistical analysis
As mentioned before, the main goal of this work is to treat statistically turbulent velocity signals either from numerical simulations or experimental observations. An important question addressed here is to look at how long of a time average is necessary to obtain well-converged statistical results. Indeed, we have looked at the difference between the time average and an ensemble average as the measure of this convergence. The flow is considered as a stochastic process given by the family of functions u = u(t,α), where α = 1, ...,N are the realizations of the process according to Figure 3.1, and in the present context, u = u(t,α) denotes the velocity of the flow. By stationary we mean that the form of the probability distribution functions does not depend on a shift of the time origin. More precisely, we say that a random process is stationary when the probability distribution of the stochastic processes u(t,α) and u(t o + t,α) are the same for any t o . For a stationary random process then, we may, in principle at least, determine the various probability distributions from the observations of u(t) for one realization of the system over a long period of time T. This time being much longer than the integral scale Θ (i.e., velocity fluctuation correlation time). This long-time record can be cut up into pieces of length T λ (where T λ is much longer than any periodicities occurring in the process), and these pieces may be treated as observations of different realizations of the system in an ensemble of similarly prepared systems. We will restrict the discussion to fluctuating quantities that are statistically steady, so that their mean values and variance are not function of time. Only under this condition does the idea of a time average make sense. In a T. F. Oliveira et al. 9  typical laboratory situation, the fluctuating velocity u(t,α) should be the streamwise velocity component measured in a wind tunnel behind a cube. In particular, the relative amount of time that u(t,α) spends at various levels is measured.
The underlying assumption here is the so-called ergodic hypothesis which states that for a stationary random process, a large number of observations made on a single system at N arbitrary instants of time have the same statistical properties as observing N arbitrarily chosen systems at the same time from an ensemble of similar systems. In dealing with general stochastic process, there are two types of mean values that can be evaluated. One is the probability average obtained by a number sufficiently larger (N) of observations at some fixed time t, denoting this average by u(t) , and the other is the time average made for a function u(t), denoting this average by u. The requirement that a time average should converge to a mean value, that is, that the error should become smaller as the integration time T increases, and that the mean value found this way should always be the same, is the ergodicity [19]. This point is discussed in more details in Section 3.2.
In the case of an ergodic stationary random process, both averages yield the same result, provided that the function u is finite and continuous in mean square [20]. The time average u over a sufficiently long realization α o (i.e., for time much longer than the velocity fluctuation correlation time) of the flow is defined as [21] The use of time averages corresponds to the typical laboratory situation, in which measurements are taken at fixed locations in a statistically steady, but often inhomogeneous, flow field. For a time average to make sense, the integral (3.1) has to be independent of t o . In other words, the mean flow has to be stationary ∂u/∂t = 0, and consequently the mean value of the velocity fluctuations, u = u(t) − u, itself is zero by definition. Here, the averaging time T needed to measure mean values depends on the accuracy desired as discussed in Section 3.2. Now, whether each realization has the same probability to occur, a statistical (or probability) average is defined as being [22] A fluctuation about the probability average is defined as being The variance (or the turbulent kinetic energy (1/2) u (t) 2 ) is calculated from the probability average as being the probability average of the square of the velocity fluctuation. While the time average of the square of the velocity fluctuations is given by

Correlation function and the spectral density.
The random processes that do occur often in flow applications are those where u (t) and u (t ) will be correlated at least for small values of τ = |t − t|. There are two more functions associated with a continuous stationary random process that is central to a statistical description. These two functions are the velocity fluctuation autocorrelation function and the spectral density. The normalized velocity fluctuation autocorrelation function of a continuous stationary random process is defined as [23] Now, since a shift in the origin of time does not affect any of the statistical properties of a stationary random process, the probability density functions simplify from f (u ,t) and g(u t ,t;u t ,t ) to f (u ) and g(u t ,u t ,t − t ), respectively. Here, f (u ,t)du is called the first probability distribution and g(u t ,t;u t ,t )du t du t , the second probability distribution, is the joint probability of finding u (α,t) between u t and u t + du t at time t and between u t and u t + du t at time t . In this particular case, when both the probability average and the normalized autocorrelation function do not vary with a shift in the origin of time, we simply write u , u 2 , and R(τ) and the process is said to be statistically stationary, with the time shift The other central function for the statistical analysis here is the spectral density of u . It is well known that the variance of the process u 2 corresponds to the average power dissipated in the interval (−T,T). Then T. F. Oliveira et al. 11 where E(ω) is the power spectrum of u (t). Thus, E(ω)dω/2π is the average power dissipated with frequencies between ω and ω + dω. The velocity fluctuation autocorrelation function and the spectral density are connected by the Wiener-Khintchine theorem which states that [7] Note that the autocorrelation function is related with the autocorrelation coefficient by the variance, that is, C τ (τ) = u 2 R(τ). Thus, according to the above theorem, the correlation function and the spectral density are simply Fourier transforms of each other.

Ergodicity.
As stated before in this work, we look at how long of a time average is necessary to obtain well-converged statistical results. For this end, we need to evaluate the mean-square difference between the time average and an ensemble average as the measure of convergence. This is an interesting question since the validity of the ergodic hypothesis is implicitly assumed in every turbulent flow simulation and its analysis. Thus, using the definition of correlation function given in (3.5), one obtains [7,24] is an important result that relates the correlation function with the variance σ 2 . Note that if T → ∞ leads to σ 2 → 0, this implies the ergodicity condition. Thus, the mean value of the fluctuating quantity can be determined by a time average with accuracy defined by the size of the integral time scale Θ. In fact, this is the requirement that a time average should converge to a mean value with the error becoming smaller as the integration time T increases. An ergodic variable not only becomes uncorrelated with itself at large time step, but it also becomes statistically independent of itself. An equation like (3.8) may be also used to evaluate the mean-square error of the difference between the average value of u(t) in the laboratory (evolving finite integration time) and the true mean value (requiring integration over an infinitely long time). Usually, in a process in which the correlation function decays rapidly for a relatively short time τ, the ergodic condition is verified. In particular, for a turbulence in which the correlation function has the same exponential decay of the one corresponding to a random walk process, we have where Θ is a correlation time associated to an interval in which the events are weakly correlated. For large time intervals compared to Θ, the flow u(t) becomes statistically independent of itself, so that Θ is a measure for the time interval over which u(t) remembers its past history. In a Markovian diffusion process, the variance in the displacement from the starting position increases linearly in time. The coefficient of the linear growth is defined to be the diffusivity of the random walk Γ = (1/2)d/dt( x 2 (t) ). Again, the symbol · denotes an average over several experiments. Now, proceeding formally, Taylor's calculation of the eddy diffusivity in turbulence is given by The diffusivity attains its constant value only after several correlation times. Thus, Γ = u 2 Θ, with u 2 being the mean square of the velocity fluctuations and Θ being the integral time scale, By the present analysis, we can show that when the above integral fails to converge, due to the slow decay of the velocity fluctuations autocorrelation, the diffusion becomes anomalous with x 2 (t) ∼ t n , with n = 1 [25,26]. This behavior is characteristic in large-scale region of the turbulent flow investigated where these scales are strongly correlated, and the turbulence does not loose its memory. It is instructive to remind the reader that the velocity correlation needs to be computed as seen by a particle moving with the fluid, and hence Θ is called the Lagrangian integral-correlation time. The length scale introduced in Section 2 as the size of the eddies can now be defined as = u Θ, that is, the distance one would move at ( u 2 ) 1/2 (i.e., the root mean square of the velocity fluctuation or simply the RMS) during the correlation time Θ. Now, the magnitude of the error associated with the nonergodicity (i.e., the convergence of time average) of the process is defined as being Considering an exponential decay for the correlation function, the integral in (3.8) may be performed to give an estimation of ε, namely, where I = u 2 / u is the turbulence intensity of the flow that is a measure of the relative importance between the velocity fluctuation and the mean flow. It is clear from (3.13) that the convergence of the time average to a mean value can be determined to any accuracy desired if the integral scale Θ is finite. In particular, (3.13) gives a good estimation of the long time T needed to verify the ergodicity of the flow. The time average should also be used as a correct approach to describe the process from a statistical point of view. In this work, the result expressed in (3.13) is tested by direct evaluations of the variance σ 2 .
As mentioned before, an important quantity to quantify flow memory is the Taylor time scale [23]. Using a Taylor series to expand u (t + τ) in a neighborhood of t and supposing the process to be stationary, the correlation function based on an ensemble T. F. Oliveira et al. 13 average may be written as (3.15) From (2.17) and (3.15), it is clear that λ T is a short time scale of the correlation process. So, using a second-degree polynomial function in order to fit the correlation function for short times, λ T may be estimated. Typically, the Taylor scale is larger than a dissipative time scale, but is not related to the integral scale observed in the macroscopic flow, that is, 2 k /ν λ T /U , where k is the Kolmogorov dissipative length scale as defined in Section 2. Effectively, the Taylor scale is a memory characteristic time of the flow. If t is the present time, we can say that the flow has a strong dependence on the events that occur in the interval (t − λ T ,t).

Numerical simulations
We now briefly summarize the sequence of steps that is necessary to perform our numerical simulations. Large eddy simulations admittedly require denser grids and more computer time than Reynolds averaged approaches, but in certain cases, such as the one under consideration here, it seems that the Reynolds averaged approaches fail to predict important statistical aspects of the flow. The purpose of this section is to provide an accurate simulation of the turbulent flow past a backward-facing step using a large eddy simulation in two dimensions. We have therefore obtained information on the statistics of the flow from the velocity time series generated from the simulations according to the procedure described in Section 3. We will show that certain region of the flow corresponding to the formation of coherent large-scale structure cannot be described by the well-converging time average approach.
The flow investigated was numerically simulated under a two-dimensional large eddy fashion. We use a finite-volume method on a two-dimensional Eulerian grid to solve hydrodynamic and energy equations of the flow in Cartesian coordinate for a structured mesh with colocalized variables. The governing equations are solved simultaneously. The Reynolds and Mach numbers based on the step height were Re = 38000 and Ma = 0.03, respectively. Figure 4.1 shows the flow domain and the location of the analyzed points in the flow and typical streamlines of the mean turbulent field. The streamwise and spanwise lengths of the computational domain were 20H and 2.5H, respectively. The filtered governing equations described in Section 2 were discretized by using the explicit Mac-Cormack method written for a finite volume formulation [27]. The numerical method uses a standard explicit predictor-corrector Euler algorithm to carry the temporal march. The set of governing equations was discretized using forward first-order differences in the predictor steps and backward first order differences in the corrector steps. In both 14 Mathematical Problems in Engineering steps, the strain rate and the temperature gradient components were evaluated at the finite volume faces by central differences. We summarize next the basic procedure of the discretization process.

Discrete approximation of the balance equations.
The general framework for the filtered balance equations described in Section 2 can be written in the form of the following vector equation: where Π = E e 1 + F e 2 , and the vectors U, E, and F are defined in the same way as given by Anderson et al. [28], Here the subscripts 1 and 2 in the above vectors denote the components x 1 and x 2 of a flow quantity, e 1 and e 2 are the unit basis vectors in directions 1 and 2, respectively, and the parameters C 1 = 2Re −1 and C 2 = [(γ − 1)PrMa 2 Re] −1 . Now, consider that the flow domain Ω is subdivided in N small regions called control volumes, so that Ω = N i=1 i . The volume average of the quantity U over a single control volume is defined as being U = (1/) UdS. Now, integrating (4.1) over and using the divergence theorem with the volume average definition, it is found that where Ꮿ is the contour of the elementary region and n is a unit vector directed outward from the enclosed control volume. The flow domain Ω was discretized by using quadrangular control volumes in a structural mesh as illustrated in Figure 4.2. The surface integral (4.3) is evaluated by line integrals over the edges of the control volumes. The value of Π on an elementary volume edge is taken as being the volume average of Π over one of the control volumes bounded by this edge. In particular, it is considered that Π is always constant with respect to the integral over Ꮿ, so that Ꮿ Π · ndᏯ ≈ 4 β=1 Π β · β , where β is a vector normal to the edge β with absolute value equal to the length of this edge. Under this condition, (4.3) reduces to the following approximation: Since Π is a function of U, (4.4) can be solved by an Euler method. In the MacCormack method, a predictor-corrector algorithm corresponding to a second-order Runge-Kutta procedure is applied. In the predictor step, the vector Π was taken as being equal to the vector of the control volume downstream to the edge, whereas in the corrector step, Π was associated to the volume upstream to the same edge. For instance, following the notation in Figure 4.2, consider the discrete quantity evaluated on the edge i + 1/2. At the predictor step, Π is calculated in terms of the quantities of the volume (i + 1, j). Subsequently, for the same edge, at the corrector step, Π is calculated in terms of the quantities of the volume (i, j). The components of the velocity gradient tensor and the temperature gradient on the faces of the control volume were evaluated by central difference. For instance, on the edge i + 1/2, the term ∂u 1 /∂x 1 is approximated by where x 1 (i, j) and x 1 (i + 1, j) denote the center coordinates of the control volumes (i, j) and (i + 1, j), respectively. The compressible formulation of the flow allows the use of a state equation for the pressure. Consequently, additional velocity-pressure coupling algorithm was not required. This numerical procedure leads to a second-order precision discretization in both space and time derivatives. The accuracy and robustness of the numerical method are defined in the context of the classical MacCormack method described in details by [27]. No extra upwinding feature was implemented and the method is stable if the CFL number (i.e., Courant-Friedricks-Lewey number) is less than unity for all grid volumes, where CFL = Δx/(u + c s ) with u and c s being the largest velocity norm and the largest sound velocity at the volume boundary, respectively, and Δx is the local grid spacing [28]. In the present simulation, CFL ≈ 0.7, which has provided a stable condition for all simulations. The initial transient, corresponding to 608 flows through (i.e., 10 6 iterations), that is, 2.3 seconds of physical time, was neglected. The typical convective time scale of the flow H/U ∞ ≈ 0.1 second and the longest simulation time was approximately 21 seconds. This wide interval was necessary because the velocity signal was fragmented into smaller temporal series when defining the stochastic set.
An equally spaced Cartesian grid was used to discretize the flow domain and to resolve the velocity, pressure, and temperature fields of the flow. The spatial resolution employed was 36 volumes/cube edge. This resolution was sufficient to resolve the turbulent length scales of interest. Additionally, a 100H stretched grid region was generated with 188 volumes in the streamwise direction, downstream of the regular grid. This region was necessary to dissipate the turbulent structures and provide a smooth condition at the outgoing section of the domain. The computations were carried out for a typical rectangular grid with (90 × 1088) control volumes. A uniform velocity profile at the inlet section was imposed with no turbulence intensity (i.e., laminar flow). No slip boundary conditions were employed in the spanwise and normal directions. The initial condition is formed by stagnated fluid with constant pressure and temperature. The dimensionless time step used was approximately Δt = 5 × 10 −3 . The simulation time was T = 48800H/U ∞ . The velocity samples at probe positions were stored in intervals of 25Δt. In postprocessing algorithm, the velocity time series have been used to compute the various statistical results to be presented next. The time separation between two sequential data fields was large enough compared to the integral time scale of the turbulent fluctuation Θ for the data fields to be considered are nearly independent realizations of the flow. The flow statistics have been obtained by or ensemble averaging over all stored data fields. The simulations all were carried out on a PC of 2.0 GHz processor and 1.0 GB of physical memory. The total CPU time required to perform the 2D LES simulation was about 27 days. region, the velocity signals were collected in order to build the stochastic set to be analyzed. In addition, we can see the flow separation and that a reattachment occurs more suddenly. The converged solutions were checked against the experimental results for the average reattachment point position given by Eaton and Johnston [29]. This is a typical parameter that has been commonly used to validate numerical simulations. The test simulation was carried out on the same conditions of the experimental setup, including flow domain geometry, thermodynamic properties of the fluid, and the imposed flow.  [29]. The error was less than 1%. Halving the grid size produced a change that was not greater than 1% in this computed quantity. The first step for the statistical characterization of the flow was to define the stochastic set in the probes positions. In order to build a stochastic process from the numerical simulations, a large temporal series was dropped into smaller temporal series corresponding to the realizations of the process. The resulting temporal series were then independent events of the turbulent flow, since the time scales involved were long enough for a complete decay of the correlation function. So, the velocity fluctuations become statistically independent with respect to events that have occurred in their past history. This procedure was equivalent to start a new simulation from a different initial condition which is uncorrelated with the past one. The dimensionless integral time scale ΘU ∞ /H may be determined by direct integration of the velocity fluctuation normalized autocorrelation function (3.11). In this work, however, we have estimated this parameter by direct inspection of the velocity fluctuation autocorrelation function as being approximately two times the time for the full decay of the autocorrelation function. This of course overestimates the integral scale value of the correlation time, but it corresponds to a suitable time scale which guaranteed the statistical independence of the set of realizations constructed from the original velocity signal.  Table 5.1 shows the length of these intervals ΘU ∞ /H and the number of realizations associated with each probe. By this procedure, the statistics of the flow were performed in terms of the probability averages. large-scale structures located in the transition shear layer region of the flow. Actually, the probes 4 and 5 are located in a place which appears to be more critical for the ergodicity than probes 1, 2, 3 that are exactly inside the recirculating region. Inside this region, the fluid parcels seem to describe approximately a rigid body rotation since the important mechanism known as vortex stretching is absent in two-dimensional flows. The probes 4 and 5 (mainly probe 4) on the other hand are located in a transition shear layer region where the fluid particles are subject to stronger velocity gradients, and the complex interactions of smaller scales with intermediate scales cannot be neglected. The coherence of the turbulent eddies in the shear layer is much increased if no three-dimensional instabilities are present (as it is the present case). The intensity of the coherent vortices grows via an inverse energy cascade, and eventually they start producing significant feedback on turbulence. Despite many achievements in numerical simulations of two-dimensional turbulence, the underlying physical mechanism of turbulent vortex interactions still remains unclear.

Results and discussions
The mean-square difference between the probability average and the time average gives a direct measurement of the error , defined in (3.12). This error can also be estimated by the simple relation proposed in (3.13). Table 5.2 shows the values of the predicted and the ergodic deviation errors for each probe in the turbulent flow. The purpose of relation (3.13) is to give an estimation of the order of magnitude of the error introduced when a probability average is replaced by a temporal average. The results presented in Table 5.2 show a good agreement between the scaling based on an exponential decay of the autocorrelation function and the ergodic deviation error calculated numerically by using (3.12) with the autocorrelation sample computed from the numerical simulations. It is seen that for those points far outside from the recirculating bubble the turbulence behaves close to a random walk. On the other hand, the results also indicates that for those points inside and around the recirculating bubble, the ergodic deviation can be very high such as that found for the probe 4. So, the time average does not produce a meaningful statistics.
For the case of the probes 4 and 5, an exponential decay does not fit the real decay behavior of the normalized correlation function. It indicates that the turbulence in this region may have a quite different behavior of a typical random walk process. The dispersion process of momentum transport by velocity fluctuations seems to characterize an anomalous diffusion in the way described in Section 3. In that case, the integral in (3.8) must be evaluated numerically. The probe 4 shows a strong nonergodic property, suggesting that in this region a time average approach fails in describing the local turbulence. It is possible to infer that the probes in the neighboring or inserted into the recirculation bubble shown in Figure 4.1 have presented a significant deviation from the ergodicity. Consequently, the flow in this regions persists strongly correlated for a long time as shown in Table 5 the turbulence memory, corresponding to the time over which the velocity fluctuation is correlated with itself (i.e., the velocity fluctuation correlation time).
In addition, the two-dimensional power spectra relative to the probes 4 and 8 are presented in Figures 5.4(a) and 5.4(b), respectively. The details of the spectra calculations were mentioned in Section 3.1. The error bars are also considered in the plots. In these T. F. Oliveira et al. 23 plots, the abscissa is the logarithm of the nondimensional frequency, whereas the ordinate is defined so that the area beneath a logarithmic plot of E(ω) is proportional to the mean square of the fluctuating signal. Both spectra show that turbulence energy at small scales is decreased, while it is increased at large scales. In particular, one can see that the spectra of present two-dimensional simulations seem to change shape for moderate nondimensional frequencies ranging form 0.1 to 1. In the case of probe 4 located in the interface of the recirculating bubble as shown in Figure 4.1, we see that the energy cascade is characterized by a −3 spectral exponent which is different from the famous ω −5/3 valid for three-dimensional small scales of a local turbulence. The decay turbulence with ω −5/3 is observed only at higher frequency (smaller scales) as a result of the subgrid model used in numerical simulations which has been based on the Kolmogorov inertial equilibrium subrange of turbulence described in Section 2.3. The interval of the ω −5/3 spectrum is shown in both inserts of the mentioned figures. The presence of a spectrum interval with an ω −3 decay may be attributed to the mechanism of the inverse cascade which arises in two-dimensional turbulence inside the recirculating bubble. The resulting spectrum of probe 8 located outside the recirculating bubble ( Figure 5.4(b)) exhibits virtually the same characteristics shown in the plot in Figure 5.4(a). However, the decay turbulence given by this spectrum is closer to ω −4 than to ω −3 . So, outside the recirculating region it is found that the turbulence decay corresponding to intermediate scales with dimensionless frequencies ranging form 1 to 10 gets steeper than −3. Most of these results seem to be in qualitatively agreement with the predictions given by the dimensional analysis of two-dimensional turbulent flow presented by Nazarenko and Laval [5]. The probability density functions associated with the turbulent velocity fluctuations in probes 4 and 8 are shown in Figures 5.5(a) and 5.5(b), respectively. The non-Gaussian deviation of the distribution function is clearly noticeable in probe 4, whereas in probe 8, the behavior of the probability density function is closer to a normal distribution. In particular, the behavior of the statistical distribution is quantified by the skewness and flatness 24 Mathematical Problems in Engineering factors, defined as ϕ = u (t) 3 /ξ 3 and κ = u (t) 4 /ξ 4 , where ξ 2 = u (t) 2 , respectively. These factors and the turbulence intensities for each probe are listed in Table 5.3. It is well known that a normal distribution has ϕ = 0 and κ = 3. It is possible to infer that all processes display some non-Gaussian behavior. The turbulent intensity at probe 4 has the order of 2 × 10 4 %, whereas between others, the greater value of this parameter has the order of 10 2 % (see, e.g., probe 1, Table 5.3). This characteristic is an indication of the strong nonergodic property of the velocity fluctuations in probe 4. The interval of confidence represented by the error bars shown in the plots and by the associated errors to the quantities in Table 5.3 is relatively large. It suggests that for a complete characterization of these parameters, more realizations are required, and consequently it would demand additional computational effort in order to simulate much larger time intervals.

Concluding remark
In this paper, a rigorous statistical approach for the treatment of turbulent velocity fluctuations has been presented. We have looked at how long of a time average is necessary to obtain well-converged statistical results. For this end, we evaluate the mean-square difference between the time average and an ensemble average as the measure of convergence. From the numerical simulations, ten different points in the flow domain have been statistically treated using a probabilistic approach. The realizations of the statistical ensemble were defined by the cut up of a long-time velocity record into pieces of length much longer than the characteristic correlation time of the velocity fluctuations. Based on the velocity fluctuations correlation time, a statistical analysis of long time has been performed. The ergodicity assumption of the turbulence was investigated. The deviation ε of this criterium was evaluated and compared with theoretical predictions given by scaling arguments. The results have suggested that the deviation due to ergodicity assumption may be used as a criterion in order to predict the upper boundary simulation time required to have a convergence of the statistical characterization of the flow. We have found from the two-dimensional large eddy simulations that inside and outside a recirculating region, the turbulence decaying corresponding to intermediate scales (with dimensionless frequencies ranging form 1 to 10) gets steeper than the classical scaling ω −5/3 spectrum. We found ω −3 and ω −4 for a point inside and another outside the recirculating bubble, respectively. The decaying turbulence like ω −5/3 was seen only at the smaller scales attributed to the subgrid model used which was based on the Kolmogorov inertial equilibrium subrange of the turbulence. Probability functions, skewness and flatness coefficients have shown a deviation from the Gaussian behavior at all investigated positions. We have seen that the break of flow ergodicity property of the flow is directly related with the large-scale structures of the turbulence. In this flow regime, the integral time scale required to define the most commonly used time average is almost always unpredictable.
T. F. Oliveira et al. 27 As future works, it would be important to think about developing more robust threedimensional LES to test the ergodicity of turbulent flows in the presence of vortexstretching and multistructure interactions. Two-dimensional LES calculations are clearly inferior to three-dimensional ones since certain important features of three-dimensional turbulence (vortex stretching) are not resolved. Moreover, two-dimensional large-scale structures are subject to three-dimensional instability which results in counter-rotation vortices. Behind this, there is a question of how the nonlocal interactions of large-scale vortices with intermediate and small-scale structures in three-dimensional turbulence affect the local ergodicity of the flow. Clearly, this fundamental topic requires further attention.