The Navier–Stokes-αβ equations as a platform for a spectral multigrid method to solve the Navier–Stokes equations

Computers & Fluids 44 (2011) 102–110

Contents lists available at ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

The Navier–Stokes-ab equations as a platform for a spectral multigrid method to solve the Navier–Stokes equations Tae-Yeon Kim a, John E. Dolbow b, Eliot Fried a,⇑ a b

Department of Mechanical Engineering, McGill University, Montréal, Québec, Canada H3A 2K6 Department of Civil and Environmental Engineering, Duke University, Durham, NC 27708-0287, USA

a r t i c l e

i n f o

Article history: Received 20 July 2010 Received in revised form 6 December 2010 Accepted 15 December 2010 Available online 21 December 2010 Keywords: Navier–Stokes-ab equations Spectral multigrid method Strained spiral vortex model

a b s t r a c t This paper describes a spectral multigrid method for spatially periodic homogeneous and isotropic turbulent flows. The method uses the Navier–Stokes-ab equations to accelerate convergence toward solutions of the Navier–Stokes equations. The Navier–Stokes-ab equations are solved on coarse grids at various levels and the Navier–Stokes equations are solved on the ‘‘nest grid’’. The method uses Crank– Nicolson time-stepping for the viscous terms, explicit time-stepping for the remaining terms, and Richardson iteration to solve linear systems encountered at each time step and on each grid level. To explore the computational efficiency of the method, comparisons are made with results obtained from an analogous spectral multigrid method for the Navier–Stokes equations. These comparisons are based on computing work units and residuals for multigrid cycles. Most importantly, we examine how choosing different values of the length scales a and b entering the Navier–Stokes-ab equations influence the efficiency and accuracy of these multigrid schemes. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction Numerical simulations of turbulent fluid flows require the solution of large systems of equations on highly refined grids. However, direct solution methods are often impractical for realistically large problems and standard iterative methods have very slow convergence rates. Multigrid schemes can accelerate the convergence of iterative methods by considering low frequency errors on a fine grid as high frequency errors on a coarse grid. Spectral multigrid methods combine both the accuracy of spectral discretizations and the efficiency of multigrid algorithms. Spectral methods were first combined with multigrid techniques by Zang et al. [1,2] for self-adjoint elliptic equations with periodic or Dirichlet boundary conditions. Zang and Hussaini [3] later extended this approach to the time-dependent Navier–Stokes (NS) equations. This was followed by the work of Erlebacher et al. [4], who used spectral multigrid methods to study turbulence via three-dimensional large-eddy simulations. Improvements to the basic approach have been proposed for elliptic problems by Brandt et al. [5] and Heinrichs [6,7] and for the incompressible NS equations by Zhang et al. [8]. Preconditioning has been investigated by Heinrichs [9–11] and Phillips [12]. More recently, researchers have demonstrated that spectral and pseudo-spectral multigrid algorithms can be applied to various

⇑ Corresponding author. E-mail addresses: [email protected], [email protected] (E. Fried). 0045-7930/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2010.12.016

research fields such as steady and incompressible flows [13,14], viscoelastic turbulence [15], quantum mechanics [16], and optimal control [17]. The purpose of this paper is to develop improved spectral multigrid methods using recently developed regularizations of the NS equations for isotropic and homogeneous turbulent flow. Various regularizations [18–23] of the NS equations have been recently studied as subgrid models which resolve motions only above some critical scales. These models allow for accurate and efficient simulations of high Reynolds-number flows, while resolving only large-scale flow structures. One such model is the Navier– Stokes-a (NS-a) equations obtained by Chen et al. [18] by adding a viscous term to the Euler-a equations of Holm et al. [24]. In the context of Lagrangian averaging, a is the statistical correlation length of the excursion taken by a fluid particle away from its phase-averaged trajectory. Heuristically, a is the characteristic linear dimension of the eddy scale that the model can resolve accurately. The microscale of the NS-a model is much smaller than the length scale a, but, as Foias et al. [25] show, it is larger than that of the NS equations. Direct numerical simulation (DNS) by Chen et al. [26] has demonstrated that the NS-a equations captures the properties of flows for eddy scales greater than a. Importantly, Foias et al. [27] proved that solutions to the classical NS equations are recovered in the limit a ? 0. Fried and Gurtin [28,29] derived a slight generalization of the NS-a equations, called the Navier–Stokes-ab (NS-ab) equations. For a fluid with density q and kinematic viscosity m, these equations constitute a system

103

T.-Y. Kim et al. / Computers & Fluids 44 (2011) 102–110 @v @t

þ ðgrad v Þu þ ðgrad uÞ> v ¼ grad qp þ mð1  b2 DÞDu;

v ¼ ð1  a

2

DÞu;

div u ¼ 0;

) ð1Þ

