Direct Numerical Simulation in Two Dimensional Homogeneous Isotropic Turbulence Using Spectral Method

Direct numerical simulation (DNS) in two-dimensiona l homogeneous isotropic turbulence is performed by using the Spectral method at a Reyn olds number Re = 1000 on a uniformly distributed 128 128× grid points. The Reynolds number is low enough tha t the computational grid is capable of resolving all the possible turbulent scales. The statistical properties in the computed flow field show a good a greement with the qualitative behavior of decaying turbulence. The behavior of the flow st ructures in the computed flow field also follow the classical idea of the fluid flow in turb ulence.


Introduction
Understanding the structure in space of a turbulent flow as well as its statistical properties remains a challenge both for the experimentalist and the theoretician.Direct simulation, i.e. resolution of the basic fluid dynamics equations using the powerful computers has proven to be a valuable additional tool for the study of fully developed turbulence.For the range of parameters in which they are feasible, the direct simulations allow measurement of many quantities inaccessible in the laboratory.Also, visualization of the small-scale vortical structures in the computed flow fields becomes easier.Now-a-days, high resolution simulations in two and three space dimensions at Reynolds numbers of several thousand or more are possible [1][2], and have revealed new and important properties of two and three-dimensional turbulence.
Direct numerical simulation (DNS) is considered as the most exact approach to turbulence simulation but it is too expensive and a large scale of computational resources is required to carry out the DNS [3][4][5][6][7][8].If Reynolds number of flow is very high and computational grid is very large then it becomes more and more expensive and time consuming.Turbulent fluid motion is fully three-dimensional and complex.However, in computational point of view, two-dimensional (2D) turbulence is easier than the fully developed three-dimensional turbulence to compute.Discretization method is another issue to conduct the numerical simulation in turbulence.A literature review suggests that the numerical method widely used for DNS is either spectral method or the conventional finite difference method with structured grids.
The overall aim of our present research is to develop a numerical code based on spectral method in order to simulate and analyze the physics of turbulence.However, before its uses to the complicated flow fields in turbulence, it is essential to examine the performance and effectiveness through some simple and benchmark problems in turbulence.Therefore, in this study, we performed the direct numerical simulation in two dimensional homogeneous isotropic turbulence at a Reynolds number Re = 1000 with 128 128 × grid points.Our interest is in whether the present numerical solver can simulate and resolve the turbulent scales in homogeneous isotropic turbulence.The Reynolds number is low enough that the computational grid is capable of resolving all possible turbulent scales in homogeneous isotropic turbulence.We also discuss about the flow structures in the computed flow field by contour and vector plots of the flow.

