Huiyuan Li

State Key Laboratory of Computer Science,Institute of Software Chinese Academy of Sciences,Beijing 100190,China.

Received 16 January 2014;Accepted 14 February 2014

Hexagonal Fourier-Galerkin Methods for the Two-Dimensional Homogeneous Isotropic Decaying Turbulence

Huiyuan Li?

State Key Laboratory of Computer Science,Institute of Software Chinese Academy of Sciences,Beijing 100190,China.

Received 16 January 2014;Accepted 14 February 2014

.In this paper,we propose two hexagonal Fourier-Galerkin methods for the direct numerical simulation of the two-dimensional homogeneous isotropic decaying turbulence.We first establish the lattice Fourier analysis as a mathematical foundation.Then a universal approximation scheme is devised for our hexagonal Fourier-Galerkin methods for Navier-Stokes equations.Numerical experiments mainly concentrate on the decaying properties and the self-similar spectra of the two-dimensional homogeneous turbulence at various initial Reynolds numbers with an initial flow field governed by a Gaussian-distributed energy spectrum.Numerical results demonstrate that both the hexagonal Fourier-Galerkin methods are as efficient as the classic square Fourier-Galerkin method,while provide more effective statistical physical quantities in general.

AMS subject classifications:65M70,65T50,76F05,76F65,76M22

Fourier-Galerkin methods,hexagonal lattices,homogeneous isotropic turbulence,direct numerical simulation.

1 Introduction

The study of two-dimensional homogeneous isotropic decaying turbulence presents several interests,because of not only its applications to geophysics and astrophysics,but also the basic understanding to hydrodynamic turbulence.Far from being a simplified version of the three-dimensional problem,two-dimensional turbulence presents a rich panorama of new phenomena[3].There are many remarkable characteristics in 2D turbulence fields,such as coherent structures[1,13,19],inverse energy cascade and direct enstrophy cascade[2,14,15].The inverse energy cascade indicates thatthe energy is transferred in the inviscid limit from small scales to large scales instead of from large scales to small scales as in the three dimensions.While the enstrophy exhibits a direct cascade process as energy cascade in 3D.The double cascade makes the study of 2D turbulence even more complicated and challenging than that in 3D.

The Navier-Stokes equation for the two-dimensional turbulence may be written in the velocity-vorticity form

where the kinematic viscosityνcan be interpreted as the reciprocal of the Reynolds number.Letψbe the stream function.Then

Thus,the Navier-Stokes equation for the two-dimensional turbulence can also be written as the stream function-vorticity equation,

Navier-Stokes equation in two dimensions has the appealing feature to be less demanding on a computational level than the three-dimensional case,allowing to reach relatively highRenumbers in direct numerical simulation(DNS).As a powerful tool,DNS provides some useful theory check and inspires the deeper thoughts of the statistical theory.The spectral method has been becoming very popular in the research of highly accurate numericalsimulations since the pioneer work ofOrszag and Patterson[21].DNS of 2D turbulence followed an explosive trend in the early stage with an increasing resolution from 2562to 40962or even higher[2,4,5,11,18,23].Up to the present,DNS of 2D homogeneous isotropic decaying turbulence are usually been carried out with a millennial resolution using variants of Fourier spectral or pseudospectral methods with the tensorial Fourier basis functions{ei(k1x1+k2x2)}-n/2≤k1,k2≤n/2-1subject to periodic boundary conditions[5,18,23].

From a general view,the classic Fourier basis functions are just samples of the complex exponentialeiξ·xon the rectangular lattice Z2in the frequency space,i.e,with the wave vectorsξ=k∈Z2.Nevertheless,they actually form a complete orthogonal system on the Voronoi cell{x∈R2:-π≤x1,x2≤π}of the dual lattice 2πZ2in the physical space,which can represent the solution of(1.1)-(1.2)by a Fourier series with its coefficients to be determined.Inspired by the success of the plane-wave method in quantum physics,one can also choose a proper latticeL?(e.g.2πA-1Z2with certain nonsingular matrixA)in the frequency space such that{eiκ·x:κ∈L?}form a complete orthogonal system on the Voronoi cell Ω of the dual latticeL(e.g.,AZ2).Then a Fourier spectral method analogue to the well-known plane wave method on Ω can be established.