for unfiltered and filtered velocities v and u and a filtered pressure p. Like a, b > 0 has dimensions of length. Setting b equal to a in (1) yields the NS-a equations. Using (1)2 to eliminate v from (1)1 yields a regularized NS equation for u; this equation contains a dispersive term, of energetic origin, with coefficient a, and a dissipative term with coefficient b. Phenomenologically, a represents eddy scales in the inertial range, where the dissipationless transfer of kinetic energy between intermediate scale eddies occurs, and b represents eddy scales in the dissipation range, where viscous damping converts the kinetic energy contained in the smallest eddies into heat. Accordingly, one expects that b < a. Since the specialization of the NS-ab equations to the NS-a equations requires equating b to a, the NS-a equations involve an implicit equating of disparate length scales. This is a vestige of the conventional derivation of the NS-a equations, a derivation which associates viscosity with the stretching tensor determined by the unfiltered velocity. Lundgren [30] introduced a strained spiral vortex model as an assembly of asymptotic solutions to the NS equations. For Lundgren’s model, it was found that the energy spectrum associated with axial vorticity obeys E(k)  k5/3, which conforms to Kolmogorov’s [31] similarity theory. Specifically, Lundgren’s model [30] consists of axial vorticity helically wound up around the Burgers’ vortex tube. This model was extensively discussed subsequently by Lundgren [32], Pullin and Saffman [33], and Pullin et al. [34]. Motivated by the results obtained for the NS equations, Chen and Fried [35] applied the spiral vortex model to the NS-a and NS-ab equations. The extension of the spiral vortex model to the NS-ab equations begins with the corresponding Helmholtz vorticity equation for the NS-ab equations. Specifically, taking the curl of (1)1 and using the identities

ðgrad v Þu þ ðgrad uÞT v ¼ gradðu  v Þ  u  curl v : curlðDuÞ ¼ Dðcurl uÞ

) ð2Þ

yields @q @t

þ ðgrad qÞu  ðgrad uÞq ¼ mð1  b2 DÞDx;

x ¼ curl u; q ¼ curl v ¼ ð1  a2 DÞx;

) ð3Þ

2. The spiral vortex model for the NS-ab equations In this section, we briefly describe a vorticity stream-function equation of the spiral vortex model for the NS-ab equations and the three-dimensional energy spectrum which can be calculated from solutions of the two-dimensional unstrained NS-ab equations. Further details are provided by Chen and Fried [35]. The vorticity stream-function for the spiral vortex model for the NS-ab equations in polar coordinates r and h is written as @Q @T

þ 1r

@ W

@Q @h @r

ð4Þ

Numerical studies by Kim et al. [23] show that the NS-ab equations at lower resolutions can sufficiently approximate the highly-resolved NS equations for simulations of three dimensional, isotropic, turbulent flows. These studies suggest that the NS-ab equations can be used to smoothen solutions at coarse grid levels in multigrid schemes. Furthermore, the spiral vortex model provides a computationally efficient means to compute the three-dimensional energy spectrum from solutions to the two-dimensional unstrained NS equations. In this paper, we present a spectral multigrid method for the spiral vortex models for the NS and NS-ab equations. The method involves solving the NS-ab equations on coarse grid levels and solving the NS equations on the finest grid level and uses Crank–Nicolson time-stepping for the viscous terms, explicit timestepping for the remaining terms, and Richardson iteration to solve linear systems encountered at each time step and on each grid level.

@Q @h

 @@rW



¼ mð1  b2 SDÞDx;

)

x ¼ DW; Q ¼ ð1  a2 SDÞx;

ð5Þ

where Q and x are the unfiltered and filtered vorticities, W is the stream-function, and S is related to the strained time T via S = 1 + bT with frequency b. For convenience, the vortex equations (5) are non-dimensionalized by introducing a characteristic length scale L, a characteristic velocity scale U, and a characteristic time scale L/U. For convenience, the symbols used for the original dimensional variables are retained. The dimensionless viscosity m then represents the inverse Reynolds number of the flow. Furthermore, the dimensionless length scales a and b can be obtained by dividing their dimensional counterparts by 2pL. Approximating the unfiltered velocity Q by a finite collection I P of Fourier modes Qk via Q ðx; TÞ ¼ I Q k eikx and substituting into (5) yields the pseudo-spectral equations dQ k dT

¼



@ W @Q @x @y

9 1þb2 k2 S 2 =  m 1þ 2 k Q k; 2 a k S k ; 2 Q k ¼ ð1 þ a2 k SÞxk ; @Q @x

 @@yW

xk ¼ k2 Wk ;

where q and x are the unfiltered and filtered vorticities. Based on these equations, Chen and Fried [35] derived approximate analytical spiral vortex solutions for both the NS-a and NS-ab equations. Furthermore, the rotational form of (1) using the identity (2)1 can also be obtained as

@v p  u  q ¼ grad þ mð1  b2 DÞDu: @t q

This paper is organized as follows: In Section 2, we present a vorticity stream-function formulation of the spiral vortex model for the NS-ab equations. This formulation is used to develop the spectral multigrid schemes. In Section 3, a time-stepping algorithm and a relaxation scheme are described. In Section 4, spectral multigrid algorithms using the vorticity stream-function equation of the spiral vortex model are explained in detail. In Section 5, the accuracy and efficiency of these multigrid schemes are investigated. Finally, in Section 6, we present a summary and concluding remarks.



ð6Þ

