Variable scale filtered Navier-Stokes Equations. A new procedure to deal with the associated commutation error
aa r X i v : . [ phy s i c s . f l u - dyn ] N ov Variable scale filtered Navier-Stokes Equations. A new procedure to deal with theassociated commutation error
M. Iovieno
Doctorate Program in Fluid Dynamics at the Politecnico di Torino,Corso Duca degli Abruzzi 24, 10129 Torino, Italy.
D. Tordella
Dipartimento di Ingegneria Aeronautica e Spaziale,Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy.
A simple procedure to approximate the noncommutation terms that arise whenever it is necessaryto use a variable scale filtering of the motion equations and to compensate directly the flow solutionsfrom the commutation error is here presented. Such a situation usually concerns large eddy simu-lation of nonhomogeneous turbulent flows. The noncommutation of the average and differentiationoperations leads to nonhomogeneous terms in the motion equations, that act as source terms of in-tensity which depend on the gradient of the filter scale δ and which, if neglected, induce a systematicerror throughout the solution. Here the different noncommutation terms of the motion equationare determined as functions of the δ gradient and of the δ derivatives of the filtered variables. It isshown here that approximated noncommutation terms of the fourth order of accuracy, with respectto the filtering scale, can be obtained using series expansions in the filter width of approximationsbased on finite differences and introducing successive levels of filtering, which makes it suitable touse in conjunction with dynamic or mixed subgrid models. The procedure operates in a way whichis independent of the type of filter in use and without increasing the differential order of the equa-tions, which, on the contrary, would require additional boundary conditions. It is not necessary tointroduce a mapping function of the nonuniform grid in the physical domain into a uniform gridin an infinite domain. A priori tests on the turbulent channel flow ( Re τ = 180 and 590) highlightthe approximation capability of the present procedure. A numerical example is given, which drawsattention to the nonlocal effects on the solution due to the lack of noncommutation terms in themotion equation and to the efficiency of the present procedure in reducing the commutation erroron the solution. I. INTRODUCTION
The problem of the non commutativity of the fil-tering operation has been considered by Ghosal andMoin (1995)[1], van der Ven (1995)[2], Fureby and Ta-bor (1997)[3], Vasilyev et al. (1998)[4] and reviewed byGhosal (1999)[5] and Sagaut in his monography on largeeddy simulation (LES) for incompressible flows (2001)[6].Ghosal and Moin[1] showed that the commutation er-ror is of the second order in the filter width. Introducinga filter definition built on the mapping function of the nonuniform grid, they proposed a procedure that can be usedwhenever numerical schemes based on pseudo-spectralmethods or on finite differencing of an order higher thanthe second order are employed. The non commutationterm are expanded in a Taylor series in the filter width,where the coefficients depend on the spatial derivativesof the filtered field and the mapping function. In thisway, extra terms appear in the filtered equations, whichincrease the differential order of the equations. The au-thors suggested both the use of additional boundary con-ditions to mantain the well-posedeness of the problem orthe use of asymptotic expansions of the filtered variables,in terms of the square of the filter width, which howeverrequires the solution of additional non homogeneous per-turbative problems.A family of one parameter filters commuting with dif- ferentiation up to any given order in the filter width,which is assumed nonuniform in the integration domain,was constructed by van der Ven [2]. If a discretizationscheme of a given order is adopted in a LES, one mayselect a filter inside the family so that the lack of com-mutation between differentiation and filtering can be ne-glected.A general formulation of the commutation error, due toa non uniform filter and to the presence of boundaries,has been proposed by Fureby and Tabor[3]. Boundarydomain terms are explicitly formulated inside this repre-sentation. A detailed numerical analysis of the field dis-tribution of the intensity of the non commutation terms isgiven in this paper by comparing LES and direct numer-ical simulation (DNS) data obtained from simulations ofthe incompressible turbulent channel flow at Re τ = 180and 395. It has been found that the local intensity canbe as high as 21% of the local advection, with a fieldvolume averaged relative intensity of about 8%. Lowervalues apply if reference is made to the sum of the localadvection, pressure gradient, and molecular viscous fluxterms. An interesting result is that the use of differentsubgrid scale models negligibly affects the local and av-erage values of the sum of the commutation error terms.As expected, the high relative intensities of the commu-tation terms are concentrated in the flow regions wherethe gradient of δ is high. However, the question of thepossible extension of the effects on the macroscopic scaleof the flow is still left open.A generalization of the procedures proposed by Ghosaland Moin and van der Ven is presented in the paper byVasilyev et al. [4]. A minimization of the commutationerror is achieved by using a class of filters with n − n is the order of the employed nu-merical discretization scheme. The authors also supply agroup of rules to construct discrete filters that commuteswith differentiation up to any given order inside complexdomains.The method here proposed relies on an approximationof the specific non commutation term that corresponds tothe different terms of the motion equations. A commuta-tion approximation of the fourth order in the filter widthcan be obtained thanks to the introduction of successivelevels of average. See Sec. II for the basic formulationrelevant to an isotropic grid stretching and the Appendixfor the more general anisotropic case, which also speci-fies the formulation that is appropriate to wall-boundedflows.While performing large eddy simulations, the presentapproach can conveniently be used together with sub-grid models based on analog multi-level filtering, e.g.models which apply the dynamic procedure, Germano et al. (1991)[7], Germano (1992)[8] or the Bardina mixedmodel (1980)[9]. The filtering approach we use in thispaper is that of the very fundamental volume average,first applied by Smagorinsky in 1963[10]. The volumeaverage formulation is advantageous because it does notintroduce an error associated to domain boundaries, thusavoiding the problem of the addition of further non com-mutation terms in the equations.The variable scale filtered Navier-Stokes equations, in-cluding the commutation terms approximation, are givenin Sec. II A.A priori tests on the turbulent channel flow data basesby Alfonsi et al. (1998)[11] and Passoni et al. (1999)[12], Re τ = 180, and by Moser et al. (1999)[13], Re τ = 590,are presented in Sec. III A. An example of application ofthe numerical procedure is presented in Sec. III B, whichfocuses attention on the fact that this systematic erroris important throughout the entire flow and not only inthe regions where the non homogeneous terms of the mo-tion equation, which originate from the lack of commuta-tion of the operations of differentiation and filtering, aredifferent from zero. The capacity of the present proce-dure to reduce the relevant absolute and relative errorsis shown.Before proceeding to the other sections, it is necessaryto open a digression on the terminology adopted in whatfollows. Since, among the points being discussed in thepaper, there are: the structure of the motion equation,once a variable scale filtering is used, and the role playedby the terms which originate from the non commutationfilter-differentiation, we have had to clearly distinguishthe concept of commutation error on the flow solutionfrom that of non commutation term in the equations.By commutation error we mean the error which affects the flow solution, when a variable scale filter is used, butthe equations are used as if the operations of filtering anddifferentiation commute. By non commutation term wemean any of the terms which originate in the equationof motion when the filter scale is a function of the point.It is necessary to recall that in previous literature thelatter was called commutation error, since the equationwere always used as if they were commutative and theomissiom of the terms, which should make them completewhen the filter length varies, introduced the error in thesolution. II. NON COMMUTATION TERMS AND THEIRAPPROXIMATION
The loss of the commutation between the spatial fil-tering and the differentiation operations is related to theuse of a variable filter, which in the more general config-uration is anisotropic. In this case the filter width is thevector δ ( x ) = ( δ ( x ) , δ ( x ) , δ ( x )).For reading convenience and as the isotropic stretchingconfiguration is conceptually non reductive, a scalar filterscale δ ( x ) is assumed in what follows. However, the gen-eral anisotropic configuration of stretching is dealt within the Appendix, which also specifies a filtering formula-tion that is suitable for wall-bounded flows.Let us suppose we have chosen a given class of integra-tion volumes V δ = (cid:8) η ∈ IR : k η k h δ (cid:9) and an average operation for the variable f ( x ): h f i δ = 1 V δ Z V δ f ( x + η ) d η = 1 V Z V f ( x + δ ξ ) d ξ , (1)where the transformation ξ = η /δ has been used and, asa consequence, V = V δ /δ . Please note that, with thischoice, the width of the averaging volumes is twice thefilter scale.A variable filter scale is introduced by allowing δ to bea function of point, δ = δ ( x ). In this case ∂∂x i h f i δ = ∂∂x i [ 1 V Z V f ( x + δ ( x ) ξ ) d ξ ]= h∇ f · δ i i δ + ∂δ∂x i ( x ) h∇ f · ξ i (2)By virtue of the fact that ∂∂δ h f i δ = ∂∂δ (cid:20) V Z V f ( x + δ ( x ) ξ ) d ξ (cid:21) = 1 V Z V ∇ f · ξ d ξ , (3)recalling that h∇ f · δ i i δ = h ∂f∂x i i δ , it results that the filterof the derivative is a differential operator acting on thefiltered field: h ∂f∂x i i δ = ∂∂x i h f i δ − ∂δ∂x i ∂∂δ h f i δ (4)The non commutation term C ′ i , which is defined as C ′ i ( h f i δ ) = h ∂f∂x i i δ − ∂∂x i h f i δ , (5)can be represented through (4) by the product of the filterspace derivative and the filter derivative of the filteredvariable: C ′ i ( h f i δ ) = − ∂δ∂x i ∂∂δ h f i δ (6)The here proposed method is based on an approxima-tion of relation (6). The problem could be faced adoptinga truncated series expansion of h f i δ in terms of powersof δ (Ghosal and Moin, 1995)[1]. However, this wouldincrease the order of the equations, and thus require ad-ditional boundary conditions. Here a numerical approx-imation of the δ first derivative is used in conjunctionwith truncated δ expansions. Let us write the secondorder finite difference approximation ∂ h f i δ ∂δ = 12 h ( h f i δ + h − h f i δ − h ) + O ( h ) (7)Choosing h = δ∂ h f i δ ∂δ = 12 δ ( h f i δ − h f i ) + O ( δ ) (8)Now, the problem to face is that of the approximation of h f i = f and h f i δ in terms of relevant averaged quan-tities. Using a Taylor expansion of the integrating func-tion in (1), we obtain the following expression for h f i δ ,in terms of f , and the filter width: h f i δ ( x ) = f ( x ) + 12 a , , ∇ f ( x ) δ ++ 14! (cid:2) a , , ( ∂ + ∂ + ∂ ) f ( x )+ 6 a , , ( ∂ ∂ + ∂ ∂ + ∂ ∂ ) f ( x ) (cid:3) δ + O ( δ )= f ( x ) + F [ f ] δ + F [ f ] δ + O ( δ ) (9)where coefficients a ijk are defined as a ijk = 1 V Z V ξ i ξ j ξ k d ξ . (10)and the operators F , F as: F [ · ] = 12 a , , ∇ · (11) F [ · ] = 14! (cid:2) a , , ( ∂ + ∂ + ∂ ) · +6 a , , ( ∂ ∂ + ∂ ∂ + ∂ ∂ ) · (cid:3) (12) From (9) it follows that f = h f i δ − F [ f ] δ + O ( δ ) (13)and then, averaging (13) on a volume of linear dimension2 δ , h f i δ = hh f i δ i δ − h F [ f ] δ i δ + O ( δ ) . (14)However, from (9) it can be observed that h F [ f ] δ i δ = F [ f ] δ + 4 F [ F [ f ] δ ] δ + ... == δ F [ f ] + O ( δ )so that h f i δ = hh f i δ i δ − F [ f ] δ + O ( δ ) (15)When the expressions (13) for f and (15) for h f i δ areintroduced into (8), one obtains ∂ h f i δ ∂δ = 12 δ ( hh f i δ i δ − h f i δ ) + O ( δ ) (16)When using (16), the non commutation term C ′ i [see (6)]can be approximated by˜ C ′ i ( h f i δ ) = − ∂δ∂x i δ ( hh f i δ i δ − h f i δ ) (17)which implies C ′ i ( h f i δ ) = ˜ C ′ i ( h f i δ ) + ∂δ∂x i O (cid:0) δ (cid:1) . (18)In order to analize the approximation error ( ∂ i δ ) O ( δ )and give a true estimate of it, let us write: δ ( x ) = ∆ ϕ ( x ) , (19)where ∆ is a reference value of the filter width whichis usually associated to the portion of the domain whereconditions of near homogeneity of the flow hold. Function ϕ ( x ) is a positive non dimensional function which belongsto the interval [ δ min / ∆ , ϕ ( x ) is constant and equal to 1. Function ϕ ( x ) varies in the inhomogeneous regions, though it keepsvalues that are greater than δ min / ∆, which is a value thatmust correspond: (i) – when the local scale invariancemay be supposed – to a convenient minimum value of thefilter width still inside the inertial range, and (ii) – whenthe local scale invariance does not hold – to a scale of theorder of the scale which characterises the local turbulencestructure, as, in case of wall flows, is the scale of thequasi-streamwise vortices peculiar to the viscous sublayer(see Moin and Kim, 1982[14], and Ghosal, 1999[5]).Introducing (19) and ∂ δ = ∆ − ∂ ϕ into (6) and deduc-ing the δ derivative from expansion (9), written up to thefourth order of accuracy, the non commutation term canbe estimated C ′ i ( h f i δ ) = − ∂ϕ∂x i ∂∂ϕ h f i ∆ ϕ == − ∂ϕ∂x i ϕ ( x ) (cid:0) F [ f ]∆ +2 ϕ ( x ) F [ f ]∆ + O (∆ ) (cid:1) (20)Using expantion (9) twice, after having inserted (19), ap-proximation (17) can be estimated as˜ C ′ i ( h f i δ ) = − ∂ϕ∂x i ϕ ( x ) (cid:0) F [ f ]∆ + (4 ϕ ( x ) F [ f ]+ F [ ϕ F [ f ]])∆ + O (∆ ) (cid:1) (21)To keep the validity of estimates (20) and (21), care mustbe taken to select a function ϕ ( x ) which, in the region offilter variation, also posseses first and second derivativesof O (1). Possible examples are trascendental functionssuch as arctan( x ) , tanh( x ).A comparison of (20) and (21) yields C ′ i ( h f i δ ) − ˜ C ′ i ( h f i δ ) ≈ O (∆ ) . (22)Consequently, when (17) is used, a fourth order noncommutation term in (4) is produced instead of the sec-ond order error, which would be obtained by totally ne-glecting the lack of commutation [1]. Introducing a finitedifference approximation of a higher order than (7), and,consequently, further levels of average, it could be possi-ble to increase the accuracy of the approximation of thenon commutation term, leading to a higher order errorin (22).This analysis pertains to differential operators of thefirst order. The analysis is similar for the second orderdifferential operators. The structure of the correctionterms remains the same and to reach the fourth order ofaccuracy the same number of levels of average must bemantained. The approximation of the non commutationterm now includes the filter of the variable spatial firstderivatives.The non commutation term of the second derivatives,being defined by C ′′ ii ( h f i δ ) = h ∂ f∂x i i − ∂ ∂x i h f i δ , (23)can be obtained by taking the derivative of the firstderivative (2) as C ′′ ii ( h f i δ ) = − ∂ δ∂x i V Z V X j =1 ξ j ∂f∂x j ( x + δ ξ )d ξ − ∂δ∂x i V Z V X j =1 ξ j ∂ f∂x j ∂x i ( x + δ ξ )d ξ − (cid:18) ∂δ∂x i (cid:19) V Z V X j,k =1 ξ j ξ k ∂ f∂x j ∂x k ( x + δ ξ )d ξ (24)that is, C ′′ ii ( h f i δ ) = − ∂ δ∂x i ∂∂δ h f i δ − ∂δ∂x i (cid:18) ∂ ∂δ∂x i h f i δ (cid:19) − (cid:18) ∂δ∂x i (cid:19) ∂ ∂δ h f i δ (25) and can consequently be approximated using the finitedifference for the δ -derivatives, as performed for the noncommutation term of the first derivatives. The use of thestandard three-point formula for the second derivative of h f i δ with respect to δ and of relation (16) for ∂ δx i h f i δ [second term on the right hand side of equation (25)],yields˜ C ′′ ii ( h f i δ ) = − ∂ i δ ( x ) δ ( x ) [ h ∂ i h f i δ i δ − ∂ i h f i δ ] − ( ∂ i δ ( x )) + δ ( x ) ∂ i δ ( x )2 δ ( x ) [ hh f i δ i δ − h f i δ ] , (26)whose order of accuracy can be determined by again us-ing the trasformation δ ( x ) = ϕ ( x )∆. If the expansionsof the non commutation terms C ′′ is compared with itsapproximantion ˜ C ′′ , it is seen that, also in this case, theerror is O (∆ ).Even if the present analysis is based on the use ofthe volume averages, it can be observed that it remainsvalid in the case where a more general kind of filtering isadopted. A weight function g ( ξ ), introduced in (1), onlymodifies the coefficients a ijk [Eq. (10)], which shouldnow be defined as a ijk = Z g ( ξ ) ξ i ξ j ξ k d ξ , However, the non commutation term (6) and its approx-imation remain unchanged, provided the weight functionhas a compact support. It is always possible to choose δ ( x ), so that the actual integration domain of the filterlies inside the flow domain (for instance setting a valuethat is lower than the distance from the wall of the firstlayer of grid points for the minimum of δ ( x )). In thisway, the compactness of the support prevents the errorlinked to the presence of finite boundaries[3]. It should benoted that this procedure operates in the physical spaceand does not rely on the use of a mapping function ofthe non uniform grid. Centered volumes of average havebeen adopted [see Eq. (1)], even though they are notstrictly necessary as far as the average process is con-sidered. However, this choice shows two advantages: (i)physically, when the flow is incompressible, the center ofthe volume of average is also the center of gravity andthus the point of application of the average momentum,(ii) analytically, it allows for a compact, and second or-der in ∆, representation of the non commutations terms,which in turn are approximated by the present procedurewith an accuracy of the fourth order. The choice of noncentered volume of averages, which, in principle, is math-ematically feasible, yields a first order in ∆ representa-tion of the non commutation terms, which would also leadto a much more cumbersome analytical structure[1, 6]. A. Non commutation terms in the averagedincompressible Navier-Stokes equations
Let us consider the incompressible Navier-Stokes equa-tions written in the form: ∂ i u i = 0 (27) ∂ t u i + ∂ j ( u i u j ) + ∂ i p − ν ∂ jj u i = 0 (28)If a filter operator is applied the system becomes h ∂ i u i i δ = 0 (29) h ∂ t u i i δ + h ∂ j ( u i u j ) i δ + h ∂ i p i δ − ν h ∂ jj u i i δ = 0 . (30)When δ , the linear scale of filtering (see for instance thedefinition proposed in Sec. II) is not uniform in the flowdomain, the averaging and differentiation operations nolonger commute. By introducing the subgrid turbulentstresses R ( δ ) ij = h u i i δ h u j i δ − h u i u j i δ and the non commu-tation terms C i ′ , C ii ′′ , for the first and second derivatives,as discussed in Sec. II for the isotropic filter configura-tion and in the Appendix for the general anisotropic andthe wall-bounded flow configuration [see the isotropic re-lations (5), (6), (23) and (25); the anisotropic relations(A.5)–(A.6), (A.17), (A.10), (A.19); the wall anisotropyrelations (A.16), (A.20)), the averaged equations arewritten as ∂ i h u i i δ = −C ′ i ( h u i i δ ) (31) ∂ t h u i i δ + ∂ j ( h u i i δ h u j i δ ) + ∂ i h p i δ − ν ∂ jj h u i i δ − ∂ j R ( δ ) ij == −C ′ j ( h u i i δ h u j i δ ) − C ′ i ( h p i δ ) + ν C ′′ jj ( h u i i δ ) + C ′ j ( R ( δ ) ij ) . (32)Together with what has been explainded in detail in Sec.II (see (17) and (26)) the present procedure approximatesthe non commutation term terms on the right hand sideof (31, 32) with an accuracy of the fourth order. Thecorrection terms can thus be represented by the follow-ing group of relations which are determined from the fieldinformation obtained through two successive average lev-els (the second being computed over a linear scale 2 δ ):˜ C ′ i ( h u i i δ ) = − ∂ i δ ( x )2 δ ( x ) [ hh u i i δ i δ − h u i i δ ] (33)˜ C ′ j ( h u i i δ h u j i δ ) = − ∂ j δ ( x )2 δ ( x ) [ hh u i i δ h u j i δ i δ − h u i i δ h u j i δ ](34)˜ C ′ i ( h p i δ ) = − ∂ i δ ( x )2 δ ( x ) [ hh p i δ i δ − h p i δ ] (35)˜ C ′′ jj ( h u i i δ ) = − X j =1 (cid:26) ∂ j δ ( x )2 δ ( x ) [ h ∂ j h u i i δ i δ − ∂ j h u i i δ ] ++ ( ∂ j δ ( x )) + δ ( x ) ∂ j δ ( x )2 δ ( x ) [ hh u i i δ i δ − h u i i δ ] ) (36) ˜ C ′ j ( R ( δ ) ij ) = − ∂ j δ ( x )2 δ ( x ) h h R ( δ ) ij i δ − R ( δ ) ij i , (37)which must be accordingly modified in the case ofanisotropy of the stretching of the computational grid,see the previous comments and the Appendix [(A.10)–(A.11), (A.15)–(A.16) and (A.19)–(A.20)].The adoption of the volume average allows the fil-tered variables to be fully supported inside the physi-cal domain. As a consequence, a peculiar property ofthe present procedure is that there is an absence of noncommutation terms associated to a finite or semi-infinitecomputational domain in the filtered equations. Suchterms arise when the filtering operator requires the ex-tension of the dependent variables beyond the rim of thedomain [see Fureby and Tabor (1997)[3]].It can be observed that the use of the double level ofaverage highlights the convenience of coupling this pro-cedure to subgrid models which also employ it. Thesesubgrid models are the mixed model, by Bardina et al. (1980)[9] and, in general, all models that apply the dy-namical procedure, by Germano et al. (1991)[7] and Ger-mano (1992)[8].The filtered equations (29, 30) are invariant underGalileian transformations. Under transformation t ′ → t, x ′ → x + ct , a spatial variation of the filter scale δ ( x ),in the x, t reference system, becomes a spatio-temporalvariation δ ( x ) = δ ( x ′ − ct ′ ) = δ ( x ′ , t ′ ) in the x ′ , t ′ system.The temporal dependence of the filter scale yields to thepresence of a non commutation term which is also asso-ciated to the non stationary term. The trasformed noncommutation terms, and the relevant approximations, re-leased by the unsteady and the convective terms, howevercancel each other. Furthermore, the terms obtained fromthe divergence of the stress tensor in a system like (31,32) and the corresponding approximations (35) - (37) areall Galileian invariants, which assures that the variablescale filtered equations and their approximations are alsosuch.Another general implication, linked to the presence ofa finite domain, is that the boundary conditions for thefiltered variables should be different from those for theunfiltered variables. The problem of wall boundary con-ditions for the filtered field could be treated with this pro-cedure by adopting one of the classical approximated con-ditions, which rely on the introduction of a special sub-grid model, that is, the wall model, which is apt to rep-resent the inner layer dynamics and which puts the firstgrid point inside the logarithmic layer[15] , [16] , [17] , [18].It could also be treated by placing grid points well in-side the viscous sublayer to resolve the near-wall dynam-ics and by assuming no slip and impermeability bound-ary conditions. It should be recalled that the latterconditions, which in theory should not be used for fil-tered velocities, introduce an error of O (∆ ), indepen-dently of the filter shape. Their use requires the sub-grid model, which should represent the non homogeneousand anisotropic structure of the viscous and buffer lay-ers, to be altered in the inner region. The employementof anisotropic models based on a tensorial turbulent vis-cosity would be opportune, see Horiuti (1990)[19], Caratiand Cabot (1996)[20], the review monography by Sagaut(2001, Chap. 5.3)[6] and also the differential angular mo-mentum model (Iovieno and Tordella, 2002)[21], which,being based on the representation of the turbulent vis-cosity through the moment of momentum vector, is wellsuited to assume an anisotropic formulation. III. NUMERICAL TESTSA. A priori tests on the turbulent channel flow
In this section a set of a priori tests is presented, whichprovides information on the field distribution of the noncommutation terms, their relevant approximations andtheir ratios with respect to the physical terms from whichthey arise. The data correlating the approximated andexact non commutation terms have been determined byfiltering the direct numerical simulation of the turbulentplane channel flow at Re τ = 180, as performed by Al-fonsi et al. (1998)[11] and Passoni et al. (1999)[12], andat Re τ = 590, as performed by Moser et al. (1999)[13].The longitudinal momentum balance, which implies azero pressure non commutation term, is considered. Rep-etition on two different grid levels has been performed.All the data here presented have been averaged over aninterval of 1.2 revolution times.Figure 1 shows the distributions of the exact non com-mutation terms and of their approximation according tothe present procedure, at two filtering levels. The convec-tion term of the longitudinal momentum balance, its gridand subgrid scale decomposition and the diffusive termare shown in parts (a), (b) and (c), respectively. Thefilter δ ( x ) = (∆ x, ϕ ( y )∆ y, ∆ z ), with constants ∆ x, ∆ y and ∆ z , varies along the transversal non dimensional y direction according to ϕ ( y ), where ϕ ( y ) ∈ [0 , y ∈ [ − ,
1] is a function of at least class C . The vari-ation of ϕ ( y ) has been laterally arranged (in 20% of thechannel width along the walls) as follows ϕ ( y ) = tanh a ( y + 1) tanh a (1 − y )tanh a , a = 4 , (38)where a is the parameter that controls the gradient ofthe filter scale at the wall. The non commutation termson the first and second derivatives have been determinedfor such an anisotropic structure of the filter through theuse of relations (A.16) and (A.20). The data in Fig. 1 (a,c) show that the present procedure yields, on average, forthe Reynolds stress | C ′ | ∆ | C ′ | = 0 . , | C ′ − ˜ C ′ | ∆ | C ′ − ˜ C ′ | = 0 . , y + -0.003-0.002-0.00100.0010.0020.003 C y ( < u v > ) , C y ( < u v > ) ~ ∆ = 9.84 10 -2 , Re τ =180∆ = 4.92 10 -2 , Re τ =180 .. ∆ = 4.92 10 -2 , Re τ =590 . ++++++ y + -0.003-0.002-0.0010.0000.0010.0020.003 C ’ y ( < u >< v > ) , - C ’ y ( R xy ) ~~ ∆ = 9.84 10 -2 ∆ = 4.92 10 -2 .. +++++ y + -0.006-0.004-0.00200.0020.0040.006 C " y y ( < u > ) , C " y y ( < u > ) ~ ∆ = 9.84 10 -2 ∆ = 4.92 10 -2 . . ++++++ FIG. 1. (a): ( —– ) Average values of the exact non commu-tation convection term and ( – – – ) correspondent averagevalues predicted by the procedure. (b): Noncommutation ap-proximated terms for the resolved Reynolds stresses ( – . –. – ) and for the subgrid scale stresses ( – – – ). (c): ( —– )Average values of the noncommutation diffusive term and ( –– – ) correspondent values predicted by the procedure. | C ′ − ˜ C ′ | ∆ | C ′ | ∆ = 0 . , | C ′ − ˜ C ′ | | C ′ | = 0 . , and for the viscous stress | C ′′ | ∆ | C ′′ | = 0 . , | C ′′ − ˜ C ′′ | ∆ | C ′′ − ˜ C ′′ | = 0 . , | C ′′ − ˜ C ′′ | ∆ | C ′′ | ∆ = 0 . , | C ′′ − ˜ C ′′ | | C ′′ | = 0 . , where ∆ = 4 . · − .These figures have been obtained by neglecting thedata that belong to the first 5% near the wall, where thenumerical uncertainty due to the spatial discretization ishigh - especially as regard the exact non commutationterms computations - and where these results deteriorateby nearly 25%.The results of the numerical test at Re τ = 590 as com-pared with those at Re τ = 180, see Fig. 1 (a), show agood invariance of the procedure accuracy with respectto variations of the flow control parameter.It should be remarked, that a procedure capable topredict at worst, using the rather large value of δ =4 . · − , the 90% of the value of the non commutationterms may be considered accurate. In fact, the analysis ofnumerical errors in LES of Turbulence with cutoff in theinertial range (errors due to spatial discretization: finite-differencing errors and aliasing errors) shows that the re-sulting errors are very large, of the same order and evenlarger than the magnitude of the subgrid term over mostof the wavenumber interval, for finite difference schemesup to eighth-order accurate, irrespective of the grid reso-lution (cf. Ghosal S., 1996, pp. 201-202)[22]. In such ageneral situation, it may be considered it a success thatthe procedure is capable of predicting nine-tenth of thevalue of the non commutation terms. The relevant aver-age error cannot spoil the overall numerical reliability ofthe simulations since it is about one order of magnitudelower than the errors due to the spatial discretation.A warning is necessary regarding the numerical compu-tation of the exact values of the non commutation terms.Direct computation through the definition is not recom-mended, as, even with the implementation of a numericaldifferentiation of the sixth order of accuracy, it artificiallyamplifies the fluctuations that are naturally present inthe data field. The exact values of the non commutationterms should be correctly evaluated by using the integralrepresentation of the derivatives ∂/∂δ or ∂/∂δ j , such as(3) in Sec. II or (A.3) in the Appendix. No such numeri-cal problems affect the computation of the approximatednon commutation terms.Figure 2 provides information on the relative impor-tance of the exact non commutation terms with respectto the physical terms which causes them. In part (a)the average ratio |C ′ y ( h uv i ) /∂ y ( h uv i ) | has been plotted for∆ = 4 . · − , ∆ = 9 . · − and also with an increaseof the wall value of the stretching factor a = ∂ y ϕ = 8. Itcan be seen, that close to the point where the Reynoldsstress reaches its maximum ( y + ≈ .
38 are reached with the coarsergrid. The comparison of these results with the resultsof Fig. 1(b) indicates that the average exact value of theconvection non commutation term is of the same orderthan the average divergence of the subgrid stresses. Fur- thermore, it is interesting to observe that the doublingof the stretching factor a increases the relevance of thenoncommutation terms in this region nearly as much asthe doubling of the grid coarsening does, see in Fig. 2(a)the near wall region where y + ≤
14. Figure 2(a) alsoshows a positive comparison of the field integral value of |C ′ y ( h uv i ) /∂ j ( h uu j i ) | given by Fureby and Tabor[3] withthe distribution of the same ratio that has been yieldedby the database used here.[11] − [12]With respect to the average ratio |C ′′ yy ( h u i ) /∂ yy ( h u i ) | ,Fig. 2(b) yields maxima local values of about 100% forthe coarser resolution, and about 60% for the finer one,close to where the relevant non commutation terms reachtheir local maxima near to the wall. Leaving aside localmaximum values detached from the wall and relevant tothe coarser resolution, this ratio, in the central part ofthe field, settles to constant lower values close to 0 . ± . B. Commutation error on an analytical solution
As seen in Sec. II A, there are four type of non commu-tation terms in the variable scale filtered incompressibleNavier-Stokes equations (31), (32), which are all sourceterms. If one limits the analysis of their influence on theflow solution to the determination of the field distribu-tion of the values they take with regards to the valuestaken by the original terms of the equations, one wouldmainly find the foreseen result that the non commutationterms are not negligible where the gradient of the filterscale is high. This is however not sufficient to under-stand the way, localized rather than extended, in whichthe commutation error affects the flow solution.For this purpose, it has here been considered useful tostudy the behaviour of an extremely simple flow model, asort of conceptual model, which has two characteristics:– just one type of commutation source term is present, –its exact filtered solution, that is, the variable scale fil-tered solution not affected by the commutation error, isknown and thus could be used as the reference solution.The second characteristics can only be obtained by fil-tering the exact solution of the unfiltered equation of themotion.On the other hand, to prove the efficiency of any givenprocedure for the correction of the commutation error, itis also necessary to know the exact filtered solution of atest flow, which, in turn, requires the knowledge of theexact flow unfiltered solution.Such a reference state cannot be found in a turbu-lent configuration of flow, for which no exact solution isavailable. Reference is therefore made to a laminar flow,which has an exact solution. One should note that, insuch a case, the filtered equation of the motion, whenthe filter length is a function of the point but the com-mutation error correction is not considered, is identicalto the unfiltered equation.The steady laminar incompressible channel flow hasbeen selected as the test flow, since it has only one com- -1.0 -0.90.00.10.20.30.40.5 C y ( < u v > ) / ∂ y < u v > , C y ( < u v > ) / ∂ j < uu j > y ∆ = 4.92 10 -2 , a = 4 ∆ = 9.84 10 -2 , a = 4 ∆ = 4.92 10 -2 , a = 8Fureby and Tabor, Ref.[3]
18 54 90 126 162 y + -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 -0.0 y C " yy ( < u > ) / ∇ ( < u > ) ∆ = 4.92 10 -2 , a = 4 ∆ = 9.84 10 -2 , a = 4 ∆ = 4.92 10 -2 , a = 8 y + FIG. 2. Average of the absolute value of the ratio betweenthe exact non commutation terms and the physical correspon-dent terms in the motion equation. (a) Convection, ( —– ) |C ′ y ( h uv i ) /∂ j ( h uu j i ) | , (b): diffusion. Adimensionalization bymeans of channel semi-width and u τ . mutation term, the diffusion one, see (36), and its solu-tion is analytically known.The non dimensional momentum equation for thesteady incompressible channel flow is written as ∂ y u = Re ∂ x p = − ( Gρ ) / d / µ , ∀ y ∈ (0 ,
1) (39)where G is the modulus of the dimensional longitudinalpressure gradient and the adimensionalization is basedon the channel width d , the gradient G and the density ρ (the non dimensional pressure gradient ∂ x p results equalto − u (0) = 0 , u (1) = 0 (40)The corresponding filtered equation is ∂ y h u i δ ( y ) + Re = −C ′′ ( h u i δ ) , ∀ y ∈ (0 ,
1) (41)where C ′′ ( h u i δ ) is set to zero to determine the solu-tion which neglects the commutation error and where y ϕ ( y ) , E r / E r , m a x -0.06-0.04-0.020 0.020.040.060.08 C " ( < u > ) / < u > " C" / " E r / E r,max ϕ (y) FIG. 3. ( —– ) Filter scale across the channel ( δ = ϕ ( y )∆, a = 8, see (44), III B); ( - - - - ) local values of the com-mutation term ˜ C ′′ , see (42), referring to the diffusion term;( —– - —– ) R = E r /E r max , relative commutation er-ror for the solution of the filtered non corrected equation( C ′′ = 0) normalized with respect to the field peak value, E r = [( h u i − h u i exact ) / h u i exact ]. y E a , E r B A E rr o r R ec o ve r y FIG. 4. Absolute E a = (8 /Re )( h u i−h u i exact ) ( —– ) and rela-tive E r = [( h u i−h u i exact ) / h u i exact ] ( - - - - ) error distributionsof the filtered velocity distributions: A – without the com-mutation correction, B – with the commutation correction.Curve ( —– - —– ): error recovery ( E a ( A ) − E a ( B )) /E a ( A )with the distance from the wall. The exact filtered velocityreference distribution is h u i exact = ( Re/ y (1 − y ) − δ ( y )]. C ′′ ( h u i δ ) is approximated by˜ C ′′ ( h u i δ ) = − δ ′ ( y ) δ ( y ) [ hh u i ′ δ i δ − h u i ′ δ ] − δ ′ ( y ) + δ ( y ) δ ′′ ( y )2 δ ( y ) [ hh u i δ i δ − h u i δ ](42)to determine the solution which accounts for the commu-tation error with the present procedure. The numericalsolution of equation (41) is determined by solving thecorresponding unsteady filtered equation ∂ t h u i δ − Re ∂ y h u i δ = 1+ 1 Re ˜ C ′′ ( h u i δ ) , ∀ y ∈ (0 ,
1) (43)through a fourth order Runge-Kutta time integrationscheme – carried out until the steady state is reached– coupled to a fourth order finite-difference discretiza-tion of the domain. The double level of average has beencomputed using a third order Hermitian quadrature for-mula.The filter of the exact solution gives the velocity dis-tribution h u i exact = Re y (1 − y ) − δ ( y )], which con-stitutes the reference on which the constrast betweenthe commutation corrected filtered solution and the noncorrected filtered solution is based. The filter scalevaries along the transversal direction according to δ ( y ) = ϕ ( y )∆, where ϕ ( y ) ∈ [0 , .
1, is a functionof at least class C . The variation of ϕ ( y ) ϕ ( y ) = tanh 2 ay tanh 2 a (1 − y )tanh a , (44)has been laterally arranged in 20% of the channel widthalong the wall, setting the parameter a = 4, see Fig. 3.By contrasting the filtered exact solution, the absoluteand relative errors relevant to the corrected and non cor-rected solutions are compared in Fig. 4. The correspond-ing error recovery, with the wall distance, is also shown inFigs. 2 and 3. It can be seen that the present proceduregreatly reduces the commutation errors: in the centralpart of the flow, an almost full recovery is obtained.For a comparison of the distributions of the local valueof the non commutation term and of the relative commu-tation error on the solution, see again Fig. 3. It is im-portant to observe that, by neglecting the commutationcorrection, the field results to be affected by a system-atic error not only in the region where the filter lengthvaries, but also in the region where it is constant. Thisbehaviour is due to the accumulation of errors on the ve-locity variable and its derivative, which is due to the lackof the two diffusion addenda (see (41), (42)) that shouldenter the momentum balance equation. Even throughthese terms are significantly different from zero in a lim-ited portion of the flow they affect the entire field to agreat extent. In the central part of the flow, where thenon commutation term is very small (Fig. 3), the relativecommutation error results to be of the same order as thepeak value of the field, while the absolute error reachesits maximum value (see again Fig. 4). IV. CONCLUSIONS
A procedure to explicitly insert the correction terms inorder to counteract the commutation error associated tothe use of a variable filter scale in the filtered equations ofmotion is here presented. With this procedure it is possi-ble to directly compensate for the commutation error on the filtered field. The procedure uses volume average fil-tering, but more general filter operators are also possible.Both isotropic and fully anisotropic filtering configura-tions are considered. Approximated commutation terms,with an accuracy of the fourth order in the filter width,are inserted into the motion equations, which do not in-crease their differential order. The difficulties related tothe addition of further boundary conditions are thereforeavoided. The proposed representation of the commutatoroperators is based on truncated expansions in the filterwidth of finite difference approximations, that make useof a multilevel average operation. This fact suggests thejoint use of the present procedure with subgrid modelswhich need an explicit filtering of the equations of mo-tion, such as dynamic and mixed models.A set of a priori tests, with a plane channel flow DNS( Re τ = 180) as a test field, proves the good correla-tion that the present procedure yields between the ap-proximate and ”exact” non commutation terms. It alsoprovide information on the relative importance of theseterms with respect to the original physical terms. The in-fluence of the field resolution on the general non commu-tation term is confirmed to be O (∆ ). Asymptotically,the accuracy of the present approximation is expectedto be O (∆ ). At the resolution levels corresponding to Re τ = 180, these tests show a reduction of the abso-lute errors, after halfing the reference filter scale, that isnearly O (∆ ).The filtering is a mathematical operation, that is notspecific of the Navier-Stokes equations and is indepen-dent of the solution typology. When it is varied, the filtercauses one kind of non commutation term for each dif-ferential term present in the equations of the motion. Atest, for which the analytically exact (unfiltered and as aconsequence filtered) solution is available, has been con-sidered to overcome the limitation of an analysis, basedon the a posteriori determination of the relative orderof magnitude of the non commutation terms, with re-gard to the original physical terms of the motion equa-tions, and to analyze the effects of the commutation er-ror on a flow solution. The chosen test flow is the twodimensional incompressible laminar channel flow, whosedynamics consist of the balance between the constantlongitudinal pressure force and the viscous diffusion. Inthis case, only one type of commutation source term –diffusion – is present, and this is only of relative impor-tance in the lateral part of the flow, according to thefilter gradient dependence. The error, due to the lack ofa non commutation term in the motion equation, is how-ever also transferred to the central part of the flow. Theresult is a biased filtered velocity distribution where therelative error in the central part of the flow, where thefilter gradient is zero, is of the same order of magnitudeas the local maximum error of the field, which is situatedat a distance from the wall of about 15% of the channelwidth. It has been shown that the present procedure canreduce the commutation error by one order of magnitudein the central part of the field.0 ACKNOWLEDGMENTS
The authors would like to thank Professor G.Passonifor helpful discussions and for making the turbulent chan-nel flow database available.
APPENDIX: NON COMMUTATIONAPPROXIMATION FOR ANISOTROPICFILTERS
When the geometry of the flow domain requiresstretching each direction independently, which implies δ ( x ) = ( δ ( x ) , δ ( x ) , δ ( x )), it is opportune to adopt aclass of integration volumes of the kind V δ = (cid:26) η ∈ IR : k ( η δ , η δ , η δ ) k h (cid:27) (A.1)and an average operation for the variable f ( x ) = f ( x j + δ j ξ j ): h f i δ = 1 V δ Z V δ f ( x + η ) d η = 1 V Z V f ( x j + δ j ξ j ) d ξ , (A.2) where = (1 , , η j = δ j ξ j (withdet ( ∂η i /∂ξ k ) = δ δ δ ) has been introduced and thus V = V δ /δ δ δ .In such a situation, by virtue of the fact that ∂∂δ k h f i δ = ∂∂δ k " V Z V f ( x j + δ j ( x ) ξ j ) d ξ = X j =1 V Z V ∂f∂x j ( x j + δ j ( x ) ξ j ) ×× ∂∂δ k ( x j + δ j ( x ) ξ j ) d ξ = 1 V Z V ξ k ∂f∂x k d ξ , (A.3)it is obtained ∂∂x i h f i δ = ∂∂x i [ 1 V Z V f ( x j + δ j ( x ) ξ j ) d ξ ] = X j =1 V Z V ∂f∂x j ( x j + δ j ( x ) ξ j ) ∂∂x i ( x j + δ j ( x ) ξ j ) d ξ = X j =1 V Z V [ δ K ji ∂f∂x j + ∂δ j ∂x i ξ j ∂f∂x j ] d ξ = 1 V Z V ∂f∂x i d ξ + X j =1 ∂δ j ∂x i V Z V ∂f∂x j ξ j d ξ = h ∂f∂x i i δ + X j =1 ∂δ j ∂x i ∂ h f i δ ∂δ j , where δ K ji is the Kronecker unit tensor. As in Sec. II[see relation (4)] for the isotropic configuration, it resultsthat the anisotropic filter of the derivative is a differentialoperator that acts on the filtered field: h ∂f∂x i i δ = ∂∂x i h f i δ − X j =1 ∂δ j ∂x i ∂∂δ j h f i δ (A.4)The anisotropic non commutation term C ′ i , which is de-fined as C ′ i ( h f i δ ) = h ∂f∂x i i δ − ∂∂x i h f i δ , (A.5)can now be represented through (A.4) as the sum ofthe products of the filter space derivatives and the fil-ter derivatives of the filtered variable: C ′ i ( h f i δ ) = − X j =1 ∂δ j ∂x i ∂∂δ j h f i δ (A.6) Proceeding in strict analogy with what has been donein Sec. II, the anisotropic non commutation term can beapproximated through second order centered finite dif-ferences:˜ C ′ i ( h f i δ ) = − X j =1 ∂δ j ∂x i δ j (cid:16) h f i δ + δ j e j − h f i δ − δ j e j (cid:17) + O ( δ j )(A.7)Again, using a Taylor expansion of the integrating func-tion in (A.2), where only even derivatives appear sincethe domain of integration (A.1) is symmetric with re-spect to all the integration variables, expressions areobtained for quantities such as h f i δ and h f i δ ± δ j e j , interms of f , and the filter width – which is now definedas δ j = ∆ ϕ j ( x ) , ∀ j , where ∆ is a reference value for | δ | and 0 ≤ ϕ j ( x ) ≤ , ∀ j, x : f ( x ) = h f i δ ( x ) −
12 ∆ X j =1 a j ϕ j ( x ) ∂ f∂x j ( x )1+ O (∆ ) , (A.8)where a j = 1 V Z V ξ j d ξ and h f i δ ± δ j e j = hh f i δ ( x ) i δ ± δ j e j − ∆ X j =1 a j h ϕ j ( x ) ∂ f∂x j i δ ± δ j e j + O (∆ )= hh f i δ ( x ) i δ ± δ j e j − ∆ X j =1 a j ϕ j ( x ) ∂ f∂x j + O (∆ ) , (A.9)since from (A.8) h g i δ = g + O (∆ ) , ∀ g .The basic approximation for the anisotropic non com-mutation term C ′ i is derived through these expansions as˜ C ′ i ( h f i δ ) = − X j =1 ∂δ j ∂x i δ j (cid:16) hh f i δ i δ + δ j e j − hh f i δ i δ − δ j e j (cid:17) (A.10)= − X j =1 ∂ϕ j ∂x i ϕ j (cid:16) hh f i δ i δ + δ j e j − hh f i δ i δ − δ j e j (cid:17) , (A.11)while the accuracy of the anisotropic first derivative com-mutation error and its approximation can be verified tobe C ′ i ( h f i δ ) = ∆ X j =1 a j ϕ j ∂ϕ j ∂x i ∂ f∂x i + O (∆ ) (A.12)˜ C ′ i ( h f i δ ) = C ′ i ( h f i δ ) + O (∆ ) . (A.13)By filtering only in the j direction, expansion (A.8) be-comes hh f i δ i δ j e j i = f ( x ) + 2∆ ϕ j ( x ) ∂ f∂x j ( x ) + O (∆ ) , which implies hh f i δ i δ j e j −h f i δ = hh f i δ i δ + δ j e j −hh f i δ i δ − δ j e j + O (∆ ) , and, though keeping the same order of accuracy, theequivalent representation for both relations (A.10, A.11) becomes˜ C ′ i ( h f i δ ) = − X j =1 ∂δ j ∂x i δ j (cid:16) hh f i δ i δ j e j − h f i δ (cid:17) (A.14)= − X j =1 ∂ϕ j ∂x i ϕ j (cid:16) hh f i δ i δ j e j − h f i δ (cid:17) . (A.15)For flow fields where the domain grid needs to bestretched along only one direction, say y and whose typ-ical examples are two dimensional wall-bounded flows,the last representation yields very simple approxima-tion formulae. By adopting the widely used notation δ ( x ) = (∆ x, ϕ ( y )∆ y, ∆ z ) in such a case with constants∆ x, ∆ y and ∆ z , the anisotropic approximation for thefirst derivative non commutation term results to be˜ C ′ y ( h f i δ ) = − ∂ϕ∂y ϕ ( y ) (cid:16) hh f i δ i ϕ ( y )∆ y − h f i δ (cid:17) . (A.16)The anisotropic non commutation term of the secondderivatives, being defined by C ′′ ii ( h f i δ ) = h ∂ f∂x i i δ − ∂ ∂x i h f i δ , (A.17)can be obtained as in Sec. II [see (23)–(25)] C ′′ ii ( h f i δ ) = − X j =1 ∂ δ j ∂x i ∂∂δ j h f i δ − X j =1 ∂δ j ∂x i ( ∂ ∂δ j ∂x i h f i δ ) − X j,k =1 (cid:18) ∂δ j ∂x i ∂δ k ∂x i (cid:19) ∂ ∂δ j ∂δ k h f i δ . (A.18)Before deriving the approximated form for (A.18), whilewishing to mantain its fourth order of accuracy, it is use-ful to observe that the terms on the r.h.s. which contain ∂ ∂δ j ∂δ k h f i δ , k = j , by (A.8), are of the same order asthe remainder term. Therefore, they do not enter theapproximation, which is˜ C ′′ ii ( h f i δ ) = − X j =1 δ j ∂ i δ j + ( ∂ j δ j ) δ j h hh f i δ i δ j e j − h f i δ i − X j =1 ∂ i δ j δ j h h ∂ i h f i δ i δ j e j − ∂ i h f i δ i (A.19)Again, in analogy with what has been done for the firstderivatives, Eq. (A.16), the following is obtained˜ C ′′ yy ( h f i δ ) = − ϕ∂ y ϕ + ( ∂ y ϕ ) ϕ ( y ) h hh f i δ i ϕ ( y )∆ y − h f i δ i − ∂ y ϕϕ ( y ) h h ∂ y h f i δ i ϕ ( y )∆ y − ∂ y h f i δ i (A.20)2 [1] Ghosal S., Moin P. The basic equations for large eddysimulation of turbulent flows in complex geometries. J.Comp.Phys. , 24–37 (1995).[2] Van der Ven H. A family of large eddy simulation (LES)filters with non uniform filter widths.
Phys. Fluids , (5),1171–1172 (1995).[3] Fureby C., Tabor G. Mathematical and physical con-straints on large-eddy simulations Theoretical and Com-putational Fluid Dynamics (2), 85–102 (1997).[4] Vasilyev O.V., Lund T.S. and Moin P. A general classof commutative filters for LES in Complex Geometries. J.Comp.Phys. , 82–104 (1998).[5] Ghosal S. Mathematical and physical constraints onlarge-eddy simulation of turbulence.
AIAA Journal , (4), 425–433 (1999).[6] Sagaut P. Large eddy simulation for incompressible flows ,Springer Verlag, Berlin-Heidelberg (2001).[7] Germano M.,Piomelli U.,Moin P. and Cabot W.H. A dy-namic subdrid-scale eddy viscosity model.
Phys. FluidsA (7), 1760–1765 (1991).[8] Germano M. Turbulence, the filtering approach. J. FluidMech. , 325–336 (1992).[9] Bardina J., Ferziger J.H., Reynolds W.C. Improved sub-grid models for large eddy simulation.
AIAA Paper –1357 (1980).[10] Smagorinsky J. General circulation experiments with theprimitive equations. Monthly Weather Review , , 99–164 (1963).[11] G.Alfonsi, G.Passoni, L.Pancaldo and D. ZampaglioneA spectral-finite difference solution of the Navier-Stokesequations in three dimensions, Int. J. Numer. Meth. Flu-ids , , 129–142 (1998).[12] G.Passoni, G. Alfonsi, G.Tula and U.Cardu A wavenum- ber parallel computational code for the numerical inte-gration of the Navier-Stokes equations, Parallel Comput-ing , , 593–661 (1999).[13] Moser R.D., Kim J., Mansour N.N., Direct numericalsimulation of turbulent channel flow up to Re τ = 590, Phys. Fluids , (4), 943–945 (1999).[14] Moin P. and Kim J. Numerical investigation of turbulentchannel flow, J. Fluid Mech. , , 341–377 (1982).[15] Deardorff J.W. A numerical study of three-dimensionalturbulent channel flow at large Reynolds numbers, J.Fluid Mech. , , 435–452 (1970).[16] Schumann U. Subgrid scale model for finite differencesimulations of turbulent flows in plane channels and an-uli. J.Comput.Phys. , , 376–404 (1975).[17] Gr¨otzbach G., in Encyclopedia of Fluid Mechanics , vol.6,ed. N.P. Cheresimov, Gulf West Orange, New York(1987).[18] Piomelli U., Ferziger J. & Moin P. New approximateboundary conditions for large eddy simulations of wall-bounded flows.
The Physics of Fluids A (6), 1061–1068(1989).[19] Horiuti, K. Higher order terms in the anisotropic rep-resentation of Reynolds stresses, Phys. Fluids A , (19),1708–1710 (1990).[20] Carati D. and Cabot W. Anisotropic eddy viscosity mod-els. Proceeding of the summer program - Center for Tur-bulence Research, Stanford, 249–259 (1996).[21] Iovieno M., Tordella D. The angular momentum equationfor a finite element of fluid: a new representation andapplication to turbulent flows, Phys. Fluids , (8), 2673–2682 (2002).[22] Ghosal S. An analysis of numerical errors in Large-EddySimulations of Turbulence J. Comput. Phys. ,125