Among various admissible lattices,the hexagonal lattice is undoubtedly one of the most favorable choices.It has been known for decades that isotropically band-limited signals are sampled 13.4%more efficiently by hexagonal lattices than by rectangular lattices[22].And it was also shown that the processing algorithms for hexagonal systems are similarly 25-50%more efficient than those for rectangular systems with the same frequency resolution[20].These advantages initiate us to adopt the hexagonal lattice in frequency space and propose the hexagonal Fourier spectral methods which may promote the efficiency of the classic rectangular Fourier spectral methods both theoretically and numerically.

As it can be predicted,the evaluation of the nonlinear term is the main difficulty in the implementation process of our new methods since it dominates the whole time expense.Generally speaking,a computational cost in O(n4)arithmetic operators is required in a direct computation.To conquer this difficulty,we follow the idea in a classic Fourier-Galerkin spectral method.At first,a discrete Fourier analysis and the corresponding hexagonal fast Fourier transform(FFT)are introduced,which actually give an efficient routine for the transformation between the function values in the physical space and the Fourier coefficients corresponding to the frequency space.Then a fastalgorithm in a computational complexity of O(n2logn)is devised for the nonlinear Jacobian term through the convolution formula and our FFTs developed.

In all,the primary goal of this paper is to propose the hexagonal Fourier spectral methods and develop their corresponding fast implementation algorithms with a potential promotion in efficiency.In the sequel,we shall concentrate ourselves to prove the effectiveness and efficiency of our new Fourier spectral methods with a series of numerical experiments appeared in the recent physical study on the DNS of the 2D homogenous isotropic turbulence.We would like to emphasize that all our experiments here are carried out only to recover some interesting results in a classic DNS study by the rectangular Fourier-Galerkin method,without any intention of new mechanism exploration or new phenomena discovery.

The paper is organized as follows.In the next section we present some background materials on planar lattices,Fourier analysis and fast Fourier transform,which set up a mathematical foundation to our further study.Both the dual and the uniform hexagonal Fourier-Galerkin spectral methods are proposed in a unified framework in Section 3 for the DNS of the 2D homogeneous isotropic turbulence.Efficient implementation features,including initial conditions and the fast evaluation of the nonlinear term,are described in details.The performances of our new methods in CPU and GPU environments are also analyzed.Section 4 is then devoted to a numerical comparison between our new methods and the classic Fourier spectral method for the study on the decay of homogeneous isotropic turbulence at both low Reynolds numbers and high Reynolds numbers.The self-similar behaviors of the energy and the enstrophy spectra are examined numerically in Section 5 through our Fourier methods proposed.Finally,conclusions and remarks are drawn in Section 6.

2 Mathematical foundation

2.1 Planar lattices and Fourier series

A planar lattice is a discrete subgroup ofR2which represents a regular,periodic array of points in the plane.One can specify each planar latticeLwith a nonsingular matrixAsuch that

In this case,Ais called a generator matrix of the latticeL.Any plane lattice has a dual lattice(reciprocal lattice)L?given by

wherextdenotes the transpose ofx,andxty=x·yis the usual Euclidean inner product ofxandy.The generatormatrix ofthe duallattice ofAZ2isA?:=2πA-t,i.e,(AZ2)?=A?Z2[16].

A bounded domain Ω?R2is said to(strictly)title the planeR2with the latticeLif

whereχΩdenotes the characteristic function of Ω.We write this as Ω+L=R2,and call such an Ω a primitive cell(fundamental domain)ofL.There is no unique way to choose a primitive cell.The parallelogramAx:x∈[1,0)2is an obvious primitive cell ofAZ2;however,such a primitive cell does not necessarily reveal the underlying symmetry of the lattice in which it is embedded.Fortunately,the Wigner-Seitz cell(Voronoi cell)gives a common choice for a primitive cell with the fully symmetry of the given lattice.A Wigner-Seitz cell about a lattice point is the region which is closer to that point than any other lattice points.Hereafter,whenever Ω is used as a Wigner-Seitz cell of a latticeL,we fix it as the primitive Wigner-Seitz cell ofLabout 0,i.e.,

It is worthy to note that the area of a primitive cell Ω is uniquely determined by the latticeLitself.This invariant,independent of the choice of cell,is denote by|L|.For a given latticeL=AZ2,|Ω|=|det(A)|=|L|.