where k = jkj is the wavenumber. Setting b equal to a in (6)1 gives the pseudo-spectral equations for the NS-a equations. The spiral vortex model makes it possible to calculate of the energy spectrum from solutions to the two-dimensional unstrained NS-ab equations. The expression for the energy spectrum is given by

Ea ðk; TÞ ¼

Z

C 2

2 k2 Þ2

k ð1 þ a

Tc

pffiffiffiffiffiffiffiffiffi SðTÞF 2d ðSðTÞ1=2 k; TÞdT;

ð7Þ

0

where C is a constant and F2d represents the two-dimensional unfiltered enstrophy spectrum of the unstrained flow, defined by

F 2d ðk; TÞ ¼ k

Z 2p

^ ðk cos h; k sin h; TÞj2 dh jQ

ð8Þ

0

with

b ðkx ; ky ; TÞ ¼ Q

1 ð2pÞ2

Z

Q ðx; y; TÞeiðkx xþky yÞ dx dy

A

the Fourier integral of Q.

ð9Þ

104

T.-Y. Kim et al. / Computers & Fluids 44 (2011) 102–110

3. Time-stepping algorithm Applying Crank–Nicolson for the viscous terms and forward Euler1 for the remaining terms yields a simple time discrete version of (5), which can be written as

Q nþ1  Q n ¼ DT

 n @ W @Q @ W @Q 1  þ ½mð1  b2 SDÞDxn 2 @x @y @y @x

þ mð1  b2 SDÞDxnþ1 ;

ð10Þ

where DT represents a time increment. (This use of the symbol D should not be confused with its use to denote the Laplace operator.) Upon introducing

F n ¼ Q n þ DT

 n @ W @Q @ W @Q DT þ mð1  b2 SDÞDxn :  2 @x @y @y @x

ð11Þ

cycle multigrid process for solving only the NS-ab equations by considering the interplay between two grids. This method can be represented by the following procedure: 1. Solve (or relax) (12) n times using the Richardson iteration (15) on Xh to obtain the filtered vorticity xh. 2. Compute the residual rh using (16) and restrict it to the h coarse grid residual r2h ¼ R2h h r . 3. Solve (12) by replacing its right-hand side with r2h on X2h and obtain the filtered vorticity x2h. 4. Update the fine grid approximation using xh xh þ Ph2h x2h . 5. Solve (or relax) (12) n times with initial guess xh from the previous step.

(10) becomes

Q nþ1 

DT mð1  b2 SDÞDxnþ1 ¼ F n : 2

ð12Þ

Importantly, the Crank–Nicolson time-stepping algorithm for the NS equations is recovered by taking a = b = 0 in (12) to yield

Q nþ1 

DT mDQ nþ1 ¼ F n ; 2

ð13Þ

where

F n ¼ Q n þ DT

 n @ W @Q @ W @Q DT  þ m DQ n : 2 @x @y @y @x

ð14Þ

We use the Richardson iteration scheme as a relaxation strategy to solve the linear systems encountered at each time step and on each grid level. Applied to either (12) or (13), the solution after one smoothing step becomes

Q nþ1

Q n  sr n ;

ð15Þ

where s is the relaxation parameter and the residual rn is defined as

rn ¼ F n  Q n þ

DT mð1  b2 SDÞDxn 2

ð16Þ

for the NS-ab equations and as

rn ¼ F n  Q n þ

DT m DQ n 2

ð17Þ

for the NS equations. Both stationary (fixed) and non-stationary (variable) relaxation parameters can be considered. Experiments here were mostly performed with a non-stationary Richardson smoothing algorithm [4]. In this case, the jth relaxation parameter of a n-parameter cycle is

snj ¼

2=kmin ; ðj  1Þ cosðð2j  1ÞpÞ=2nÞ þ j þ 1

j ¼ 1; . . . ; n;

ð18Þ

where j is the condition number

j¼

kmax kmin

Secondly, we describe a new multigrid algorithm using both the NS and NS-ab equations for a V-cycle with three grid levels. This method involves solving the NS equations on the finest grid and the NS-ab equations on coarse grids at various levels. For convenience, we refer to this method as the coupled spectral multigrid (CSM) algorithm. The CSM algorithm can be represented by the following procedure:

ð19Þ

and kmax and kmin are, respectively, the maximum and minimum eigenvalues of the linear operator underlying the system of Eq. (12). The special case of a stationary relaxation parameter arises on taking j = 1 in (18). 4. Spectral multigrid schemes In this section, the spectral multigrid algorithms that we developed are described. For convenience, we use Xh to denote the computational domain with the grid size h. First, we describe the V1 The second-order Adams–Bashforth (AB2) scheme can also be effective for the nonlinear term in (10).

1. Solve (or relax) (13) n times using the Richardson iteration (15) on Xh to obtain the unfiltered vorticity Qh. 2. Compute the fine residual rh using (17) and restrict it to the h coarse grid residual r4h ¼ R4h h r . 3. Solve (12) by replacing its right-hand side with r4h on X4h and obtain the filtered vorticity x4h. 4. Compute the residual r4h using (16) and restrict it to the 4h coarse grid residual r8h ¼ R8h 4h r . 5. Solve (12) by replacing its right-hand side with r8h on X8h and obtain the filtered vorticity x8h. 6. Update the fine grid approximation using 8h x4h x4h þ P4h x . 8h 7. Solve (or relax) (12) n times with initial guess x4h of the previous step. 8. Update the fine grid approximation using Qh Q h þ Ph4h x4h . 9. Solve (or relax) (13) n times with initial guess Qh from the previous step.

