Asymptotic states of decaying turbulence in two-dimensional incompressible ﬂows

We investigate the relaxation of a strongly turbulent ﬂuid to metastable states, in times much shorter than the dissipation time scale. Several simulations of decaying two-dimensional Navier-Stokes ﬂows were performed, which show the relaxation to organized states dominated by coherent vortex structures of length scales com-parable to the size of the system. In the case of periodic boundary conditions, the organized state is characterized by a strong correlation between vorticity and stream function, the second of which satisﬁes a sinh-Poisson equation quite accurately. However, in the case of free-slip boundary conditions the relaxed state does not display any signiﬁcant correlation between its vorticity and its stream function. Notwithstanding, in both cases the role of nonlinearities is found to be essential even at these late stages of the evolution. However, we show that even severe truncations of a few (cid:126) short wave number (cid:33) nonlinearly coupled Fourier modes provide an accurate description of the long-term dynamics of the ﬂuid. Therefore the dynamics of the ﬂow in these metastable states is somewhere in between a strong turbulent regime and a (cid:126) mostly linear (cid:33) dissipative relaxation stage. (cid:64) S1063-651X (cid:126) 96 (cid:33) 03506-4 (cid:35)


I. INTRODUCTION
Thanks to the increasing computational resources available in recent years, it has become possible to study longstanding problems like the decay of two-dimensional turbulent Navier-Stokes flows at very long times. Numerical simulations of this process show the formation of coherent vortex structures ͓1͔ at times much shorter than the viscous dissipation time scale. This metastable state relaxes slowly toward the final trivial state, with zero fluid velocity everywhere.
In the coherent, nontrivial state, the system is dominated by the largest scales, with concentration of vorticity in a few cores. The simulations performed in ͓2͔, which correspond to periodic boundary conditions, show a state with two cores that can be adequately described by the sinh-Poisson equation. This behavior has been related to the negative temperature state of the discrete inviscid line vortex model ͓3͔ by a two-fluid formulation of the continuous Navier-Stokes flow ͓4͔.
The temporal decay of two-dimensional turbulence has been studied with high spatial resolution by Brachet et al. ͓5͔,Mc Williams ͓6͔,and Montgomery et al. ͓4͔. The simulations in ͓5͔ display the early formation of macrovortices together with small-scale structures, such as vorticitygradient sheets, associated with an initial inverse cascade process ͑to large scales͒ for the energy and a direct cascade ͑to small scales͒ for the enstrophy. The direction of these cascades is consistent with those obtained for an externally driven two-dimensional turbulent flow. However, these simulations do not extend to times long enough to study the relaxation of the system to metastable states.
Mc Williams ͓6͔ claims that, when vorticity concentrations are sufficiently intense, the cascade processes are virtually suppressed: there is a cessation of spectral transfer and nonlinearities can be neglected. However, as reported by Montgomery et al. ͓7͔, the long-time behavior of the flow, characterized by the formation of two macrovortices, still displays features characteristic of a turbulent flow, like a broadband spectrum and cascading processes.
Nonetheless, even though the energy spectrum at these long times follows a power law that varies as k Ϫ7 or steeper, it is sharply peaked at small wave numbers. It is therefore arguable whether these steep spectra can be considered broadband, as expected for strong turbulent regimes. This suggests that we can describe the long-time behavior of the system with relatively few Fourier modes, taking, however, their nonlinear interactions into account.
Our simulations show that a Fourier truncation to a few long-wavelength modes describes very well the dynamics of the system in this coherent nontrivial stage of the turbulence decay. In the case of periodic boundary conditions, the sinh-Poisson equation is also verified for the truncated system, and the parameters governing the equation evolve in time in very good agreement with the full, nontruncated system. The organization of the paper is as follows. In Sec. II we present the dynamic equations and describe the characteristics of the numerical code. The properties of the metastable states described by the sinh-Poisson equation are summarized in Sec. III. The results for periodic boundary conditions and free-slip boundary conditions are shown in Sec. IV. Section V is dedicated to the study of long-wavelength truncations. The conclusions are listed in Sec. VI.