Figure 1:The square lattice(left),the hexagonal lattice(center)and the dual hexagonal lattice(right).The polygonal regions are the Wigner-Seitz cells about the origin.

Hexagonal lattice.Taking

Tiling and Fourier analysis are closely related(cf.[12])due to the following theorem.

Theorem 2.1(Fuglede[12]).LetΩbe a bounded domain inR2and L be a lattice ofR2.ThenΩ+L=R2if and only if?φκ(x)=eiκtx:κ∈L??is an orthonormal basis,

A pointxis said to be congruent toy∈R2with respect to the latticeLifx∈y+L.Further a functionfdefined onR2is called periodic with respect toLif and only iff(x)=f(y)for anyx∈y+Landy∈R2.If no confusion would arise,we simply callfperiodic.

Let Ω be any primitive cell of the latticeL.Assumefis a periodic function with respect toL,which is square-integrable on Ω.Thenfhas the following Fourier series,

2.2 Discrete Fourier transform

Figure 2:The set of nodes/wavenumbers 2nπ ΛSn/ΛSn in physical/reciprocal space(left),the index set ΛSn=(center)and the congruence index set Λn=Z2∩[0,n)2(right)in the classic DFT.

k∈Z2:-n≤2k1-k2,2k2-k1,k2+k1＜n=k∈Z2:H?k∈nΩ?H.The DFT on dual hexagonal lattices coincides with the classic DFT up to a permutation,thus can be efficiently evaluated through classic FFT just as the DFT on square lattices.

Figure 3:The set of wavenumbers in reciprocal space(left),the index set (middle)and its congruent index set Λn(right side).The congruence relation of the index sets provides a storage scheme for the discrete Fourier coefficients.

Figure 4:The set of quadrature/interpolation nodes Γn(middle)and its congruent index set Λn(right side).The congruence relation of the index sets provides a storage scheme for the function values in physical space.

for which the hexagonal FFT is designed by a generalized Cooley-Tukey algorithm based the periodicity matrix factorization

The hexagonal FFT is decomposed into three classic FFTs of sizen×n,a series of twiddle factor multiplications and a data permutation,which can be fulfilled eventually in O(n2logn)arithmetic operations[9,17,24].

Figure 5:The set of nodes/wavenumbers in physical/reciprocal space (left), the index set and its congruent index set

3 Hexagonal Fourier-Galerkin methods

In this section,we shall propose a lattice Fourier spectral method and its implementation for the two-dimensional homogeneous isotropic turbulence.This lattice Fourier spectral method serves a canonical framework for both the hexagonal Fourier spectral methods and the classic square Fourier spectral method.

3.1 Fourier-Galerkin approximation schemes

The semi-discrete Fourier-Galerkin spectral approximation for the direct numerical simulation of 2-dimensional turbulence is to find

such that

As for time integration,we propose the third-order semi-implicit Adams-Bashforth(AB3)scheme,

Two types of hexagonal Fourier-Galerkin methods are proposed in this paper for the direct numerical simulation of the two-dimensional homogeneous isotropic turbulence.

3.2 Initial conditions

Taking the Fourier transform on(1.2),we have

which gives

Thus,by applying the Parseval’s theorem on the formula of the kinetic energy,one obtains

By using the polar coordinatesξ=(kcosθ,ksinθ)t,one further derives

We specify initial conditions for our flow field by assuming an energy spectrum of wave number magnitudekof the general form

The initial vorticity fieldω0is generated in frequency space with random phases and with amplitude corresponding to(3.3),i.e.,

withζa different uniform deviate for eachκsubject to the requirement of complex conjugate symmetry of the Fourier components.

Given the initial energy spectrumE(k,0),a numerical simulation is uniquely identified by its Reynolds numberR(t)att=0,where the Reynolds number at timetis de fined by

hereafter we use the notations〈·2〉for〈·,·〉.In the case of the initial energy spectrum(3.3),one has

3.3 Spectral coefficients of Jacobian

and get the function valuesw(z)=u(z)v(z),z∈Ym.Further by a forward FFT,we obtain the discrete Fourier coefficients

3.4 Performance analysis