Here, n represents the number of relaxations for the Richardson iteration. The CSM algorithm is graphically described in Fig. 1. Moreover, R2h h is the operator that yields the restriction from the computational domain Xh with grid size h to the computational domain X2h with coarser grid size 2h. On the other hand, the prolongation operator P h2h represents the grid transfer from X2h with h grid size 2h to Xh with grid size h. The operators R4h h and P 4h have analogous interpretations. In the spectral multigrid method, the restriction of a variable from a fine grid to a coarse grid can be achieved by representing the variable in Fourier space, removing the highest modes not resolvable on the coarse grid, and then transforming back to physical space on the coarse grid. Similarly, the prolongation of a variable from a coarse grid to a fine grid can be performed by representing the variable in Fourier space, adding the highest modes resolvable on the fine grid into the Fourier series with zero coefficients, and then transforming back to physical space on the fine grid. The algorithms presented above can be easily extended to spectral multigrid schemes with multiple V-cycles involving more than

105

T.-Y. Kim et al. / Computers & Fluids 44 (2011) 102–110 2

X0 ðrÞ ¼ C

ð1 þ hr =r 20 Þ expðr2 =r 20 Þ ð1 þ hÞpr 2

C

1 þ h  ð1 þ h þ hr =r 20 Þ expðr2 =r 20 Þ : ð1 þ hÞpr 3

ð21Þ

The vorticity Q0 in the central core is given by 2

Q 0 ðr; tÞ ¼ C

ð1 þ hr =r 20 Þ expðr 2 =r 20 Þ ð1 þ hÞp 4

 4a2 SðtÞC

ðhr =r40  ð3h  1Þr2 =r 20 þ h  1Þ expðr2 =r 20 Þ : ð1 þ hÞp ð22Þ

Fig. 1. Graphical description of a CSM multigrid V-cycle for four grid levels. The choices of a and b for the NS-ab equations can be different at different coarse-grid levels.

The requirement that h obey 0 < h < 1 ensures that X is a monotonically decreasing function of its argument. Following Lundgren [32], h is set to 1/2 and C and r0 are set to unity. Importantly, the spiral solution for the NS equations can be recovered by taking b/a = 1 in (20)1 and a = 0 in (22) and is used as an initial condition for the CSM algorithm. 5.1. Multigrid scheme based on the NS-ab equations

two or three grid interplays. The main difference between the two algorithms is evident from Steps 3 and 4. In these steps, the CSM method solves the NS-ab equations on the coarse grids and considers the interplay between X4h and Xh, instead of between X2h and Xh. For example, if we consider a V-cycle with the finest grid 322, the NS-ab or NS algorithm solves the NS or NS-ab equations at 322,162, 82, and 42 resolutions. On the other hand, the CSM algorithm solves the NS equations at 322 and then the NS-ab equations at 82 and 42 resolutions without solving the NS or NS-ab equations at 162 resolution. This choice is based on the knowledge that the NS-ab equations can efficiently represent the highly-resolved NS equations at lower resolutions granted appropriate choices of the length scales a and b, as discussed by Kim et al. [23]. The computational time required for the CSM method is obviously less than that required for either the NS or the NS-ab algorithm.

5. Numerical study We examine the computational efficiency of our multigrid schemes for homogeneous and isotropic turbulent flows. A square domain with side length 2p and periodic boundary conditions is considered. Aliasing errors are removed using full dealiasing, as described by Canuto et al. [36], and all simulations use the spiral solution (20) as an initial condition. The spiral vortex model analogous to that provided by Lundgren [32] for the NS equations has the form of a two-sided roll-up with two halo spirals of amplitude f, with separation p, plus an additional central core. The dimensionless functions Q, f, and X take the form

Q ðr; h; tÞ ¼ Q 0 ðr; tÞ þ 2f ðrÞ

1 X

cosf2nðh  XðrÞtÞg exp

n¼1

( ) 4b2 n2 X0 ðrÞ2 mt3 ; 3a2

2

hr =r20 expðr 2 =r 20 Þ f ðrÞ ¼ C ; 2pð1 þ hÞr 20