Flow Governing Equations
The governing equations are the unsteady incompressible Navier-Stokes equation and the equation of continuity as follows: ( ) Here u is the velocity field, p is the pressure and Re is the Reynolds number of the flow.The Eq. ( 1) can be written in the rotational form as follows: where, , is the vorticity. ( It is noted that the rotational form of Navier-Stokes is commonly used in high Reynolds number flow simulations [9].Rotational form of Navier-Stokes equations gives better physical properties in terms of conservation laws of fluid and the result is more stable than the convective form.It is also less expensive from the computational point of view. Hence, in this study we solve Eq. ( 3) with continuity Eq. (2).

Discretization Method
The discretization of governing equation is performed in three steps of time integration based on spectral method [10].Some description of this method is given in the previous studies by Tanahashi et al. [4][5][6].It is stated earlier that we intend to develop spectral DNS code to simulate three-dimensional flow fields so the details of the spectral method will be given elsewhere.However, the three steps of time integration performed in this study are described as follows: First step: ( ) where, and ∆t is the time step increment .
Second step: , where, , where, Here the first step is for the nonlinear term performed by second order Adams-Bashforth method.The second step is for pressure adjustment performed by backward Euler method, and the third step is for viscosity term performed by Crank-Nicolson method.

Initial Conditions
It is the difficult part to produce initial flow fields for homogeneous turbulence.In this study, in the procedure to make the initial flow fields, the initial flow condition is assumed to have the decay of energy spectrum [11].We assume ( ) k φ , k is the wave number, is a vector potential and then the initial velocity field is determined (in Fourier space) by ( ) ( ) is the amplitude of the vector potential and ( ) determined by an uniform random number within the range [ ] π 2 , 0 .The amplitude ( ) is determined so that the following relation is satisfied: Here we should note that the energy does not depend on ( ) The above operation (Eq.( 7)) corresponds to determining the energy distribution within a spherical shell, whose radius is k , in the wave-number space.However, in the current study we uniformly distribute the energy as where k w N , is the number of the wave number vectors whose length is within the range of ( ) . Here the actual velocity condition for homogeneous turbulence in Fourier space ( ) ( ) should be fulfilled.Using this procedure we have calculated the initial velocity for DNS in which the decay of initial energy is fully monitored.In addition, to obtain a reliable result for DNS in homogeneous isotropic turbulence, we are concerned about two necessary conditions such as; (i) grid spacing is less than three times of Kolmogorov microscale, and (ii) dimension of the computational domain is larger than four times of integral length scale.According to these conditions the maximum possible Reynolds number (Re) of the flow is considered for this computation.

Computational Parameter
The computational domain of the mesh was selected to a periodic domain 2π×2π.The computation has been performed using 128×128 computational grids and the possible Reynolds number of the flow is 1000.The computation has been done with non dimensional time increment, ∆t = 0.0001 and is executed up to time, t = n ∆t, where n is the time step.

Turbulence Statistics
In this section we discuss some statistics in 2D homogeneous isotropic turbulence.The total resolved energy, E k versus time is presented in Fig. 1, where In Fig. 1 we can observe that the total energy or, equivalently the kinetic energy decreases with time until the end of the calculation.The trend of the profile is always decaying which is in agreement with the result of Yokota et al. [12], where the kinetic energy spectrum of a decaying homogeneous isotropic turbulence is calculated using a pure Lagrangian vortex method.Fig. 2 shows the decay of resolved enstrophy, where the enstrophy is defined by In Fig. 2, we can observe that the enstrophy decreases with increase of time, which is also in agreement with the result of Yokota et al. [12] for a decaying homogeneous isotropic turbulence.The turbulence dissipation, ε versus time is presented in Fig. 3, where the dissipation is defined by where, Like the total energy and enstrophy given in Figs. 1 and 2, the dissipation also decreases with time, and with the increase of time of the computation the dissipation tends to be zero, that is, we will get dissipated flow fields.However, observing the profiles of total energy, entropy and dissipation in Figs. 1,  The decay of root mean square (r.m.s.) of velocity fluctuations is presented in Fig. 4, where r.m.s. is calculated by using the definition given as It is revealed that the r.m.s. of velocity fluctuation decreases with the increase of time and it reaches about 0.80 at the end of the calculation (t = 0.30).The development of Kolmogorov microscale, η is presented in Fig. 5, where η is defined as In Fig. 5, we can observe that the Kolmogorov microscale increases with the increase of time.At the end of the calculation it reaches a value of about 0.0066.Since the computation is performed for small time and also in 2D, so we can say that the calculated result in our study is reasonable.
Fig. 6 shows the development of integral length scale, l which is defined as where u(x) is the two-point correlation co-efficient.The integration is evaluated using Simpson's 3/8 rule in the space of computation.It is revealed that the profile of integral scale turns up from the initial stage to near about middle stage then it turns down up to final stage of the computation.At the end of the computation the result of integral scale reaches about 0.19 (t = 0.30).The time dependence of Taylor microscale is presented in Fig. 7, where Taylor microscale is defined as: Fig. 7 shows that the profile of Taylor microscale in 2D homogeneous isotropic turbulence increases with the increase of time.It is revealed that in the initial stage the increasing rate is very slow and tends to be constant.Then the increasing rate turns up gradually up to final stage and finally it reaches about 0.10 at t = 0.30.Since the computation in our study is performed for small time and 2D homogeneous turbulence, so our calculated result is reasonable.
The development of Enstrophy based length scale, l d is presented in Fig. 8, which is defined as , where In Fig. 8, we can observe that the enstrophy based length scale tends to be constant at the initial stage.Then it decreases with time up to the middle stage and after that it increases gradually with time up to the final stage.At the end of the calculation it reaches about 0.0105.Like the other statistical values this result is also considerable for this computation in the homogeneous isotropic turbulence.
Skewness and flatness factors of velocity and their derivatives are important statistical properties that represent characteristics of turbulence.The production of the rate of dissipation of turbulent kinetic energy, or equivalently, the production of enstrophy is directly related to skewness in isotropic turbulence [13].Skewness and flatness of a velocity component, u i are defined as follows: Fig. 9 shows the profile of skewness of velocity component u.At the beginning of the computation it decreases, and then for few time steps it is almost constant.Then it increases rapidly within a short time steps and finally, it decreases till the end of the computation (t = 0.30).The development of derivative of skewness of velocity component u is given in Fig. 11 which is about 0.08 at the beginning of the computation (t = 0) and it reaches about 0.32 at the end of calculation.The behavior of the derivative of the skewness in the whole calculation domain is in good agreement with the result of Yokota et al. [12].The flatness and derivative of flatness of velocity component u (Figs. 10 and 12) are almost 3.4 at the beginning of the computation and shows almost similar pattern up to the final stage of the calculation.Since it is well established by the researchers that the flatness should reach the values somewhere in between 3 to 4, so our computed results in 2D simulation show an excellent agreement with them.

Flow structures
We have calculated the velocity components, vorticity and pressure at non-dimensional time t = 0.30 by using spectral method, and using these computed data we have shown different contour and vector plots of the flow field.We have also investigated the relation among velocity, vorticity and pressure by examining their flow structures.
We know from the properties of fluid flow that there is a relation among velocity, vorticity and pressure.At an instantaneous time, where the velocity is positive (i.e.anticlockwise), there the vorticity is positive (ω > 0) and the pressure is low (p < 0).On the other hand, where the velocity is negative (i.e.clockwise), there the vorticity is negative (ω < 0) and a high pressure (p > 0) should exist.
In Fig. 13, we observe the behavior of velocity, vorticity and pressure of the flow field at a non-dimensional time t = 0.30.We can see from these figures that, where the velocity is positive (i.e.anticlockwise) (Fig. 13(a)), the vorticity is positive (Fig. 13(b)) and the pressure is low (Fig. 13(c)) in the corresponding region.On the other hand, where the velocity is negative (i.e.clockwise), the vorticity is negative and high pressure exist in the corresponding region.
The highest value of ω and p appear at red region and lowest value at blue region in the corresponding contour plots.Whereas, ω takes a value in the range 11.92 ≤ ω ≤ 80.7133 in most of the regions and p takes a value in the range 0 ≤ p ≤ 1 in most regions of the corresponding contour plots.
By examining the vector plot of velocity, contour plot of vorticity and pressure field, we can observe that the computed results agree well with the classical idea of the fluid flow in turbulence.Fig. 13(a).Vector plot of velocity at t = 0.30.Fig. 13(b).Contour plot of vorticity (ω) at t = 0.30.Fig. 13(c).Contour plot of pressure at t = 0.30.By examining the vector plot of velocity, contour plot of vorticity and pressure field, we can observe that the computed results agree well with the classical idea of the fluid flow in turbulence.

Conclusion
Direct numerical simulation in two dimensional homogeneous isotropic turbulence have been successfully performed by using spectral method at a Reynolds number, Re = 1000 with 128 128 × grid points.To our knowledge with this Reynolds number our selected computational grid points hopefully resolve possible turbulent scales.The statistical properties such as energy, enstrophy, dissipation, root mean square of the velocity decrease with the increase of time which reveal that the computed flow fields agree well with the qualitative behavior of decaying turbulence.From the development of kolmogorov microscale, Taylor microscale, integral length scale, enstrophy based microscale, skewness, flatness and their derivatives, it is also revealed that the computed flow fields agree well with the qualitative behavior of fully developed and decaying turbulence.By examining the flow structures in our computation, we have observed that our results also maintain the relation among the fluid velocity, vorticity and pressure.
2 and 3, we can conclude that our DNS t produce a real and fully developed turbulence fields at the Reynolds number Re = 1000. u