Numerical experiments are performed in a heterogeneous system equipped with GPU devices.Each node has a dual Hexacore Intel X5660 CPU(12M Cache,2.80GHz,6.40 GT/s Intel QPI)and 32GB Memory.Each node is connected with Nvidia Tesla M2070 GPU(1.15GHz,448 cores,6GB 1.5GHz DDR5 Memory).For classic FFTs,the FFTW Library 3.0 is adopted in CPU programs,while the CUFFT Library 4.0 is used in GPU programs.CUFFT has a FFTW compatible data layouts but is optimized only for the input data sizen=2a3b5c7din each direction.

The algorithm described in Section 3.3 withm=3/2n,which is usually referred to as the 3/2-rule for de-aliasing,is adopted for the evaluation of the Fourier coefficients of the nonlinear Jacobian term.To avoid the expensive,time-consuming data transfer between the host(CPU)memory and the devise(GPU)memory in the GPU programs,we would like to keep the entire datum residing on the devise memory all over the time.

At first,the average elapsed time per time step in the CPU programs for Method I,II and III are reported in Table 1,respectively.Roughly speaking,the elapsed time for the Jacobian calculation takes more than 80%of the total elapsed time,so that the total time expense is dominated by the cost of FFTs for the computation of the Jacobian(also refer to Fig.6).Another obvious observation is that Method II becomes more efficient as the resolution increases,and it is about 2 times faster than other methods forN≥40962.This is due to our hexagonal FFT algorithms,which decompose the entire FFT into 3 small classic FFTs,and thus make full use of CPU cache to improve the hit rate of memory access.We also find that time cost for Method I is only slightly higher than that of Method III,since Method I uses the same classic FFT as Method III,and differs from Method III only in a less regulardata structure(referto Figs.2-4),which finally results in a slowdown in the program speed.

Figure 6:The elapsed time for the Jacobian evaluation per time step in comparison to the total elapsed time in three methods with different resolutions on CPU(left),and M2070 GPU(right).

Moreover,we present the average elapsed time per time step for the GPU programs in Table 2.Significant acceleration has been demonstrated in comparison to the corresponding CPU programs.For a high resolutionN=40962,a 45-fold speedup is obtained for Method I,and even higher speedup(49x)for Method III.However,the GPU performance of Method II is relatively low owing to the complicated non-uniform data structure and the corresponding conditional branches.On a conditional branch where the threads diverge in which path to take,the threads taking different paths have to run serially.In spite of the serious performance degradations caused by thread divergence,Method II still attains a speedup about 20 times on the Nvidia Tesla M2070 GPU.

We shallconclude this section with Fig.7,which indicates the Gflops outputdecreasesevidently in our CPU programs-especially in those for Method I and Method III,asNincreases.Fortunately,however,there is a continuing uptrend in the Gflops output in the GPU programs as the resolutionNincreases,which partially demonstrates that GPU implementation of our Fourier-Galerkin methods has a better scalability and thus is adequate for the DNS of 2-D homogeneous isotropic decaying turbulence.

Table 1:The elapsed time(ms)for CPU programs.

Table 2:The elapsed time(ms)for GPU programs.

Figure 7:Gflops statistics for three kinds of Fourier-Galerkin spectral methods versus different scale N.

4 Decay of homogeneous isotropic turbulence

4.1 Decay at low Reynolds numbers

We first consider the evolution of the flow field at relatively low initial Reynolds numbers.Chasnov postulated the existence of a critical initial Reynolds number above which the Reynolds number of the turbulence increases asymptotically,and below which it decreases to small values eventually attaining the final period of decay[8].

Just as in[8],we define the logarithmic derivative in time of the energy and enstrophy as follows,

The advantage of the above definition is obvious-the logarithmic derivatives are just the power-law exponents ifthe energy and enstrophy decay as powerlaws in time.Moreover,we shall use the following normalized time

which can be considered a measure of the number of eddy turnover times undergone by the flow at timet.The normalized timeτbest represents the time interval over which one expects significant changes in the power-law exponent.

Numerical results are first reported for the energy and enstrophy decay of the dual hexagonalFourier-Galerkin method(Method I)and the uniform hexagonalFourier-Galerkin method(Method II)with initial Reynolds numberR(0)=8 and spatial resolutionN=20482.The logarithmic derivativespandq,are plotted versusτin Fig.8,in the left and the middle columns for Method I and Method II,respectively.It is obvious that,at large times,pandqapproach-2 and-3,respectively,and the final period of decay solution[7,8]