We first begin by studying the spectral multigrid scheme of the NS-ab equations. In Fig. 2, vorticity contours with different resolutions are provided for a = 1/48, b = 1/64, and C/m = 25,000. We emphasize that the values of a and b are normalized by 2pL hereinafter. The results indicate that the NS-ab vorticity contours at both 1282 and 2562 resolutions are qualitatively quite close to the NS vorticity contour at 10242 resolution. Fig. 3 shows the normalized energy spectra with different resolutions for three cases of a and b: a = b = 0 at 5122 resolution, a = 1/ 48, b = 1/256 at 1282 and 642 resolutions. The energy spectra are computed via the enstrophy spectrum (8) using a fast Fourier transform algorithm and summing over time. We note that the energy spectra show high peaks due to high enstrophy peaks at low wavenumbers, as evidenced in Fig. 6. The NS-ab energy spectra follow Kolmogorov’s 5/3 law in the inertial range and are seen to be consistent with the NS spectrum in the dissipation range at resolutions of both 642 and 1282. These results are consistent with those obtained by Kim et al. [23] for three-dimensional isotropic and homogeneous turbulence flow in a periodic cube. Both results imply that the highly-resolved NS results can be well approximated by NS-ab equations at lower resolutions, provided that appropriate values of length scales a and b are selected. Furthermore, these results suggest that the NS-ab equations can be used to smooth results obtained at coarse grid levels in multigrid schemes. We now consider the sensitivity of the energy spectrum to variations of a and b in Fig. 4. Plots of energy spectra are provided, for 1282 resolution and Reynolds number C/m = 2500. The results show that for a fixed, the change of Kolmogorov’s 5/3 law in the inertial range is almost indistinguishable for different choices of b. As b is increased, the energy spectrum in the dissipation range becomes steeper. On the other hand, for b fixed, as a is decreased, Kolmogorov’s 5/3 law is better preserved in the inertial range, leaving the slope of the energy spectrum in the dissipation range unchanged. These results clearly verify that decreasing or increasing a influences the inertial range while variations of b influence the dissipation range.

2

XðrÞ ¼ C

1 þ h  ð1 þ h þ hr =r 20 Þ expðr 2 =r 20 Þ ; 2pð1 þ hÞr 2

ð20Þ

where r0 is a dimensionless radial value and h is a constant. The derivative of the average angular velocity X with respect to r is given by

5.2. Multigrid scheme based on both the NS and NS-ab equations In this section, we study the efficiency of the CSM algorithm combining the NS and NS-ab equations. More importantly, we examine how choosing different values of a and b influences the

106

T.-Y. Kim et al. / Computers & Fluids 44 (2011) 102–110

Fig. 2. Vorticity contours of the NS-ab at lower resolutions and the NS equations at higher resolution at T = 100 (a = 1/48, b = 1/64, and C/m = 25,000).

solving the NS equations at 2562 resolution and the NS-ab equations at coarse grid levels 642,322, 162,82 and 42 for C/m = 25,000. The enstrophy spectrum is obtained from (8) by taking b/a = 1 at 2562 resolution. For the NS-ab equations, we take length scales a = 1/48 and b = 1/68. These results are qualitatively close to those obtained by Lundgren [32] and confirm that our multigrid algorithm can be used to simulate the spiral vortex model. Now, we study the efficiency and accuracy of the CSM algorithm using the spectral multigrid method based on the NS equations as a basis for comparison. There are many different ways to measure the efficiency of a multigrid algorithm. Ultimately, the user is interested in the total CPU clock-time a code takes to complete execution. However, this timing is strongly computer-dependent, and on a given system, programming skill can greatly influence the results. We study the efficiency of our algorithm based on work of Erlebacher et al. [4]. They calculated the ratio of norms of the residual after and before a single V-cycle, and defined the convergence rate l 1 per V-cycle as Fig. 3. Comparison of energy spectra of the NS-ab equations at lower resolutions with the NS equations simulation at higher resolution (C/m = 2500, a = 1/48, and b = 1/256).

efficiency and accuracy of the CSM scheme. To do so, we first define the relative norm of the residual as

krk ¼

krexit k ; krentry k

ð23Þ

l 1 ¼ krk1=Nf ;

ð24Þ

where Nf is the total number of the fine-grid sweeps during on V-cycle. This measure of efficiency does not incorporate CPU time. Another way to measure of the overall algorithm efficiency is via the convergence rate

l 2 ¼ krk1=Nt ;

ð25Þ

where Nt is defined as

where rexit and rentry represent the residual vector after and before a single V-cycle. Unless stated otherwise, we use the maximum norm and 1014 as tolerance for the convergence criteria. To verify our algorithm, vorticity contours and the enstrophy spectrum are determined. Figs. 5 and 6 show results obtained by

(a)

Nt ¼

CPU time per V  cycle : CPU time per fine grid relaxation

ð26Þ

In (24) and (25), krk is as defined previously in (23). The conver 2 can measure the decrease in the residual norm per gence rate l

(b)

Fig. 4. Influence of (a) b with a = 1/48 (b) a with b = 1/128 for the NS-ab equations at 1282 resolution and C/m = 25,000.

107

T.-Y. Kim et al. / Computers & Fluids 44 (2011) 102–110

Fig. 5. Vorticity contours obtained from the CSM algorithm with a = 1/48, b = 1/68, and C/m = 25,000 at 2562 fine grid-resolution.

Fig. 6. Two-dimensional enstrophy spectrum versus wavenumber obtained from the CSM algorithm with a = 1/48, b = 1/68, and C/m = 25,000 at 2562 fine grid-resolution.