II. MODEL EQUATIONS AND NUMERICAL SIMULATIONS
Given a two-dimensional incompressible flow in the (x,y) plane, the evolution equation for the scalar vorticity is which is obtained by taking the curl of the Navier-Stokes equation. In Eq. ͑1͒, is the kinematic viscosity, (x,y) is the stream function, which is related to the two-dimensional velocity field by v x ϭ‫ץ‬ y , v y ϭϪ‫ץ‬ x , and w(x,y)ϭ‫ץ‬ x v y Ϫ‫ץ‬ y v x is the scalar vorticity field, which satisfies The Poisson bracket in Eq. ͑1͒ is defined as We consider the flow in a square box of sides Lϭ2 and assume periodic boundary conditions in both directions. Due to the periodicity, the fields can be expanded in a Fourier series, where k is a two-dimensional vector of integers and the asterisk denotes complex conjugation. We numerically solve ͑1͒ and ͑2͒ with a dealiased pseudospectral method, i.e., working in Fourier space, with implicit time integration in the linear term ͑the dissipative term͒ and an explicit, twostep Adams-Bashforth ͓8͔ technique in the nonlinear term. The initial conditions were chosen so as to simulate a highly turbulent state. The modal energy for the initial state is where C and k 0 are constants. Thus the energy spectrum ͓2kE(k)͔ is initially broad, varying as k Ϫ3 at small wavelengths, which in a stationary turbulent regime corresponds to the energy spectrum in the enstrophy inertial range. The modulus of the initial stream function in Fourier space is obtained from The complex phase is chosen randomly with a uniform distribution law in ͓0,2͔. The Reynolds number for this flow is where ͗v 2 ͘ 1/2 is the rms velocity and l is a typical length, both depending on initial conditions. Working with dimensionless units, and choosing a mean velocity and a typical length of order 1, the Reynolds number is of order 1/, which is also the time by which the system decays to the trivial (vϭ0) state. We have performed simulations for Reynolds numbers of the order of 1000.