is approached asymptotically in two hexagonal Fourier spectral methods.To make a deep comparison,we plotthe logarithmic derivativespandqofthe classic square Fourier-Galerkin method(Method III)in the right column.The plots of two hexagonal Fourier-Galerkin methods agree well with those of the classic Fourier-Galerkin method,so that no significant differences can be observed.

Figure 8:Evolution of the logarithmic derivative of the energy(p)and enstrophy(q)for R(0)=8 with N=20482.Left:Method I-the dual hexagonal Fourier-Galerkin method;middle:Method II-the uniform hexagonal Fourier-Galerkin method;right:Method III-the classic square Fourier-Galerkin method.

Figure 9:Time evolution of the Reynolds number R for R(0)=14,15.73 and 18 with N=20482.Left:Method I;middle:Method II;right:Method III.

Next,we consider some larger initial Reynolds numbers.The numerical experiments are performed corresponding tokc=300,ˉu0=1 andR(0)=14,15.73,18 in spatialresolutionN=20482.Just as in the classic square Fourier-Galerkin method,the numerical results presented in Fig.9 demonstrate the existence of a critical Reynolds numberRc≈15.73 in the two hexagonal Fourier-Galerkin methods such that forR(0)＜Rcthe Reynolds number decays monotonically in time and forR(0)＞Rcthe Reynolds number decreases initially,and then increases asymptotically.

Results for the logarithmic derivatives of the energy and enstrophy withR(0)=Rcversusτare shown in Fig.10 which yield an approximate power-law decay of the energy ast-1,and an approximate enstrophy decay ast-2.An analytical derivation of these power-law exponents in[8]gives

whereR′c≈ 12.5 is a time-independent constant which the Reynolds numberR(t)approaches at large times.The energy and enstrophy decay forR(0)=15.73 is compared to the analytical prediction(4.3).The simulation results of each Fourier-Galerkin method and analytical solution are in good agreement at large times.

Figure 10:Evolution of the logarithmic derivative of the energy(p)and enstrophy(q)for R(0)=14,15.73,18 with N=20482.left:Method I;middle:Method II;right:Method III.

Figure 11:Time evolution of the energy and enstrophy for R(0)=15.73 with N=20482.The straight lines indicate the analytic prediction in(4.3).Left:Method I;middle:Method II;right:Method III.

4.2 Decay at high Reynolds numbers

WhenR(0)＞Rc,the Reynolds numberR(t)increases asymptotically,and some unique similarity states exist since all flows with initial Reynolds numbers greater thanRcpresumably approach infinite Reynolds numbers ast→∞.

Figure 12:Time evolution of the kinetic energy for R(0)=32,64,···,4096 with N=40962.Left:Method I;middle:Method II;right:Method III.

Figure 13:Time evolution of the enstrophy for R(0)=32,64,···,4096 with N=40962.Left:Method I;middle:Method II;right:Method III.

Figure 14:Time evolution of the palinstrophy for R(0)=32,64,···,4096 with N=40962.Left:Method I;middle:Method II;right:Method III.

The logarithmic derivative of the energyp,plotted in Fig.15,is observed to be an increasing function of the initial Reynolds number,which reveals that the energy decay becomes less steep with increasing initial Reynolds numbers.This observation is also rather obvious from the energy decay itself in Fig.12.

Further,different qualitative behaviours ofpare found in[8]for low and high initial Reynolds numbers.ForR(0)≤256,a rapid initial decay of the energy is observed in Fig.12,which subsequently becomes less steep as time evolves due to the increasing Reynolds number of the turbulence.This meanspdecreases to a minimum and then increases in time,just as shown in Fig.15.However,the flows with high initial Reynolds number(R(0)≥512)are nearly inviscid in so much that the energy,Figs.12 and 15,decays very little over the times simulated.Moreover,the magnitude ofpis quite small,and the slow decrease ofpin time implies that the energy decay is steepening as time evolves.

Meanwhile,no universal decay exponent of the enstrophy,Fig.16,at large times for all initial Reynolds numbers greater thanRcis observed.However,as Chasnov discovered,the long-time decay exponent forR(0)＞1024 appears to change only slightly and the asymptotic decay law for these large Reynolds numbers behaves approximately ast-0.8.Besides,since the energy dissipation is proportional to the enstrophy,the change in the qualitative behaviour of the logarithmic derivative of the energy coincides with the power-law exponent of the enstrophy increasing from less than to greater than negative one.