fine sweep taking the total multigrid overhead into account. Furthermore, the performance of the algorithm can be efficiently ascertained by considering the CPU time and the residual norm together per V-cycle.  2 for the NS and CSM algoTable 1 shows the convergence rate l rithms. Grid resolutions of 2562, 5122, and 10242 are considered for C/m = 25,000. For the CSM algorithm, a = 1/48 and b = 1/64. The difference between the NS and CSM algorithms is explained in Section 4. The results show that the convergence rates for the CSM method are faster than those of the NS algorithm. The CSM method is therefore more efficient than the NS multigrid method. The values of the convergence rates are quite small because the residuals converge rapidly to machine precision. This is because the coefficients of the system of Eqs. (12) and (13) are constant and thus they are almost exact with a small number of relaxations.  1 with the change Tables 2 and 3 display the convergence rate l of a for a given b and the change of b for a given a, respectively. We solve the NS equations at 2562 resolution and the NS-ab equations at coarse-grid levels 642, 322, 162, 82, and 42 for C/m = 25,000 with n = 2. The results show that the convergence rate decreases as a decreases. This suggests that our algorithm becomes more efficient with decreasing a. On the other hand, Table 3 shows that decreasing b does not alter the convergence rate at fixed a. Fig. 7 shows CPU time units versus a change of grid-resolution for one V-cycle of the NS and CSM algorithms. We tested fine grid resolutions 1282, 2562, 5122, and 10242. For the CSM algorithm, a = 1/48 and b = 1/64. The CPU time units include the time required for iterations on the fine and coarser grids, as well as for inter-grid transfers. The CSM algorithm reduces the CPU time units of the NS algorithm by a factor of approximately 1/3 and can significantly