III. FORMATION OF COHERENT VORTEX STRUCTURES: SINH-POISSON EQUATION
For freely decaying two-dimensional turbulence, energy and enstrophy are inviscid invariants, i.e., conserved quantities if no dissipation is present. Indeed, they satisfy the equations ͓9͔ where ⍀ is the enstrophy (⍀ϭ ͚k 4 ͉ k ͉ 2 ) and P the palinstrophy ( Pϭ͚k 6 ͉ k ͉ 2 ). Also, it can be demonstrated ͓10͔ that the ratio ⍀/E decays monotonically in time. This ratio can be interpreted as a mean square wave number of the flow, since expanding ⍀ and E in their Fourier components, it becomes Thus the monotonic decrease of this ratio can be interpreted as an increase of the average wavelength. In other words, the system evolves in time toward states dominated by the largest scales available ͑the largest allowed by the boundary con-ditions͒.
A first atempt at understanding these organized states has been the mechanism of selective decay ͓11͔. Because enstrophy decays faster than energy, a variational principle is invoked, minimizing enstrophy with the constraint of constant energy. Writing enstrophy as ⍀ϭ 1 2 ͐w 2 d 2 x and energy as Eϭ 1 2 ͐wd 2 x, the variational principle ␦⍀Ϫ␦Eϭ0, with a Lagrange multiplier, leads to a linear equation, Then vorticity is proportional to the stream function in these states. The solutions of Eq. ͑12͒ are single Fourier modes, and the one which ensures minimum ⍀/E is a superposition of modes with ͉k͉ϭ1. However, Montgomery et al. ͓2͔ showed in their simulations that vorticity and stream function in the coherent states can be described by the relationship which is known as the sinh-Poisson equation. The parameters c and ␤ are obtained by a fitting procedure. This relationship can be explained in terms of an alternative variational principle ͑see, for instance, ͓12͔͒. Let us define a positive vorticity field w ϩ and a negative vorticity field w Ϫ in such a way that w Ϯ у0 and wϭw ϩ Ϫw Ϫ . Let us assume that these fields are respectively proportional to the number of positive and negative vortices n Ϯ , each of which satisfies passive scalar evolution equations. Therefore for this gas of vortices we can define the entropy as Let us further assume that the total positive and negative vortices are approximate invariants at large Reynolds numbers (W Ϯ ϭ͐w Ϯ d 2 xϭ const, and W ϩ ϭW Ϫ ), as well as the energy of the system [Eϭ 1 2 ٌ͉͉͐ 2 d 2 xϭ 1 2 ͐(w ϩ Ϫw Ϫ )d 2 x͔. Maximization of S subject to these constraints leads to Since wϭϪٌ 2 , the change →Ϫ implies that w→Ϫw and therefore c ϩ ϭc Ϫ ϭc/2. Thus wϭcsinh(␤). The above variational principle is a simplified version of the variational principle proposed by Mongtomery et al. ͓4͔. These authors define positive and negative vorticities in the following way: at tϭ0 they are the positive and negative parts of the total vorticity; for tϾ0 they are the solutions of evolution equations analogous to Eq. ͑1͒. These equations assure exact conservation of total positive and negative vorticity, for both the inviscid and the viscous cases. Moreover, these invariants are also conserved in any arbitrary truncation in Fourier space.
Other studies of equilibrium statistical mechanics applied to two-dimensional flows can also be found in ͓12-14͔. For instance, Robert and Sommeria ͓14͔ seek for maximal entropy states with the constraints of all the constants of motion of the Euler equations. Their approach provides a rela-tionship between vorticity and stream function, wϭ f (), which characterizes the steady state. This function f displays an exponential behavior and is such that its second derivative has the same sign as , therefore behaving very much like a sinh function.

A. Periodic boundary conditions
In this section we display the results of numerical simulations of Eqs. ͑1͒ and ͑2͒ for periodic boundary conditions in a square box of sides 2 and using a spatial resolution of 96ϫ96 gridpoints. This moderate resolution allowed us to perform several simulations, for different types of initial conditions.
In the early stages of the evolution ͑between tϭ0 and tϭ10) the decay is indeed governed by an inverse energy cascade associated with the formation of macrovortices, and a direct enstrophy cascade connected to the generation of vorticity-gradient sheets, as reported by Brachet et al. ͓5͔ and Mc Williams ͓6͔. At longer times, as reported by Montgomery et al. ͓1͔, the flow is increasingly dominated by the larger structures. The vortices of the same sign merge until just two vorticity cores of opposite sign are left. This behavior is clearly observed in Fig. 1, which displays the stream function contour lines at different times, and in Fig. 2, which shows the velocity vector field over a halftone of the vorticity field.
As mentioned above, in this case the coherent state is well described by the sinh-Poisson equation ͑13͒. The correlation between vorticity and stream function can be seen in Fig. 3, which shows a scatter plot of the vorticity versus stream function for different times. By tϭ60 we can perform a fitting of the two parameters of Eq. ͑13͒ and their evolution in time. We also calculated the correlation between w and sinh(␤), which becomes better than 99% by tӍ200.
Due to the periodic boundary conditions, we can tile the plane with copies of the 2ϫ2 square domain on which the simulations have been carried out ͑see Fig. 1͒. If we do this, we observe that the characteristic length scale of the vortices is such that four vortices, two with positive and two with negative vorticities, fit into a tilted square domain. The sides of this square are rotated by 45°with respect to the original axes and thus have a length equal to 2ͱ2. This arrangement is consistent with two properties that the asymptotic state must satisfy: on one hand, the conservation of total vorticity, which is exactly equal to zero inside each square; on the other hand, the tendency of the system to reach a state of maximum wavelength, compatible with the size of the domain. It is clear that the diagonal is the largest length that can fit into a square, and this is probably the reason why the vortices organize themselves along the diagonals.
The vorticity and the stream function are very small at any point on the sides of the tilted squares. Had we obtained a perfectly symmetrical state, one with four vortices of identical shapes, the vorticity and the stream function would be exactly equal to zero at these sides. For this reason we decided to perform numerical simulations on a square domain with a combination of periodic and free-slip boundary conditions on the stream function, i.e., ϭ0 at the borders.