Finally,we emphasise that both plots of the two hexagonal Fourier-Galerkin methods are almost identical to those of the classic square Fourier-Galerkin method,except for the logarithmic derivative of the enstrophy in Fig.16,where slight differences are observed among the three types of spectral methods.

Figure 15:Time evolution of the logarithmic derivative of the energy(p)for R(0)=32,64,···,4096 with N=40962.Left:Method I;middle:Method II;right:Method III.

Figure 16:Time evolution of the logarithmic derivative of the enstrophy q for R(0)=32,64,···,4096 with N=40962.Left:Method I;middle:Method II;right:Method III.

5 Self-similar energy/enstrophy spectra

5.1 Self-similar energy spectra

The energy spectrumE(k,t)oftwo-dimensionalturbulence is found to decay self-similarly,i.e.,without change of shape.

Figure 17:Evolution of the normalized energy spectrum at time τ =0,5,10,···,35 with R(0)=15.73 with N=40962.

Figure 18:Rescaling of the energy spectrum of Fig.17 with N=40962.

Figure 19:The energy spectrum at time t=0,10,20,···,120 for R(0)=64 with N=40962.

Figure 20:Rescaling of the energy spectra in Fig.19 with N=40962.

which indicates substantialnonlinear backscatterofenergy from small-to-large scale even for this low Reynolds number turbulence.

We further examine the time evolution of the energy spectra for flows with large initial Reynolds numberR(0)＞Rc.The time evolution of the spectra forR(0)=64,256,and 4096 are shown in Fig.19,Fig.21 and Fig.23,respectively.For a better validation of the self-similar energy spectrum,the self-similar spectra obtained using(5.1)are plotted in Figs.20,22 and 24,immediately below the corresponding figure for the spectra evolution.It appears from the reasonable collapse of the spectra at different times that the decay of two-dimensional turbulence at large Reynolds number is also self-similar in the energy containing scales.

In Figs.20,22 and 24,we have plotted as a dashed line the expectedk3low wave number behavior for two-dimensional turbulence.Moreover,Batchelor pointed that the energy spectrum in fully developed turbulence will acquire the self-similarE(k)∝k-3over the so-called “inertial range”,i.e.,wavenumberskextending from the “energycontaining”scale wavenumberk0to the viscous scale wavenumberkν～ν-1/2.An inertial subrange appears to have developed in Figs.20 and 22 forR(0)=256 andR(0)=2048,which is slightly steeper than the predictedk-3behavior.In particular,when the initial Reynolds number increases toR(0)=4096,this inertial subrange(in Fig.24)developed in the two hexagonal Fourier-Galerkin methods perfectly matches the predictedk-3inertial subrange behaviour for two-dimensional turbulence.

Figure 21:The energy spectrum at time t=0,10,20,···,120 for R(0)=256 with N=40962.

Figure 22:Rescaling of the energy spectra in Fig.21 with N=40962.

Figure 23:The energy spectrum at time t=0,10,20,···,120 for R(0)=4096 with N=40962.

Figure 24:Rescaling of the energy spectra in Fig.23 with N=40962.

5.2 Self-similar enstrophy spectra

Batchelor[2]predicted thatthe enstrophy spectrumin fully developed turbulence willacquire the self-similar form Ω(k)∝χ2/3k-1(withχ=ν〈|?ω|2〉)over a range of wavenumberskextending from the“energy-containing”scales,say around wavenumberk0,to the viscous scales,say aroundkν～ν-1/2,i.e.,over the rangek0?k?kν,the so-called “inertial range”.Batchelor’s theory of two-dimensional turbulence is concluded by assuming that there is a finite,non-zero enstrophy dissipationχin the limit of infiniteR.This theory has been successful in describing certain aspects of numerical simulations at high Reynolds numbers.However,Dritschelet al.[10]found that the enstrophy dissipation in fact vanishes for flows with finite vorticity,and the non-zero enstrophy dissipation assumption is not true.After a careful observation together with a mathematical analysis of vanishingχ,Dritschel finally made a simple modification of Batchelor’s theory by replacing the enstrophy spectrumχ2/3k-1with 〈ω2〉k-1(lnR)-1.

which peaks atkc≈3.2 and is more than 1036times smaller byk=32.Such a spectrum is already peaked at low wavenumbers,no significant inverse energy cascade takes place over the short duration of our simulations.

