## Abstract

Two-dimensional decaying turbulence in a square container has been simulated using the lattice Boltzmann method. The probability density function (PDF) of the vorticity and the particle distribution functions have been determined at various regions of the domain. It is shown that, after the initial stage of decay, the regional area averaged enstrophy fluctuates strongly around a mean value in time. The ratio of the regional mean and the overall enstrophies increases monotonously with increasing distance from the wall. This function shows a similar shape to the axial mean velocity profile of turbulent channel flows. The PDF of the vorticity peaks at zero and is nearly symmetric considering the statistics in the overall domain. Approaching the wall, the PDFs become skewed owing to the boundary layer.

## 1. Introduction

In the last decade, the numerical computations of decaying two-dimensional turbulence in a container have revealed some interesting results. First, Li & Montgomery [1] pointed out that the decay of two-dimensional turbulence in a circular container differs significantly from the decay process in a double-periodic domain. In particular, the enstrophy decays more slowly in wall-bounded domains because the walls act as vorticity sources. Similar observations were given by Clercx *et al.* [2], who studied decaying two-dimensional turbulence in a square container. These authors concluded that vortices approaching walls produce new oppositely signed vortices, which are then advected into the interior of the domain. Clercx & van Heijst [3], studying the one-dimensional energy spectra along the walls, found that near the walls the spectra follow a power-law *k*^{−n} with exponents *n*=−5/3 for *k*<*k*_{i} and *n*=−3 for *k*>*k*_{i} in a limited period of time during the evolution. Similar two-dimensional spectra can be observed in the case of forced two-dimensional turbulence where *k*_{i} is in direct relation to the wave number of forcing [4,5]. In a wall-bounded case, Clercx and van Heijst [3] associated *k*_{i} with an average boundary-layer thickness. The evolution of Eulerian statistics of vorticity and pressure in the case of two-dimensional decaying turbulence in a cylinder was first studied by Schneider & Farge [6]. It was pointed out that an initially Gaussian vorticity probability distribution function (PDF) peaks at zero during the evolution owing to the viscous boundary layer. In a recent study, the Lagrangian statistics of the same problem yielded the surprising conclusion that the Lagrangian boundary layer thickness extends up to the ratio of *r*_{0}/*R*=0.3, where *R* is the radius of the domain and *r*_{0} is a radius above which the influence of the boundaries on the Lagrangian statistics is important [7].

In this paper, we study regional statistics of wall-bounded two-dimensional decaying turbulence using the simulation results obtained by lattice Boltzmann simulation. In particular, we would like to learn more about the decay of turbulence in a container, since real physical experiments made, for example, in thin electrolyte layers [8] are inherently bounded by the wall of the experimental device. Therefore, a comparison of experiments with theories developed by assuming infinite domains is not straightforward at all. Obviously, before doing such a comparison one should determine the extent of the near-wall region in which the effect of walls is significant or, from another point of view, one should demonstrate that there is a central region, in which the wall effects are negligible. The outline of the paper is as follows. In §2, the numerical method used is briefly introduced; in §3, the problem description is given, including the discussion of regional statistics. In §4, our simulation results are presented and conclusions are drawn in §5.

## 2. Lattice Boltzmann method

For completeness, let us briefly recall some basic facts about the lattice Boltzmann method (for details, see the review of [9,10]), which was first applied for the simulation of two-dimensional turbulence in the pioneering work of Benzi & Succi [11]. Using the lattice Boltzmann method one can solve a discrete kinetic equation for the one-particle velocity distribution functions *f*_{i}’s,
2.1where **c**_{i} is the lattice vector, *δ*_{t} is the timestep and *Ω*_{i} is the collision operator. In this paper, we use the simplest form of the latter, i.e. the Bhatnagar–Gross–Krook operator [12,13],
2.2where *τ* is the relaxation time and is a local equilibrium distribution function, which can take the following form:
2.3where , *w*_{i} is the lattice weight, *u* is the hydrodynamic velocity and *c*_{s} is the speed of sound (repeated Greek indices imply summation). It is worth mentioning that the above form of the local equilibrium distribution function can be derived by a systematic method based on a suitable entropy maximum principle [14].

Solving equation (2.1) one can obtain the macroscopic quantities by taking the suitable moments of the distribution functions,

The pressure can be obtained through the equation of state of an ideal gas, i.e. .

For a specific model, the lattice vector and the lattice weights need to be selected. In our analysis, we will use the parameters of a two-dimensional nine-velocity model (D2Q9). The orientation between the lattice links and the simulation domain can be seen in figure 1 (bottom right).

For this model, the lattice links and the corresponding weights are defined as follows:

Using the Chapman–Enskog expansion one can show that the solution of equation (2.1) results in solutions of the incompressible Navier–Stokes equations with some errors. The errors are related to the finite lattice spacing, timestep and Mach number. For the present problem, one can keep these errors between strict limits by choosing suitable resolution and by setting sufficiently low initial velocity (for more details about the validation of the calculations see below).

In this work, equation (2.1) was solved by a simple streaming and collision procedure, and to model no-slip walls we used the bounce-back rules [15], i.e. distribution functions encountering a solid surface were reflected back in the direction they came from.

## 3. Problem description and regional statistics

We consider a square container with rigid no-slip walls and divide the domain into lattices using the resolution *N*×*N*=1024×1024. The initial vorticity field contains 64 equal-size shielded Gaussian vortices,
3.1placed on a regular lattice in the centre of the domain, with some small perturbations in their individual positions (*x*_{k},*y*_{k}) in order to break the symmetry. The corresponding velocity field of a vortex is given by
3.2where *r*_{0} is the initial vortex radius. Half of the vortices have positive and the rest have negative circulation, forming a chessboard-like initial vorticity field. Each vortex has a dimensionless radius 0.04 and the maxima of their vorticity magnitude is *ω*_{0}=1. Quantities used in this paper are normalized by using the initial r.m.s velocity and half-width of the domain. The initial eddy turnover time is normalized by *θ*=*ω*_{0}*t*/*M*, where *M*^{2} is the number of vortices at the beginning of the simulation. Two sets of simulations were performed and each set included six individual runs. The Reynolds number of the two sets based on the initial r.m.s velocity and half-width of the domain were 5000 and 35 000. The calculation started by allowing this flow field to evolve freely. During the evolution, we recorded not only the macroscopic quantities but also the particle distribution functions. Furthermore, area-averaged statistics were obtained from the simulation results by dividing the overall domain into squares with various sizes, as is shown in figure 1. According to figure 1, region 1 extends to 2/*N* distance from the walls and region 2 corresponds to a union of R1 and R2, where the border of R2 is at 10/*N*. Similarly, the union of R1, R2 and R3 defines region 3 and so on. The borders of R3, R4 and R5 are at 50/*N*, 200/*N* and 400/*N*, respectively. The area of these regions, *A*_{region}, will be used to obtain area-averaged statistics in the subsequent discussion. To study the role of a single wall, we also define left region 2 and left region 4, which are rectangles along the left wall of the container (as shown also in figure 1). The finer regional resolution used in the near-wall region is motivated by the fact that walls act as enstrophy sources. Therefore, we expect higher regional enstrophy gradients near the walls than in the outer region.

## 4. Results

### (a) Evolution of global quantities

In order to obtain some information about the performance of the lattice Boltzmann method in the case of the problem in question, first we studied simple vortex interactions using the same method. As a part of this work, we compared our simulation results with pseudospectral calculations [16]. Obtaining excellent agreement, we started the investigation of the problem outlined in the previous section. First, we studied the evolution from a global point of view, i.e. we considered the evolution of the overall kinetic energy and enstrophy and compared our simulation results with those in the literature. Details about this comparison can be found in Házi & Tóth [17]; here, we give just a short overview of the results. At the very beginning of the simulations there are very intense vortex interactions inside the domain owing to the initial dense vortex field. After this period, the vortices start to interact with the walls. At a later stage the activity in the centre of the domain relaxes and the dynamics is driven by the vortices detaching from the walls. These vortices interact with each other and form larger and larger structures or break up as they interact with those detaching from the wall. This can be seen in figure 2, in which a pair of images of the vorticity are shown at *θ*=80 and 160. The same time instants will be used to evaluate regional statistics in the subsequent analysis. White and black correspond to positive and negative vorticity, respectively. Since large amplitude vorticity is generated at the boundaries, we scaled the images in order to make the vortices in the central region clearly visible too. This means that the limits of the image maps correspond to 10 per cent of the maxima and minima of vorticity. From a qualitative point of view, our simulation results are fully consistent with the observations presented by Li & Montgomery [1], Clercx & Nielsen [18] and Clercx *et al.* [2]. That is, walls with the interacting flow field produce strong small vortices in two-dimensional wall-bounded freely decaying turbulence.

Owing to the generation of new vortices, the decay of enstrophy is much slower than in a double-periodic domain. Studying the evolutions of the normalized enstrophy and kinetic energy as an ensemble average of the simulations at *Re*=5000 and 35 000, we found that the decay of enstrophy follows a power law, *Z*(*θ*)∝*θ*^{n}, with exponent *n*=−1 in the case of *Re*=5000 and the exponent decreased as the Reynolds number increased. In our simulations, the decay exponent is higher, *n*=−1.44, at low Reynolds number, *Re*=5000, but at *Re*=35 000 it is *n*=−0.75, which is quite close to the value extrapolated from the simulation results of Clercx & Nielsen [18] for this Reynolds number. Considering the decay of the kinetic energy, a power-law behaviour is less obvious at high Reynolds number. While the decay followed a slope with *n*=−0.5 at *Re*=5000, there was only a very limited range, for which a fit of *n*=−0.15 seemed to be reasonable at *Re*=35000. The decay of the ratio between enstrophy and kinetic energy also exhibited a power law in our simulations and the exponent *n*=−0.89 was in quite good agreement with the value *n*=−0.8 given by Clercx & Nielsen [18] for *Re*=5000. At *Re*=35 000 our exponent decreased to *n*=−0.62, in a similar manner to that shown in Clercx & Nielsen [18], where the value *n*=−0.65 was achieved at the highest Reynolds number, *Re*=20 000.

A grid convergence study revealed that the discrepancies observed between our simulation results and those presented by Clercx & Nielsen [18] are due to our moderate resolution. Indeed, doubling the resolution of the domain at *Re*=5000 we obtained an exponent *n*=−1.14 for the enstrophy decay, while using as high a resolution as *N*×*N*=4096×4096 we have found *n*=−0.98, which is practically the same value as that obtained in the work cited above.

### (b) Vorticity production

Defining the area-averaged regional enstrophy as , where *ω*=∂*u*_{y}/∂*x*−∂*u*_{x}/∂*y* is the vorticity and is the regional enstrophy, it can be seen in figure 3 that small-scale vortices detaching from the wall play a significant role during the evolution. Indeed, the area-averaged regional enstrophy is higher in the near-wall region (region 3), which means that the strongest vortices are created in a narrow band in the near-wall region. As we moved the border of averaging away from the wall we found that the area-averaged regional enstrophy and its fluctuations decay rapidly.

Normalizing the regional enstrophy with the overall enstrophy we can obtain more information about the evolution. In figure 4*a*, which shows the regional enstrophies for a run at *Re*=5000, we can see that the enstrophy is concentrated at the centre of the domain at the beginning of the simulation due to the initial condition, and regions 1–3 are free of enstrophy. After starting the decay, this situation changes almost immediately and the enstrophy rises in regions 1−3 because of the interactions between vortices and walls. The production of enstrophy saturates and as a consequence the regional enstrophies fluctuate around mean values in the time interval *θ*=50…200. This process is quite robust, as demonstrated in figure 4*b*, where the normalized regional mean enstrophies are shown for six individual runs in the different regions.

Taking these mean values, we obtain a nice mean enstrophy profile as a function of position in figure 5*a*. The profile in logarithmic scale is also shown in the insert. This profile is similar to the axial mean velocity profile of turbulent channel flows including the near-wall region, where indeed a boundary layer develops. Up to 0.1, the slopes of the profiles are similar, i.e. the enstrophy increases roughly with a constant rate independently from the initial Reynolds number.

It is worth noting that an increase in the resolution of the simulation domain also has some effect on the regional enstrophy profile, as shown in figure 5*b*. Here, one can see the regional enstrophy profiles at *Re*=5000 using the resolutions *N*=1024,2048 and 4096. Obviously, at higher resolution the regional enstrophy gradient is better predicted in the near-wall region (<0.01).

### (c) Enstrophy statistics

In order to obtain further information about the statistics of vorticity in the near-wall region, we calculated the PDFs of vorticity in left regions 2 and 4, i.e. at the left-hand side of the domain (see figure 1). Figure 6 shows a pair of these PDFs in line with the images shown in figure 2. While the shapes of PDFs for the overall domain remain more or less symmetric during the simulation and peak at the value zero owing to the effect of the walls, the regional PDFs are skewed and their tails change significantly depending on the actual level of vorticity at the left wall. Note that the tails of the PDFs are shorter in the near-wall region than in the overall domain, demonstrating that the production of high-intensity vorticity is more effective at other walls at a given time instant or that vortices with higher level vorticity were formed previously and were advected to another part of the domain. The skewness of the regional PDF is due to the boundary layer, which develops at the wall. The skewness is changing because the boundary layer interacts with interior vortices and breaks up.

## 5. Conclusion

Using the lattice Boltzmann method we simulated two-dimensional decaying turbulence in a square container. We showed that, after the initial stage of decay, the regional area-averaged enstrophy fluctuates around a mean value in time. The fluctuations are strongest in the near-wall regions because of vortices detaching from the wall. The ratio of the regional mean and the overall enstrophies increases monotonously as we move away from the wall and the shape of the mean enstrophy profile is similar to the axial mean velocity profile of turbulent channel flows. We determined the regional PDF of the vorticity. We pointed out that this function peaks at zero and is skewed as we approach the wall because of the boundary layer.

## Acknowledgements

This work was supported by the National Development Agency and the European Union within the framework of project TAMOP 4.2.2-08/1-2008-0021 at the Szechenyi Istvan University, entitled ‘Simulation and optimization—basic research in numerical mathematics’. G.H. was supported by the János Bolyai Research Scholarship of the Hungarian Academy of Sciences.

## Footnotes

One contribution of 25 to a Theme Issue ‘Discrete simulation of fluid dynamics: applications’.

- This journal is © 2011 The Royal Society