B. Free-slip boundary conditions
The combination of periodic and free-slip boundary conditions ͑also called Dirichlet boundary conditions͒ is achieved by imposing the condition ϭ0 at the borders in the expansion ͑4͒, which implies that where i, j run over the positive integers. It can be shown that the subspace spanned by the functions sin(ix)sin(jy) is invariant under the evolution given by Eqs. ͑1͒ and ͑2͒.
We have performed a second set of numerical simulations of Eq. ͑1͒ with the same choice of initial conditions as in the previous section but only allowing the values and 0 for the initial phases. We had two motivations for carrying out these simulations. On one hand, we expected to obtain an evolution similar to the one observed in the tilted square of the previous section. Furthermore, we expected it to be easier to analyze in terms of analytical solutions of the sinh-Poisson equation, which are known for Dirichlet boundary conditions.
On the other hand, due to the symmetry that any function of the form ͑16͒ has on the (0,2)ϫ(0,2) domain, solving the Navier-Stokes equation on this domain amounts to solving it simultaneously in four squares of length with Dirichlet or free-slip boundary conditions ͑see, e.g., Fig. 4͒. Free-slip boundary conditions are physically relevant for inviscid fluids. Since we are interested in the evolution for times much shorter than the dissipation time scale in flows with high Reynolds numbers and we are not looking at the dynamics in the boundary layer, we think they are also relevant to our analysis.
The results of the simulations are shown in Fig. 4. The four independent boxes of sides Lϭ can clearly be seen. The mirror antisymmetry of the stream function, imposed by ͑16͒, is also apparent in this figure. The flow in this case is also increasingly dominated by the largest structures allowed by the boundary conditions, which at long times consist of four vorticity cores.
This evolution from an initially disordered state to a rather organized state is shown in the scatter plots of vorticity versus stream function of Fig. 5. This figure convincingly shows that the ordered state cannot be described by the sinh-Poisson equation, or by any other functional dependence between w and .
Given this particular behavior, which is clearly different from the one achieved in the previous section, it seems natural to ask the following two questions. First, why does this happen? Second, since these simulations are a subcase of those with periodic boundary conditions, how typical or how peculiar are these departures from the relaxation to a sinh-Poisson state? As mentioned above, the subspace spanned by the functions sin(ix)sin(jy) is an invariant subspace of the space of periodic functions. However, we have found numerically that is not attracting. Indeed, we have simulated Eqs. ͑1͒ and ͑2͒, considering an initial condition in this invariant subspace plus a small addition of white noise. The solution remained close to the invariant space for a while but finally moved away from it. Furthermore, the simulation showed that the system relaxed to a sinh-Poisson state. There are lots of other invariant subspaces in the space of periodic functions, such as those spanned by functions of the form sin(inx)sin(jny) with n an arbitrary integer, or the tilted square described above. However, they seem to be all unstable under a small perturbation. Therefore, if we consider arbitrary initial conditions in the space of periodic functions ͑with initially broad energy spectrum and random phases͒, it is highly plausible that we will end up relaxing to a final state described by the sinh-Poisson equation.
Since free-slip boundary conditions are physically relevant, it is important to understand why in this case we did not get a relaxation to a final sinh-Poisson state. As mentioned above, the relaxation to a sinh-Poisson state can be explained through a variational principle: entropy S ͓see the definition in Eq. ͑14͔͒ is maximized subject to various constraints. One of these constraints is the conservation of energy, which is only approximately conserved in the dissipative case. We find that in the case with free-slip boundary conditions, energy decays much faster than in the case with periodic boundary conditions. Without the constraint of conservation of energy, the variational principle leads only to a trivial state (wϭ0). Therefore we conclude from our simu- lations that in the case of free-slip boundary conditions the system does not reach a metastable nontrivial state described by the sinh-Poisson equation.
There is an interesting fact, though, which indicates that the relaxation to the sinh-Poisson state plays a role also in this case. At tϭ60 we observe in Fig. 4 that there are eight vortices of two different shapes, four in the middle region and four at the bottom and top of the square. On the other hand, Fig. 5 looks like the plot of two functions, as if there were two possible functional relations between w and . Moreover, one of them looks like a sinh-type curve. In fact, we have found that points coming from the vortices at the top and bottom of Fig. 4 fall into this curve, which can be fitted by one of the form c sinh(␤), at least for values of away from zero. Therefore some of the vortices that are formed in the case of free-slip boundary conditions do relax to a sinh-Poisson state. However, the dissipation enters into play much sooner than in the case with periodic boundary conditions and this effect seems to destroy the correlation between w and .