The viscosity is chosen so that the approximate viscous wavenumber occurs at 3/4 of the maximum effective wavenumbern/2 in the classic square Fourier-Galerkin method,i.e.,ν=4π/(3n/8)2.This choice was found to ensure a finite initialvorticity with kω(0)k∞≈4πand adequate dissipation at high wavenumbers to resolve the statistics of〈(Δω)2〉and other fine-scale quantities.

Numerical experiments have been performed with the spatial resolutionN=40962,and the initial Reynolds numberR(0)=5.55×104(in view of(3.5))using the two hexagonal Fourier-Galerkin methods.These resolutions are high enough to examine the nature of dissipation in two-dimensional turbulence.Numerical results are,once again,compared with those of the classic Fourier-Galerkin method.

Since the behaviour of the enstrophy dissipation depends on the palinstrophy evolution,we first concentrate on the palinstrophy spectra,which are shown in Fig.25.Increasingly,a range close to the classicalk1spectrum develops at low wavenumbers.The spectrum has filled out completely by the palinstrophy peak,and afterward the spectrum at moderate to high wavenumbers decays exponentially while preserving its basic shape.It is meaningful to say that the turbulence is“fully developed” by the time the palinstrophy reaches its maximum.Fig.25 exhibits no obvious distinction between the palinstrophy evolutions of the three cases.While the spectrum at low wavenumbers is much closer to the referencek1slop in the plot for Method II.

Figure 25:Palinstrophy spectra k2Ω(k,t)for N=40962 at t=0 and at times when the palinstrophy is half its peak value(and growing),at its peak,and half its peak(and decaying).A prediction of k1 is shown as a dashed line.Left:Method I;middle:Method II;right:Method III.

Figure 26:The scaled enstrophy spectra Ω(k)ln R/〈ω2〉at the peak enstrophy dissipation time for N=10242(grey thin),N=20482(thin)and N=40962(thick).A reference slope of k-1 is shown as a dashed line.

Figure 27:The scaled enstrophy spectra Ω(k)ln R/〈ω2〉at the final time for N=10242(grey thin),N=20482(thin)and N=40962(thick).A reference slope of k-1 is shown as a dashed line.

Further the scaled enstrophy spectra Ω(k)lnR/〈ω2〉at the peak enstrophy dissipation time and the final time forN=10242,20482and 40962have been shown in Fig.26 and Fig.27,respectively.The curves collapse together over the common inertial ranges.An interpretation of this phenomena is that dissipation is playing an increasingly benign role asR→∞.And the 〈ω2〉(klnR)-1enstrophy spectrum reflects that a fixed amount of enstrophy must spread itself ever more thinly across a widening inertial range as the Reynolds number grows.Once again,there is only a slightdifference between the simulation results of the three kinds of Fourier-Galerkin spectral methods.Meanwhile,Method IIprovides us a slop ofthe enstrophy spectra over the“inertialrange”closer tok-1,which demonstrates the modified self-similarity theory of Dritschel.

6 Conclusion

We have proposed two hexagonal Fourier-Galerkin spectral methods for the direct numerical simulation of the two-dimensional homogeneous isotropy freely-decaying turbulence to obtain an even better sampling/approximation efficiency to the band limited/isotropic functions in comparison to the classic rectangular Fourier spectral method.As the soul of the implementation of the approximation scheme,an efficient algorithm for evaluating the de-aliased Fourier coefficients of the nonlinear term is well designed.By adopting some variants of FFTs,our new Fourier spectral methods preserve the high efficiency of the classic Fourier spectral method;and a higher performance can even be attained in our methods.Numerical experiments show that our hexagonal Fourier-Galerkin spectral methods acquire a slightly better statistical results in general than those of the classic rectangular Fourier spectral method,which are in agreement with the corresponding physical theories.Anyway,hexagonal Fourier-Galerkin spectral methods provide alternative choices for the DNS of the two-dimensional isotropy turbulence indeed.

Acknowledgments

This work was partially supported by Natural Science Foundation of China(No.91130014),and also supported in part by the Guangdong Province Key Laboratory of Computational Science and the Guangdong Province Computational Science Innovative Research Team.

We would like to acknowledge Professor Jie Shen for his valuable suggestion and comment on the current research topic.And we thank Mr.Haijun Qiao to prepare the figures in this paper.