Table 1  2 (C/m = 25,000, a = 1/48, b = 1/64). Convergence rate l Resolution (# of relaxations)

NS

CSM

2562(3) 5122(4) 10242(5)

1.38e11 1.75e11 1.80e11

2.51e13 3.70e13 2.46e12

Table 2  1 for the CSM method with change of Convergence rate l a for a given b = 1/64 at 2562 fine grid-resolution (C/ m = 25,000, n = 2).

a

l 1

1/8 1/16 1/32 1/48

1.90e3 1.74e3 1.55e3 1.45e3

Table 3  1 for the CSM method with change of Convergence rate l b for a given a = 1/8 at 2562 fine grid-resolution (C/ m = 25,000, n = 2). b

l 1

1/8 1/16 1/32 1/64

1.90e3 1.90e3 1.90e3 1.90e3

save computational time with increasing grid resolution for one V-cycle. Fig. 8 shows CPU time units with increasing time T for NS and CSM algorithms for 2562 and 5122 resolutions and C/m = 25000. We take the time increment to be DT = 0.01. In contrast to the studies reported above, we now consider the computational efficiency of the CSM algorithm for multiple V-cycles. At each time step, we start computing a V-cycle with a certain number of relaxations for the Richardson iteration and check if the residual norm (23) is convergent within a tolerance of 1013. If it is not convergent, another V-cycle is computed with an increased number of relaxations. We continue this process until the residual norm is no greater than the imposed tolerance. The CSM and NS algorithms typically converge within two to four relaxations. The CPU time measures presented here include all of these processes in a multigrid cycle at each time step. The results show that the CSM algorithm reduces the CPU time of the NS algorithm by a factor of

108

T.-Y. Kim et al. / Computers & Fluids 44 (2011) 102–110

Fig. 7. CPU time units with various grid resolutions between the NS and CSM algorithms for C/m = 25,000, a = 1/48, and b = 1/64.

ative errors with changes of b for a = 1/16, 1/32, 1/48, and 1/60 are provided. In particular, the cases a = 1/48 and 1/60 are displayed in the figure. The results show that there exists an optimal choice of a and b that, based on the norm (27), exhibits the smallest deviation from solutions to the NS equations at higher resolution. The right panel of the figure shows the relative errors at the optimal values of a and b as a function of a. The plot shows that the error decreases with decreasing a. This is consistent with the observation that the NS-a equations approach to the NS equations as a decreases. Based  1 is determined for both staon these results, the convergence rate l tionary and non-stationary relaxation parameters, as shown in Table 4. The convergence rate decreases as the relative error (27) decreases. This result suggests that appropriate choices of a and b can increase the efficiency of the CSM algorithm. To check the efficiency of the CSM algorithm for different choices of a and b at different coarse-grid levels, we consider the  1 (Table 5). We solve the NS equations at convergence rate l 10242 resolution and the NS-ab equations at coarse-grid levels of 2562 and 1282. As shown in Fig. 9, optimal values of b were obtained for a = 1/16, 1/32, 1/48, and 1/60 at resolutions of 2562 and 1282. These values of a and b are used for NS-ab at coarse-grid  1 is provided for three cases of a and b: (i) levels. In particular, l a = 1/32, b = 1/76(2562) and a = 1/16, b = 1/44 (1282), (ii) a = 1/48, b = 1/128(2562) and a = 1/32, b = 1/84(1282), and (iii) a = 1/60,  1 for b = 1/164(2562) and a = 1/48, b = 1/116(1282). The value of l the same choices of a = 1/16 and b = 1/32 at both resolutions is also  1 for different choices of provided. The results show that values of l a and b at different grid levels are faster than that for fixed values of a and b at coarse-grid levels. Hence, the efficiency of the CSM algorithm can be improved by taking different values of the length scales a and b at different grid levels.

6. Summary

Fig. 8. CPU time units with a change of time between the NS and CSM algorithms for C/m = 25,000, a = 1/48, and b = 1/64.

approximately 1/3 at 5122 resolution. These results suggest that the CSM algorithm can reduce significantly the computational time of the NS algorithm as time increases. Fig. 9 shows the result of a numerical study of potential optimal choices for the parameters a and b, with the criterion for optimality being to minimize the deviation between the vorticity field calculated from the NS-ab equations and that calculated from the NS equations. Solutions to the NS equations on a grid resolution of 10242 are used as the base solution to calculate the error. The relative norm E between the coarse and fine grids is defined as

E¼

jjQ f ðiÞ  Q c ðiÞjj ; jjQ f ðiÞjj

ð27Þ

ffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi PNc 2 where jjQ f ðiÞjj ¼ i jQ f ðiÞj , with Nc being the number of coarse grid points. Here, Qc and Qf represent the unfiltered vorticity fields obtained at 2562 and 10242 resolutions, respectively. The computation is achieved by obtaining Qf at coarse grid points using bicubic interpolation of the simulation results at the 10242 resolution. Rel-

In this paper, the spectral multigrid algorithms for periodic homogeneous and isotropic turbulence in two space-dimensions that use the NS-ab equations are developed. These algorithms are based on the vorticity/stream-function equation of the spiral vortex model for the NS and NS-ab equations. First, we presented the V-cycle multigrid algorithm for the NSab equations and studied the influence of length scales of a and b on the energy spectrum. Importantly, numerical results show that the highly-resolved NS results can be well represented by the NSab equations at lower resolutions, provided that appropriate values of length scales a and b are selected. These results suggest that the NS-ab equations can be used to smooth results obtained at coarse grid levels in multigrid schemes. Second, the new spectral multigrid algorithm which solves the NS equations on a fine grid and the NS-ab equations on coarser grid levels was developed. We checked the efficiency of the algorithm by computing the convergence rates and computational time for a V-cycle. Numerical results show that the convergence rates of this algorithm for the same number of relaxations are faster than those of the conventional NS multigrid scheme. Furthermore, this algorithm can result in more computational savings than traditional multigrid schemes based on the NS equations with increasing grid resolutions and time steps. We also examined how the length scales a and b influence the efficiency and accuracy of the CSM algorithm. The results show that the CSM algorithm converges rapidly with decreasing a for fixed b, but that decreasing b does not change the convergence rate at fixed a. Numerical experiments show that our multigrid algorithms can be used to accelerate convergence of a spectral multigrid method of the NS equations. Based on these algorithms, we expect that it should be possible to develop an efficient spectral multigrid

T.-Y. Kim et al. / Computers & Fluids 44 (2011) 102–110

(a)

109

(b)

Fig. 9. (a) Relative error as a function of b for a = 1/48, 1/60. (b) Relative error at optimal values of b for a = 1/16, 1/32, 1/48, 1/60 as a function of a (2562 resolution and C/ m = 25,000).

Table 4  1 at optimal values of length scales b and a for the CSM method at Convergence rate l 2562 fine grid-resolution (C/m = 25,000, n = 2).

a

b

Stationary

Non-stationary

1/16 1/32 1/48 1/60

1/28 1/76 1/128 1/164

4.25e3 3.89e3 3.72e3 3.65e3

1.74e3 1.55e3 1.45e3 1.41e3

Table 5  1 at the same (a = 1/16, b = 1/32) and different values of b and a on Convergence rate l coarse-grid levels of 2562 and 1282 for the CSM method at 10242 fine grid-resolution (C/m = 25,000, n = 2). 2562

1282

l 1

a

b

a

b

1/16 1/32 1/48 1/60

1/32 1/76 1/128 1/164

1/16 1/16 1/32 1/48

1/32 1/44 1/84 1/116

3.02e2 2.74e2 2.55e2 2.44e2

algorithm for three-dimensional periodic homogeneous turbulence. An open question concerns the optimal choice of b given a as a function of grid resolution and Reynolds number, as noted by Kim et al. [23]. We believe this question can be addressed by studying a similarity theory of the NS-ab model like that used by Layton and Neda [37,38] to study time-relaxation and deconvolution models for turbulence. This is a subject for future research. Acknowledgements The work of EF was supported by the Department of Energy under Award Number DE-SC0004604. References [1] Zang TA, Wong YS, Hussaini MY. Spectral multigrid methods for elliptic problems. J Comput Phys 1982;48:485–501. [2] Zang TA, Wong YS, Hussaini MY. Spectral multigrid methods for elliptic problems II. J Comput Phys 1984;54:489–507. [3] Zang TA, Hussaini MY. On spectral multigrid methods for the time-dependent Navier–Stokes equations. Appl Math Comput 1986;19:359–72. [4] Erlebacher G, Zang TA, Hussaini MY. Spectral multigrid methods for the solution of homogeneous turbulence problems. Lect Notes Pure Appl Math 1988;110:177–94. [5] Brandt A, Fulton SR, Taylor GD. Improved spectral multigrid methods for periodic elliptic problems. J Comput Phys 1985;58:96–112.

[6] Heinrichs W. Line relaxation for spectral multigrid methods. J Comput Phys 1988;77:166–82. [7] Heinrichs W. Multigrid methods for combined finite difference and Fourier problems. J Comput Phys 1988;78:424–36. [8] Zhang W, Zhang CH, Xi G. An explicit Chebyshev pseudospectral multigrid method for incompressible Navier–Stokes equations. Comput Fluids 2010;39:178–88. [9] Heinrichs W. A 3D spectral multigrid method. Appl Math Comput 1991;41:117–28. [10] Heinrichs W. Spectral multigrid techniques for the Navier–Stokes equations. Comput Method Appl Mech Eng 1993;106:297–314. [11] Heinrichs W. Finite element preconditioning for spectral multigrid method. Appl Math Comput 1993;59:19–40. [12] Phillips TN. Relaxation schemes for spectral multigrid methods. J Comput Appl Math 1987;18:149–62. [13] Chou MH. A multigrid pseudospectral method for steady flow computation. Int J Numer Methods Fluids 2003;43:25–42. [14] Krastev K, Schäfer M. A multigrid pseudo-spectral method for incompressible Navier–Stokes flows. C R Mecanique 2005;333:59–64. [15] Housiadas KD, Wang L, Beris AN. A new method preserving the positive definiteness of a second order tensor variable in flow simulations with application to viscoelastic turbulence. Comput Fluids 2010;39:225–41. [16] Borzi A, Hohenester U. Multigrid optimization schemes for solving Bose– Einstein condensate control problems. SIAM J Sci Comput 2007;30:441–62. [17] Gong Q, Fahroo F, Ross IM. Spectral algorithm for pseudospectral methods in optimal control. AIAA J Guid Control Dynam 2008;31:460–71. [18] Chen S, Foias C, Olson E, Titi ES, Wynne S. The Camassa–Holm equations as a closure model for turbulent channel and pipe flow. Phys Rev L 1988;81:5338–41. [19] Foias C, Holm DD, Titi ES. The Navier–Stokes-alpha model of fluid turbulence. Phys D 2001;152–153:505–19. [20] Cao C, Holm DD, Titi ES. On the Clark-a model of turbulence global regularity and long-time dynamics. J Turbulence 2005;6:N20. [21] Geurts BJ, Holm DD. Regularization modeling for large-eddy simulation. Phys Fluids 2003;15:L13–6. [22] Geurts BJ, Holm DD. Leray and LANS-a modelling of turbulent mixing. J Turbulence 2006;7:1–33. [23] Kim TY, Cassiani M, Albertson JD, Dolbow JE, Fried E, Gurtin ME. Impact of the inherent separation of scales in the Navier–Stokes-ab equations. Phys Rev E 2009;79:045307. [24] Holm DD, Marsden JE, Ratiu T. Euler–Poincaré models of ideal fluids with nonlinear dispersian. Phys Rev L 1998;80:4273–7. [25] Foias C, Holm DD, Titi ES. The Navier–Stokes–alpha model of fluid turbulence. Phys D 2001;152–153:505–19. [26] Chen S, Holm DD, Margolin LG, Zhang R. Direct numerical simulations of the Navier–Stokes alpha model. Phys D 1999;133:66–83. [27] Foias C, Holm DD, Titi ES. The three dimensional viscous Camassa–Holm equations, and their relation to the Navier–Stokes equations and turbulence theory. J Dyn Differ Equat 2002;14:1–35. [28] Fried E, Gurtin ME. Turbulence kinetic energy and a possible hierarchy of length scales in a generalization of the Navier–Stokes-a theory. Phys Rev E 2007;75:056306. [29] Fried E, Gurtin ME. A continuum mechanical theory for turbulence: a generalized Navier–Stokes-a equation with boundary conditions. Theoret Comput Fluid Dyn 2008;22:433–70. [30] Lundgren TS. Strained spiral vortex model for turbulent fine structure. Phys Fluids 1982;25:2193–203. [31] Kolmogorov AN. The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers. Dokl Akad Nauk SSSR 1941;30:301–5. [32] Lundgren TS. A small-scale turbulence model. Phys Fluids A 1993;5:1472–83.

110

T.-Y. Kim et al. / Computers & Fluids 44 (2011) 102–110

[33] Pullin DI, Saffman PG. On the Lungren–Townsend model of a turbulent fine scales. Phys Fluids A 1993;5:126–45. [34] Pullin DI, Buntine JD, Saffman PG. On the spectrum of a stretched spiral vortex. Phys Fluids 1994;6:3010–27. [35] Chen X, Fried E. The influence of the dispersive and dissipative scales a and b on the energy spectrum of the Navier–Stokes-ab model for turbulent flow. Phys Rev E 2008;78:046317.

[36] Canuto C, Hussaini MY, Quarteroni A, Zang TA. Spectral methods in fluid dynamics. Springer-Verlag; 1988. [37] Layton W, Neda M. Truncations of scales by time relaxation. J Math Anal Appl 2007;325(2):788–807. [38] Layton W, Neda M. A similarity theory of approximate deconvolution models of turbulence. J Math Anal Appl 2007;333(1):416–29.