V. LONG-WAVELENGTH TRUNCATIONS
We return to the case of periodic boundary conditions to study a long-wavelength truncation of the system in the coherent stage of the evolution. Figure 6 shows the energy spectra at three different times. We can see that, however broad, by tϭ60 it is sharply peaked, with a decaying exponential law of k Ϫ7 at short wavelengths. This suggests that we perform a truncation of Fourier modes ͑besides the truncation imposed by the finite resolution͒, and test how it affects the dynamics of the sys- tem at this late stage of the evolution.
In order to do this comparison, we evolve the whole system up to tϭ100. At that time we truncate the system to Fourier modes with ͉k͉р4 and evolve only those modes for tϾ100. The spectra obtained for this drastically truncated system are compared in Fig. 7 to the spectra of the complete system at times tϭ150 and tϭ250. The agreement is very good. Figure 8 shows the stream function of the truncated system at different times, as well as a comparison with the stream function for the complete system. The only noticeable difference between the two evolving systems is just a slight rigid shift of the patterns of the truncated system with respect to the complete system.
The truncated system also evolves to a sinh-Poisson state as shown in Fig. 9, where we have plotted w vs for both the truncated and the complete systems at different times. Moreover, the parameters c and ␤ coincide within the fitting procedure error.
The energy transfer between modes and indirectly the degree of approximation made by the truncation can be quantitatively measured by means of the following procedure. Suppose we consider a truncated system with modes up to a certain kϭk 0 . It follows from the evolution equation ͑Navier-Stokes͒ written in Fourier space that the energy balance equation for this reduced system is where E k 0 ϭ 1 2 ͚ kрk 0 k 2 ͉ k ͉ 2 is the energy of the truncated system, ⌸ k 0 ϭRe͓ ͚ kрk 0 k 2 k *͓,w͔ k ͔ represents the nonlinear energy transfer from modes with kϾk 0 ͑the modes out of the truncated system͒, and 2⍀ k 0 is the dissipation term, with ⍀ k 0 ϭ 1 2 ͚ kрk 0 k 4 ͉ k ͉ 2 the enstrophy of the truncated system. FIG. 6. Energy spectrum for the whole system at three different times.
FIG. 7. Energy spectrum for the truncated system ͑dashed line͒, with maximum wave number of order 4 and for the complete system ͑solid line͒, with maximum wave number of order 40, at times tϭ150 and tϭ250. The truncation starts at tϭ100.
FIG. 8. Stream function contours for the truncated system and for the complete system, at different times, for the case with periodic boundary conditions. The square boxes have sides 2.
FIG. 9. Evolution of the sinh-Poisson state in the truncated system compared to the complete system, for periodic boundary conditions.
The ratio ⌸ k 0 /2⍀ k 0 is a measure of the importance of the nonlinear exchange of energy between modes of kϾk 0 and modes with kрk 0 compared to the dissipative term.
A similar equation can be obtained for the enstrophy balance, where ⌫ k 0 ϭRe͚ kрk 0 k 4 k *͓,w͔ k represents the nonlinear enstrophy transfer from modes with kϾk 0 , and 2 P k 0 is the dissipation term, with P k 0 ϭ 1 2 ͚ kрk 0 k 6 ͉ k ͉ 2 the palinstrophy of the truncated system. The importance of the nonlinear exchange of enstrophy between modes of kϾk 0 and modes with kрk 0 compared to the dissipative term can be measured by the ratio ⌫ k 0 /2P k 0 .
The energy transfer and enstrophy transfer ratios can be seen in Fig. 10 and Fig. 11. In Fig. 10 the ratios are plotted as functions of time for k 0 ϭ4. By tϭ100 both ratios are smaller than 0.06; therefore k 0 ϭ4 is a reasonable choice for truncating the system. In Fig. 11 the ratios are plotted as functions of k 0 for two different times.

VI. CONCLUSIONS
We have studied the relaxation of strongly turbulent twodimensional incompressible flows via numerical simulations of the Navier-Stokes equations. We have performed these simulations for periodic boundary conditions and for freeslip boundary conditions. We have observed a relaxation to a sinh-Poisson state in the first case. In the second case, however, the system did not reach any metastable state in which vorticity could be approximated by any function of the stream function. We think this is due to the fact that energy is dissipated much faster in this second case, thus preventing the system from reaching the above mentioned metastable state. Although the case of Dirichlet boundary conditions considered in this paper is a subcase of the more general case of periodic boundary conditions, we have found that the invariant subspace of solutions which satisfy our free-slip boundary conditions is not stable. Furthermore, an initial condition arbitrarily close to this invariant space finally relaxes to a sinh-Poisson state. For this reason we believe that most initial conditions with arbitrary periodic boundary conditions will reach the metastable state in which the stream function satisfies the sinh-Poisson equation.
For the general periodic case, we performed a truncation to a few long-wavelength modes, neglecting the interaction with the remaining modes ͑modes with kϾk 0 ). The nonlinear interactions among these long-wavelength modes were consistently taken into account. The truncated system approximates very well the dynamics of the complete system in the coherent stage of the evolution, as can be seen by comparing energy spectra and stream function. Moreover, the evolution of the sinh-Poisson state is very well reproduced, as can be seen by comparing the parameters c and ␤ at different times.
We calculated the energy and enstrophy nonlinear transfer rates between the truncated system and the discarded modes and found that, after a certain time much smaller than the dissipation time, and for a certain ͑small͒ k 0 , they are negligible compared to the dissipation terms. This result indicates that the state of the flow is intermediate between a fully developed turbulent regime and a linear dissipative stage.