[1]H.Aref and E.D.Siggia.Vortex dynamics of the two-dimensional turbulent shear layer.J.Fluid Mech.,100:705-737,1980.

[2]G.K.Batchelor.Computation of the energy spectrum in homogeneous two-dimensional turbulence.Physics of Fluids,12(12):II-233,1969.

[3]G.Boffetta and R.E.Ecke.Two-dimensional turbulence.Annual Review of Fluid Mechanics,44(1):427-451,2012.

[4]V.Borue.Spectral exponents of enstrophy cascade in stationary two-dimensional homogeneous turbulence.Physical Review Letters,71(24):3967,1993.

[5]A.Bracco,J.C.McWilliams,G.Murante,A.Provenzale,and J.B.Weiss.Revisiting freely decaying two-dimensional turbulence at millennial resolution.Physics of Fluids,12(11):2931-2941,2000.

[6]C.Canuto,M.Y.Hussaini,A.Quarteroni,and T.A.Zang.Spectral Methods in Fluid Dynamics.Springer-Verlag,1988.

[7]J.R.Chasnov.Similarity states of passive scalar transport in isotropic turbulence.Physics of Fluids,6(2):1036-1051,1994.

[8]J.R.Chasnov.On the decay of twodimensional homogeneous turbulence.Physics of Fluids,9(1):171-180,1997.

[9]J.Chen,H.Li,and X.Zhang.A CUDA-MPI algorithm for the fast Fourier transform on the hexagon and its implementation(in Chinese).Journal on Numerical Methods and Computer Applications,33(1):59-72,2012.

[10]D.G.Dritschel,C.V.Tran,and R.K.Scott.Revisiting Batchelor’s theory of two-dimensional turbulence.J.Fluid Mech.,591:379-391,2007.

[11]U.Frisch and P.-L.Sulem.Numerical simulation of the inverse cascade in two-dimensional turbulence.Physics of Fluids,27:1921,1984.

[12]B.Fuglede.Commuting self-adjointpartialdifferentialoperators and a group theoretic problem.Journal of Functional Analysis,16(1):101-121,1974.

[13]A.E.Hansen and P.Tabeling.Coherent structures in two-dimensional decaying turbulence.Nonlinearity,13(1):C1-C3,2000.

[14]R.H.Kraichnan.Inertial ranges in two-dimensional turbulence.Physics of Fluids,10:1417,1967.

[15]M.Lesieur.Turbulence in Fluids,volume 84.Springer,2008.

[16]H.Li,J.Sun,and Y.Xu.Discrete fourier analysis,cubature,and interpolation on a hexagon and a triangle.SIAM Journal on Numerical Analysis,46(4):1653-1681,2008.

[17]H.Li and J.Sun.Fast Fourier transform on hexagons.In W.Zhang,W.Tong,Z.Chen,and R.Glowinski(Eds.),Current Trends in High Performance Computing and Its Applications,pp.357-362.Springer Berlin Heidelberg,2005.

[18]M.S.I.Mallik,M.A.Uddin,and M.A.Rahman.Direct numerical simulation in two dimensional homogeneous isotropic turbulence using spectral method.Journal of Scientific Research,5(3),2013.

[19]J.C.McWilliams.The emergence of isolated coherent vortices in turbulent flow.Journal of Fluid Mechanics,146:21-43,1984.

[20]R.M.Mersereau.The processing of hexagonally sampled two-dimensional signals.Proceedings of the IEEE,67(6):930-949,1979.

[21]S.A.Orszag and G.S.Patterson.Numerical simulation of three-dimensional homogeneous isotropic turbulence.Physical Review Letters,28:76-79,1972.

[22]D.P.Petersen and D.Middleton.Sampling and reconstruction of wave-number-limited functions inN-dimensional euclidean spaces.Information and Control,5(4):279-323,1962.

[23]O.San and A.E.Staples.High-order methods for decaying two-dimensional homogeneous isotropic turbulence.Computers&Fluids,63(0):105-127,2012.

[24]J.Sun and J.Yao.A fast algorithm of discrete generalized Fourier transforms on hexagon domains(in Chinese).Mathematica Numerica Sinica,26(3):351-366,2004.

?Corresponding author.Email address:huiyuan@iscas.ac.cn(H.Li)