A modiﬁed selective structure function subgrid-scale model

. We propose a modiﬁed version of the selective structure function (SSF) model originally proposed by David (1993 PhD Thesis National Polytechnic Institute, Grenoble). As compared with the SSF model, the modiﬁed selective structure function (MSSF) model respects in a better way the energetic exchanges between the supergrid and the subgrid scales and automatically adjusts itself to the discretization thinness of the most energetic scales. Using the PRICELES industrial code, we ﬁrst simulate decaying homogeneous isotropic turbulence at high Reynolds number. The MSSF model is shown not to dissipate the supergrid scale energy when the kinetic energy is present at large scale only. It leads to a good prediction of the time decay law and to kinetic energy spectra with a well deﬁned k − 5 / 3 spectral range extending up to the cut-oﬀ wavenumber. The results of the MSSF model are found to be very close to those obtained with Smagorinsky’s dynamic model of Germano et al (1991 Phys. Fluids A 3 1760–5). For low Reynolds number isotropic turbulence, a very satisfactory comparison with the experimental data of Comte-Bellot and Corrsin (1971 J. Fluid Mech. 48 273–337) is obtained. In the turbulent channel ﬂow, the MSSF model was shown to allow for a much better prediction of the friction velocity owing to its less dissipative character in the near-wall region as compared to the SSF model.


Introduction
A subgrid-scale model which parametrizes the small scales in a Large-Eddy Simulation (LES) approach is properly designed when its action on the supergrid scales is limited to the most active turbulent flow regions.Smagorinsky's model in his original formulation (Smagorinsky [4]; see also Lesieur and Métais [5]) is known to be too dissipative when one uses a Smagorinsky's constant (called C S ) analytically determined by assuming that the cut-off scale between the subgrid and the supergrid scales is located within a Kolmogorov cascade (see Lilly [6]).For a Kolmogorov constant of 1.4, this analytical determination yields C S ≈ 0.18.Deardorff [7] in his pioneering LES of a turbulent channel flow took C S = 0.1, which represents a reduction by nearly a factor of four of the eddy viscosity.For this value, Smagorinsky's model behaves reasonably well for the turbulent channel flow with wall laws (Deardorff [7], Moin and Kim [8]) and for free-shear flows.Despite the reduction of C S , it can be shown (see, e.g., Meneveau and Katz [9]) that Smagorinsky's model does not allow for a proper decay of the turbulent viscosity near a wall.It then yields too large a dissipation in this region unless one uses an ad hoc wall function (see, e.g., Deardorff [7], Moin and Kim [8]).For this reason, it does not work for transition in a boundary layer developing upon a flat plate: it artificially relaminarizes the flow if the upstream perturbation is not high enough.Inspired by a spectral eddy-viscosity formulation, Métais and Lesieur [10] derived the so-called structure function (SF) subgrid-scale model based upon the second-order structure function of the velocity.Although this model yields a better behaviour than Smagorinsky's model for decaying isotropic turbulence (see Métais and Lesieur [10]), it also turns out to be too dissipative for transition in a boundary layer again yielding relaminarization.
These defects of the original models led to the development of new models which automatically adjust to the local flow conditions.In Smagorinsky's dynamic approach (Germano [11], Germano et al [2]), Smagorinsky's constant (which is no longer a constant) is determined in each point and at each instant through a double-filtering approach.Note that this dynamic procedure can also be applied to any eddy-viscosity model such as the SF model.The reader can, for instance, refer to El-Hady and Zang [12] for an application of the dynamic SF model to a compressible boundary layer above a cylinder.As far as spectral models are concerned, Lamballais et al [13] proposed the spectral-dynamic model consisting of a refinement of spectral eddy-viscosity models taking into account non-developed turbulence in the subgrid scales.This model is based on an instantaneous adjustment of the turbulent eddy-viscosity coefficient to the deviation, at small scales, of the spectral slope with respect to the standard Kolmogorov law.The selective (David [1]) and the filtered (Ducros et al [14]) SF models also constitute self-adjusting subgrid-scale models.In the filtered version, the SF is computed not from the raw supergrid scale velocity field but from the supergrid scale velocity field previously submitted to a high-pass filter: this pre-filtering operator removes the contribution to the eddy viscosity of the non-turbulent low-frequency oscillations.This model has allowed the successful study of the full transition to turbulence in a weakly-compressible spatially developing boundary layer over a flat plate (Ducros et al [14], Briand [15]).The present paper is focused on the selective structure function model originally proposed by David [1] for which we propose an improved version here.
The paper is organized as follows.Section 2 is devoted to the original selective structure function (SSF) model proposed by David [1].We first recall the SSF model formalism.The numerical methods used for the various tests are described and we present LES of decaying isotropic turbulence at high Reynolds number and compare the predictions of the SSF model with the classical SF model.An improved version of the SSF model is proposed in section 3: the modified selective structure function (MSSF) model.It is then applied to decaying isotropic turbulence at high Reynolds number and compared to various other models: the SF model, the SSF model and Smagorinsky's dynamic model (section 4.1).The case of low Reynolds number turbulence is next considered in section 4.2 by comparing the LES results with the experimental data of Comte-Bellot and Corrsin [3].Finally, section 5 is devoted to LES of the turbulent channel flow.

Methodology
We consider here a flow of constant density ρ 0 .Let Δx be a scale characteristic of the grid mesh.To eliminate the subgrid scales, we introduce a filter of width Δx such that the convolution of any quantity f ( x, t) by the filter function corresponds to the supergrid-scale field.The subgrid-scale field is the departure of the actual flow with respect to the filtered field: The application of the filter to the momentum equations yields the classical subgrid-stress tensor, which has to be modelled.u i stands for the ith component of the velocity field expressed in a Cartesian frame of reference.The SF model is the extension to physical space of the spectral subgrid-scale models based upon a spectral eddy viscosity scaled on E(k C )/k C where E(k C ) is the kinetic energy spectrum at the cut-off wavenumber k C (see Métais and Lesieur [10], Lesieur and Métais [5]).Since it is an eddy-viscosity model, one can then write: where Sij is the supergrid scale deformation tensor: To determine the turbulent eddy-viscosity coefficient ν SF t , Métais and Lesieur [10] use the local structure function constructed with the filtered velocity field u( x, t), where F 2 is calculated with a local statistical average of square (filtered) velocity differences between x and the six closest points surrounding x on the computational grid.Assuming a k −5/3 spectrum extending from zero to the cut-off wavenumber one obtains, through energy conservation arguments (see Métais and Lesieur [10] for details), where C K is the Kolmogorov constant.
The idea of the SSF model proposed by David [1] is to switch off the eddy viscosity in the regions where the flow is not three dimensional enough.The criterion of three-dimensionalization is defined as follows.First, a mean vorticity vector ω m is computed at each point M and at each instant.On a structured three-dimensional mesh, it is calculated by considering the arithmetic mean weighted by the inverse of the distance d(M, M i ) between the point M and the six closest points M i located around M : where ω(M i , t) is the instantaneous vorticity vector at point M i .Second, the angle α between the instantaneous vorticity vector at point M , ω(M, t), and ω m (M, t) is determined: David [1] carried out a LES of decaying isotropic turbulence at a resolution of 32 3 collocation points with a Kolmogorov k −5/3 spectrum extending from k = 0 to k = k C .It was found that the probability distribution function (PDF) of α peaks at a value of 20 • , which is thus the most probable value.Based on this observation, David [1] then proposed to cancel the eddy viscosity at points where this angle is smaller than 20 • , that is to say in regions where the flow is locally close to a two-dimensional state.As compared to the original SF model, this subgridscale model dissipates the supergrid scale energy at fewer points of the computational domain and the model constant of 0.105 (see equation (8)) has then to be increased to satisfy energy conservation.David [1] chose to impose to have the same spatially-averaged eddy coefficients for both the SF model and the SSF model, so that where the brackets stand for a spatial average over the whole computational domain.In the framework of his decaying isotropic turbulence LES, David [1] finally obtained the following eddy-viscosity formulation: where Φ 20 • ( x, t) is the indicating function based on the value of α: Although this model allows one to obtain good results for various incompressible and compressible turbulent flows (see, e.g., Lesieur and Métais [5], for a review), it presents some weaknesses which may render difficult its adaptation to very irregular meshes or to unstructured meshes.First, it seems obvious that the critical angle α c has to depend on the local numerical resolution.Indeed, for an infinitely refined resolution with a local grid size tending to zero, the vectors ω m ( x, t) and ω( x, t) are identical and the angle α c tends to zero.Second, it is well known that the global kinetic energy exchange between the supergrid and the subgrid scales for any eddy-viscosity model is given by If k C is assumed to be located in a Kolmogorov cascade, the loss of energy by the supergrid scales closely approximates the energy which cascades through the inertial subrange and is subsequently dissipated by viscosity.From the energetic point of view, a more strongly justified SSF model therefore requires that the global energy lost by the supergrid scales has to be identical when using the SF model or the SSF model.This leads to the following relation, which has no reason to be fulfilled by the SSF model satisfying equation (11).The modified SSF model that we propose in the next sections is aimed at correcting these deficiencies of the original SSF model.

Numerical algorithm
One of the motivations of the present work consists in testing the performances of our subgridscale models with an industrial code which could also be used for further applications in complex geometries of industrial interest.The simulations are then performed with the PRICELES code (Bieder et al [16]).The discretization can be either structured with the finite-difference volume version of the code or unstructured with the finite-element volume version.We only consider the structured version here.Pressure and velocity components are defined on different nodes of a staggered grid.The discrete form of the momentum and continuity equations for the supergrid scale field leads to a linear algebraic system.To solve this system, we use here a matrix projection method inspired by the SOLA algorithm originally proposed by Hirt et al [17].The reader is referred to Ackermann [18] for further details.The Poisson equation is solved by a conjugate gradient method with a SSOR preconditioning.It is now well recognized that LES cannot be properly performed with dissipative discretization schemes: a centred scheme is therefore used.Since the PRICELES code is designed to deal with very complex geometries at a computer cost as reduced as possible, the scheme is of second order only.It is well known that the differencing errors of such a scheme are large in the smallest resolved scales as compared with compact finite-difference schemes for instance and that its resolving efficiency is not very good (see, e.g., Lele [19]).Despite this deficiency, it will be shown in the next sections that quite satisfactory results are obtained concerning the decay of homogeneous isotropic turbulence and the turbulent channel flow.Note that the Laplacian operator within the Poisson equation for pressure is computed as the product of the divergence and gradient approximations used in the basic equations: this is important to ensure the numerical consistency of the operators (see Ferziger and Perić [20]).To check the conservative properties of the PRICELES code with such a numerical implementation, we have simulated the time evolution of an isotropic turbulent field in the case of zero molecular viscosity and zero turbulent viscosity, that is to say without any diffusion operator.The code was shown to conserve the total kinetic energy and to eventually yield a k 2 equipartition spectrum for the kinetic energy over almost all the simulated wavenumbers (see Ackermann [18] for further details).Finally, the temporal discretization scheme is a Runge-Kutta third-order scheme.

LES of freely-decaying three-dimensional isotropic turbulence
To get rid of the specific difficulties linked with solid boundaries, LES of decaying threedimensional isotropic turbulence constitute a classical test.Also geometrically simple, this flow configuration represents a quite severe validation for the numerical schemes and the subgrid-scale models.In the course of this paper, we will perform this classical test at moderate Reynolds number by comparing the kinetic energy decay law and the time evolution of the kinetic energy spectra with the laboratory experiment by Comte-Bellot and Corrsin [3] (see section 4.2).We also consider the large Reynolds number case by performing LES with near-zero molecular viscosity in which the supergrid scale energy dissipation is mainly attributable to the eddy viscosity.
The flow domain is a periodic cubic box of width L box = 2πL unit in each direction.Here, the molecular viscosity ν verifies ν ≈ 10 −21 L unit u unit where The computation is initially started with a three-dimensional random isotropic velocity field whose kinetic energy is concentrated around a given wavenumber located at large scale.The random field distribution corresponds to a Gaussian distribution.The initial kinetic energy spectrum is equal to and peaks at k = k i (0).Here, A is chosen such that the initial mean kinetic energy verifies The Eddy-Damped Quasi-Normal Markovian theory (EDQNM, see e.g.Lesieur [21] for details) predicts that an initial spectrum varying as k s at low wavenumber will immediately pick up a k 4 infrared behaviour when s > 4 (Lesieur and Schertzer [22]): for this reason, we chose s = 8 here.The characteristic time scale corresponds to the large-eddy turnover time defined as For an initial kinetic energy spectrum E(k, 0) ∝ k s when k → 0 with s ≥ 4, EDQNM predicts a t −1.38 decay law for the kinetic energy decay (Lesieur and Schertzer [22]).This decay law is only valid once a well defined Kolmogorov cascade has established.EDQNM-based calculations by André and Lesieur [23] have also shown that the necessary time t * for the cascade establishment is of the order of 5t ref .
For t < t * , the mean kinetic energy is conserved at high-enough Reynolds number, while the mean enstrophy D(t) grows significantly for t < t * before decreasing at a later time.Furthermore, the enstrophy maximum, which is reached at time t slightly larger than t * , increases with increasing Reynolds number (see also Lesieur and Ossia [24]).We now compare LES performed with the standard SF model given by equation ( 8) and the SSF model given by equation (12).The Kolmogorov constant C K is taken equal to 1.4.We take k i (0) = 10 and a spatial resolution of 64 3 cells with a regular mesh.Figure 1 shows the time evolution of the kinetic energy in the resolved scales defined by on a log-log plot.Since no energy is initially present at small scale, the kinetic energy is expected to remain constant until the energy spectrum starts to build up at small scale where it is dissipated.The SSF model correctly reproduces this trend with a constant level of supergrid scale kinetic energy until t * ≈ 2t ref .
During the same time, the supergrid scale enstrophy is seen to increase with a maximum which is reached at an instant slightly later than to 2t ref (figure 2).
Conversely, the SF model induces dissipation as soon as a velocity gradient is present at any scale: both the kinetic energy and the enstrophy start then to decrease from the beginning of the computation.As previously pointed out, the ability for the SSF model to dissipate only when small-scale velocity gradients are present is important in transitional flows to allow for a proper evolution from a laminar to a turbulent state.Note that the computed t * is much shorter than the value predicted by the EDQNM theory at very large Reynolds number.The work by Ossia [25] using spectral numerical methods of very high precision clearly demonstrates that t * varies with the numerical resolution and with the ratio k i (0)/k C .Ossia [25] obtained a value of t * very close to the present value in his spectral computations at similar resolution (64 3 ) and with k i (0) = 10.This constitutes a validation of our numerical code.The enstrophy maximum obtained by Ossia [25] with the spectral dynamic model of Lamballais et al [13] was, however, significantly larger, showing that the small-scale resolving efficiency of the second-order finite-difference scheme is much lower than for spectral methods.The long-time decay law of the kinetic energy is very close to the EDQNM prediction of t −1.38 for both the SF and the SSF models with a slightly faster decay for the latter.We next consider the supergrid scale enstrophy given by The smallest resolved scales are found to be located within a Komogorov cascade (see below) and, for such a spectrum, the main contribution to D comes from the smallest resolved scales.D(t) can therefore be roughly estimated by assuming a Kolmogorov spectrum extending from k = 0 to k = k C .This yields In the asymptotic decay law regime, D(t) is then expected to decay like 2/3 where is the dissipation rate of supergrid scale kinetic energy.If the supergrid scale kinetic energy E c (t) is varying like t −β , D(t) should therefore decay like t −(2/3)(β+1) (Bertoglio 2000 Private communication): for β = 1.38, one then gets ≈t −1.6 for the supergrid scale enstrophy decay law. Figure 3, which shows the time evolution of the supergrid scale enstrophy in a log-log plot, confirms this estimation.
We next consider, in figure 4, the kinetic energy spectra obtained with both the SF and the SSF model in the asymptotic decay regime: here t = 34t ref .As previously observed with spectral methods (see Métais and Lesieur [10]), the spectrum given by the SF model exhibits a well defined spectral slope ≈k −5/3 .Conversely, the SSF model yields a spectral slope which is not steep  enough indicating an insufficient dissipation at small scales.In the most energetic scales around the spectrum peak, a higher energetic level is obtained with the SF model than with the SSF model.Comparisons with the laboratory experiments by Comte-Bellot and Corrsin [3] presented in section 4.2 will indeed confirm that the SSF model tends to lead to an underestimation of the large-scale energy.

An improved version of the selective structure function model
As previously pointed out, the SSF model in its original formulation presents two main deficiencies: the model constant in formula ( 12) should be determined through energy conservation arguments given by the relation (15) and the critical angle α c should be resolution dependent.We use here the results of several LES of decaying isotropic turbulence performed with the standard SF model to determine the parameters of the MSSF model.

Critical angle α c
For any LES, one may intuitively apprehend that the critical angle α c will strongly depend upon the thinness of the numerical discretization of the large turbulent energetic scales.This resolution thinness is efficiently measured by the ratio between the cut-off wavenumber k C (characteristic of the filter width) and the wavenumber k i at which the spectrum peaks (characteristic of the large energetic scales).To demonstrate the dependence of α c with respect to the ratio k C /k i , we have performed several computations in which the resolution and the value of the initial spectrum peak k i (0) are varied.Tests with  resolutions ranging from 16 3 to 64 3 resolution points are then performed.We chose values of k i (0) going from 1 to 10.The molecular viscosity is still taken equal to a very small value.Since the results are obtained from the turbulent field at a given instant t fin during the self-similar period of decay, the pertinent parameter is then the ratio Figure 5 shows the PDF of the angle α c for different grid-point resolutions and different values of the ratio k C /k i (t fin ).It is important to note that α c depends only on the value of k C /k i and not on the total number of discretization points.Indeed, for a fixed value of k C /k i but for varying number of points, the various PDF are extremely close.As anticipated, α c decreases with increasing k C /k i (see figure 5 and table 1).To construct the analytical dependency of α c with k C /k i , the following observations can be made.First, our computations show that α c seems to exhibit an asymptotic value of ≈9 To summarize, we then propose the following analytical variation of α c with k C /k i : Figure 6 shows the variation of α c with k C /k i and the proposed fit: one can clearly see the closeness of the numerical data and the analytical approximation.

Model constant
Knowing the angle α c we may then derive the MSSF model which writes where Φ α c ( x, t) is given by and α c varies with k C /k i following relation (27).Next, the constant C M SSF has to be determined.We then perform several LES using the standard SF model at different resolutions and with different values of k i (0) (see table 1).The computed velocity field is then used to determine C M SSF by verifying the energy conservation constraint: Note that ν M SSF t in formula ( 30) is determined with α c derived from k C /k i (t fin ) through relation (27): the values of α c for the different runs are given in table 1.We have checked (see Ackermann [18] for details) that, whatever the resolution and the value of the ratio k C /k i , the different runs yield a value of C M SSF very close to 0.142.The MSSF model finally writes: with α c given by (27).Note that if α c is taken equal to 20 • , we recover the classical SSF model but with a reduced constant (0.142 instead of 0.172).This is consistent with the results of figure 1 showing a slightly too fast decay of the global kinetic energy with the SSF model, suggesting that the constant of the initial SSF model was indeed too high.

Isotropic turbulence at high Reynolds number
To check the validity of our new model, several LES of decaying isotropic turbulence are performed (still with C K = 1.4).The parameters for the various runs are gathered in table 1.
Note that k i decreases with time with a rapid decay during the early stage and a much slower decay during the self-similar stage.α c is then computed from an estimated value of k C /k i in the late decay stage (and not from k C /k i (0)).Based on the results of section 2.3, here we chose   k C /k i = 2k /k i (0) as indicated in table 1.We will see how to estimate this ratio in a more practical application to the turbulent channel flow (see section 5).Figures 7-10 compare the kinetic energy spectra obtained during the self-similar decaying regime with the classical SF model and the new MSSF model.The ratio k C /k i ranges from 3.2 to 8 and two resolutions of 32 3 and 64 3 resolution points are considered.Whatever the value of k C /k i and the resolution, the two subgrid-scale models give spectra which are very close at all scales.As compared to figure 4, the modified model yields a much more accurate behaviour at small scales than the SSF model.The slope with the MSSF model is indeed very close to the slope obtained with the standard SF model.
We now present a more detailed comparison between the three models (SF, SSF and MSSF) by focusing on the LES of decaying turbulence with k i (0) = 10 and 64 3 resolution points already described in section 2.3.The initial conditions are absolutely identical for the three computations with the three different subgrid-scale models.To determine the angle α c of the MSSF model, we take k C /k i = 2k C /k i (0) = 6.4 which yields α c = 11 • .Note that this corresponds to an angle almost twice as small as the 20 • angle of the SSF model.The MSSF model therefore dissipates the supergrid scale energy at more points of the computational domain than the SSF model and is then less selective in that respect.Conversely, as pointed out above, the constant of the MSSF model is lower.Figure 11 shows on a log-log plot the time evolution of the kinetic energy given by the three LES with the three different subgrid-scale models.As expected, the MSSF model yields an early conservation of the kinetic energy.However, the latter starts decreasing sooner than with the SSF model, which is a sign of the less selective character of the MSSF model.The modified model is also associated with a lower value of the enstrophy peak as shown in figure 12.During the self-similar period, the slightly too fast decay rate of the SSF model is   corrected by the MSSF model with a decay law closer to the EDQNM prediction.The spectra given by the three models at t = 34 t ref are plotted in figure 13.The SF and MSSF models lead to a very similar spectral behaviour at all scales.The kinetic energy build-up at small scales observed with the SSF model is now absent in the MSSF LES with a much better concordance with a k −5/3 behaviour.Conversely, more energy is present in the most energetic scales with the modified version of the model.
We finish this section devoted to high Reynolds number homogeneous isotropic turbulence by comparing the predictions of the MSSF model with those of Smagorinsky's dynamic model (Germano [11], Germano et al [2]).As previously stated, the dynamic procedure consists in determining Smagorinsky's constant C S at each point and at each instant through a doublefiltering approach.Following Lund's [26] recommendation, here the width of the second filter is taken to be equal to √ 6Δx (see Ackermann [18] for details).The allowance of negative values of C S constitutes a conceptual advantage of the model because these values represent a sort of backscatter in physical space, that is to say the kinetic energy is able to be locally transferred from the subgrid towards the supergrid scales.However, very large negative values of the eddy viscosity is a destabilizing process in a numerical simulation and may lead to a non-physical growth of the supergrid-scale kinetic energy (Lund et al [27]).Furthermore, the very sharp fluctuations of the model constant make the model formalism inconsistent.Several more or less sophisticated techniques have been proposed to overcome these large variations (see, e.g., Sarghini et al [28]).For the sake of simplicity, here we simply performed at each point a local spatial averaging of C S on the six closest points of the computational mesh.We also used the clipping technique consisting of cancelling the values of C S smaller than zero and larger  than 0.5.Figure 14 is the equivalent of figure 11, but now compares the temporal evolution of the kinetic energy as given by the SF, MSSF and Smagorinsky's dynamic models.The MSSF model performs slightly better than the dynamic model in the early stage since the energy is conserved during a slightly longer time with the MSSF model.This too dissipative nature of the dynamic model is confirmed by figure 15: the enstrophy peak obtained with the dynamic model indeed has a lower amplitude.In the self-similar decay period, the three models yield very close decay laws.This is confirmed by figure 16 showing the kinetic energy spectra obtained with the three models at t = 34t ref .The three spectra are very close, with a well defined k −5/3 range.The energy in the most energetic scales is slightly larger in the case of the dynamic model.It is important to note that the MSSF model yields, at least for decaying isotropic turbulence, slightly better results than the dynamic model at early times and very close results at later times.Smagorinsky's dynamic model requires more operations and memory than the MSSF model and it is encouraging to see that a conceptually simpler model yields very close results.

Comparison with grid-turbulence experiments
We next consider low Reynolds number decaying isotropic turbulence by comparing our LES results with the classical experiment on decaying turbulence behind a grid by Comte-Bellot and Corrsin [3].It consists in measurements in the turbulent flow downstream of a grid of spacing M = 5.08 cm oriented normal to a uniform steady flow of velocity U 0 = 10 3 cm s −1 .In a reference frame moving with the average flow velocity the flow can be considered as decaying isotropic turbulence.In the experiment, statistical data were collected at three downstream   18)) that we have estimated to be t ref = 10.4M/U0 .The origin of x and of the corresponding time t is taken at the first experimental station (location x 1 and time t 1 ), so that the correspondence between (x−x 1 )/M and (t−t 1 )/t ref is given by (x−x 1 )/M = 10.4(t−t 1 )/t ref .
The initial field is obtained by running the computation long enough to allow for the skewness factor S to build up to values typical of laboratory experiments S ≈ −0.5 and for the kinetic energy decay law to be close to the experimental law t −1. 25 .The velocity field so obtained is then normalized in such a way that its kinetic energy spectrum matches the experimental spectrum at x/M = 42.Note that a more sophisticated initialization technique has been proposed by de Bruyn Kops and Riley [29] that we have not deemed necessary to test owing to the low order of our numerical techniques.We compare here the standard SSF model and the MSSF model.Since the peak wavenumber k i (0) of the experimental spectrum is ≈2, the critical angle α c of the MSSF model is taken equal to 9 • .Here we use 64 3 discretization points.Figure 17 shows the decay of the energy in the resolved scales as a function of time.The two models give consistent results which agree well with the experiment.As previously noticed in the high Reynolds number case, the SSF model induces a too rapid decay between x/M = 98 and x/M = 171 as compared with the experimental points.This trend is corrected by the MSSF model.We next compare the experimental spectra and the LES spectra at the three experimental locations corresponding respectively to t−t 1 = 0, t−t 1 = 5.38t ref and t−t 1 = 12.4t ref .As pointed out by de Bruyn Kops and Riley [29], a quite stringent large-scale resolution requirement has to be met to avoid the energy to be removed too rapidly from the large scales.This constraint is not fulfilled by the present LES and one observes that the peak of the spectrum moves to the left more slowly than that of the laboratory flow (see figure 18).This is particularly noticeable at the third experimental station.However, at the second station, the MSSF model agrees quite satisfactorily with the experimental data and yields better results than the SSF model.Similarly to the high Reynolds number case, the SSF model induces too much energy removal from the most energetic scales.

The turbulent channel flow
Since engineering flows often occur in the presence of solid boundaries, we now test the behaviour of the MSSF model in a turbulent boundary layer.We simulate the classical turbulent channel flow, that is to say the turbulent flow between two infinite parallel plates.x, y and z are respectively the streamwise direction, the direction normal to the channel walls and the spanwise direction.u( x, t) = (u, v, w) designates the velocity vector and its components along x, y and z.In the present LES, we impose the mass flux to be constant so that the instantaneous bulk velocity U m is identical in the various runs.We then define the macroscopic Reynolds number Re = U m 2h/ν based upon the bulk velocity, the width of the channel 2h and the kinematic viscosity ν.As far as statistics are concerned, the brackets • indicate an average over planes parallel to the walls (homogeneous directions) and time.The rms values are determined from the fluctuating quantities about this average.In the following, most quantities are expressed in wall units and the superscript + indicates the conventional boundary-layer scaling by the friction velocity u τ and the viscous thickness δ v : with h + = u τ h/ν.

Determination of the critical angle α c of the MSSF model
As previously shown for the isotropic turbulence case, the value of the critical angle α c strongly depends on the resolution thinness of the most energetic scales.In spectral space and for homogeneous turbulence, it can be characterized through the ratio k C /k i .For computations carried out in physical space, k C /k i has to be related to quantities which can be easily derived in this space.The integral scale L of the turbulence can be approximated as Since For the channel flow, the mixing length l m provides a good estimation of the large energetic scales.It is well known that, close to the wall, the large turbulent eddies are limited by the distance y to the wall.Following Prandtl's linear-variation assumption, we then have l m ≈ κy where κ = 0.41 is the Von Karman constant.The channel flow laboratory experiments and numerical simulations (see Kravchenko et al [30]) show that l m stops being limited by the wall for y 0.2h.We then chose the following final expression for l m (y): It is important to note that the critical angle α c is no longer identical for all the points of the computational domain, as in the case of isotropic homogeneous turbulence, but it is now a function of the distance y from the wall.

Numerical parameters
We present LES results of the turbulent channel flow at Re = 6666.The size of the computational domain, in the streamwise, wall-normal and spanwise directions, is 2πh×2h×πh and a moderate resolution is used with 64 × 65 × 32 grid points.Two distinct LES are performed: the first with the original SSF model and the second with the MSSF model.The two simulations respectively correspond to h + = 204 and to h + = 209.We impose the physical no-slip conditions at the wall and the flow in the near-wall region is explicitly computed.To reach a better description of this region, the physical space (x, y, z) is transformed into the computational space (x, ξ, z) by a hyperbolic tangent transformation: In the transformed coordinate system, the grid points are uniformly spaced: N y being the number of discretization points in the direction normal to the wall.a is an adjustable parameter (0 < a < 1) which determines the points concentration: here a is taken to equal 0.9018.The first computational point away from the wall is located at y + ∼ 1.
For Re = 6666 (based on the bulk velocity), the regime is subcritical from the point of view of linear-stability theory.In this case, it would require a very long time integration to reach a realistic turbulent state if one starts from an initial random noise.The critical Reynolds number for a rotating channel flow with the rotation vector orientated along the spanwise direction is much lower with Re ≈ 100.To economically obtain fully developed turbulent initial conditions, we then initially impose a moderate Coriolis force.The turbulent state so-obtained is not fully symmetric with respect to the channel middle plane (see Lamballais et al [13]): the Coriolis force is then removed and a totally symmetric state is reached after some more integration in time.

Results
In figure 19, the mean longitudinal velocity profiles for the two cases (SSF and MSSF models) are compared with the LES of Piomelli [31] using Smagorinsky's dynamic model (Re = 6500, h + = 203).Note that, although a rather coarse resolution was used by Piomelli (48×65×64 grid points), precise Fourier-Chebychev pseudo-spectral collocation methods were employed and a very good agreement with experimental data was obtained.Our results are also compared with the standard theoretical wall laws: the linear law U + = y + of the viscous sublayer and the logarithmic law where κ = 0.41 is the Von Karman constant.By comparison with the LES carried out with the SSF model, the LES with the MSSF model compares better with Piomelli's data at least in the logarithmic layer.This indicates that the modified model yields a better prediction of the friction velocity at the wall.A possible explanation for this better behaviour of the MSSF model as compared to the SSF model may be the less dissipative nature of the former in the   near-wall region.This is illustrated by the ratio r tv (y) of the two averaged turbulent eddy viscosity coefficients obtained in the two LES, shown in figure 20.Close to the wall the SSF model is much more dissipative with a value of r tv close to 4. In the fully developed turbulent core of the channel, the ratio is slightly larger than unity, indicating a comparable dissipation induced by the two models.It is a well known that the discretization errors associated with low-resolution computations based upon second-order finite-difference schemes can be large.In channel flow LES, a symptomatic effect of these errors is an overprediction of the amplitude of the u rms peak and an underprediction of both the v rms and w rms peaks.Figure 21 shows the profiles of u rms , v rms and w rms compared with Piomelli's data.Considering the low resolution used, the overprediction of the u rms peak is not large with u max rms ≈ 2.99 (in wall units) for the SSF model and u max rms ≈ 2.9 for the MSSF model.As compared with the classical value u max rms ≈ 2.7, the MSSF model yields a 7.4% overprediction against 10.8% for the SSF model which indicates a slight improvement in the fluctuations statistics by the former model.Note that the pseudo-spectral LES computations carried out by Lamballais et al [13] on the basis of the spectral-dynamic model led to a value of the u rms peak very close to the current value with u max rms ≈ 2.8.Better results would be obtained by increasing the resolution of the present LES.However, second-order methods are obviously associated with a slow decrease of the discretization error with decreasing Δx.A more rapid decrease could be reached with the use of a fourth-order method as confirmed by the recent LES by Méri and Wengle [32].It is well recognized that the position of the u rms peak is very sensitive to the subgrid-scale model.Both the SSF and MSSF models predict a value very close to the expected value y + max ≈ 12: here y + max ≈ 13.5 corresponding to an overprediction of only 12.5%.

Conclusion
In recent years, the effort in LES has been largely devoted to the development of subgrid-scale models clever enough to automatically adapt themselves to the local flow characteristics and to remove (or add) energy from (to) the supergrid scales at the proper locations.This is of particular importance for transitional flows for which the transition process may be inhibited by the action of the subgrid-scale model.Starting from the SF model originally proposed by Métais and Lesieur [10], David [1] proposed the SSF model which acts selectively in the regions of highenough three-dimensional turbulent activity.The three-dimensionalization criterion consists in measuring the angle of the vorticity at a given grid point and the average vorticity on the closest neighbouring points.If this angle exceeds a critical value α c , the eddy viscosity is locally turned on.We have constructed an improved version of the SSF model (the MSSF model) which better respects the energetic exchanges between the supergrid and the subgrid scales.Furthermore, the angle α c was shown to be a function of the discretization thinness of the most energetic scales measured by the ratio k C /k i .We have then proposed an analytical variation allowing this dependency to be taken into account.Although the use of second-order schemes for LES is still controversial, successful second-order accurate LES have been performed for a wide variety of flow configurations (see, e.g., Akselvoll and Moin [33]).Here, we have deliberately chosen to test this model in the PRICELES code which is an industrial software designed to deal with complex flow geometries: only second-order finite-difference schemes have then been considered.Decaying homogeneous isotropic turbulence at high Reynolds number with very small molecular viscosity was first simulated.When the kinetic energy is present at large scale only, the MSSF model was shown not to induce any energy dissipation until the energy cascades to scales of the order of the mesh size.This property may be important for transitional flows for which the energy must not be dissipated by the subgrid-scale model before some turbulent energy is present at small scale.In the self-similar decay regime, the time decay law was found to be in good agreement with the EDQNM predictions.In this regime, the kinetic energy spectrum given by the MSSF model exhibits a well defined k −5/3 behaviour extending until the cut-off wavenumber: this represents an improvement as compared with the SSF model which leads to an energy accumulation in the smallest scales.We have checked that the results of the MSSF model are very close to those obtained with Smagorinsky's dynamic model proposed by Germano et al [2] with a slightly better behaviour of the MSSF model in the early non-dissipative stage of the evolution.We have then tested the performance of the MSSF model by comparison with the low Reynolds number decaying turbulence experiments by Comte-Bellot and Corrsin [3].The simulated decay law was found to be in good agreement with the experimental data and the MSSF model proved to lead to a better spectral representation of the large energetic scales than the SSF model.The final test concerned the turbulent channel flow; the modification of the SSF model was shown to allow for a much better prediction of the friction velocity and, consequently, for a much better reproduction of the mean longitudinal velocity profile in the logarithmic region.This is attributable to the less dissipative character of the MSSF model in the near-wall region.Despite the low resolution of the channel flow computations and the low order of the numerical scheme, the prediction of both the amplitude and the location of the maximum of longitudinal velocity fluctuations was quite satisfactory.These results are very encouraging concerning the LES of flows in complex geometries using industrial codes.Only simple geometries have been considered here.However, the implementation of the MSSF model to more complex geometries only requires a rough estimate of the variation in size of the integral scale with the distance from the obstacle to improve the results with respect to the SSF model.This will be the topic of further studies.

Figure 1 .
Figure 1.LES of decaying isotropic turbulence at high Reynolds number with k i (0) = 10 and 64 3 grid points.Time evolution of the supergrid kinetic energy with the SF model and with the SSF model on a log-log plot.

Figure 2 .
Figure 2. LES of decaying isotropic turbulence at high Reynolds number with k i (0) = 10 and 64 3 grid points.Time evolution of the supergrid enstrophy with the SF model and with the SSF model on a linear-linear plot.

6 Figure 3 .
Figure 3. LES of decaying isotropic turbulence at high Reynolds number with k i (0) = 10 and 64 3 grid points.Time evolution of the supergrid enstrophy with the SF model and with the SSF model on a log-log plot.

3 Figure 4 .
Figure 4. LES of decaying isotropic turbulence at high Reynolds number with k i (0) = 10 and 64 3 grid points.Kinetic energy spectra at t = 34t ref with the SF model and with the SSF model on a log-log plot.

Figure 6 .
Figure 6.Variation of the critical angle α c as a function of k C /k i .The points are the values obtained through the various LES and the lines represent the analytical fit given by equation (27).

Figure 7 .
Figure 7. LES of decaying isotropic turbulence at high Reynolds number with k C /k i (t fin ) = 3.2 and 32 3 grid points.Kinetic energy spectra at t = 35t ref with the SF model and with the MSSF model on a log-log plot.

Figure 8 .
Figure 8. LES of decaying isotropic turbulence at high Reynolds number with k C /k i (t fin ) = 6.4 and 32 3 grid points.Kinetic energy spectra at t = 34t ref with the SF model and with the MSSF model on a log-log plot.

Figure 9 .
Figure 9. LES of decaying isotropic turbulence at high Reynolds number with k C /k i (t fin ) = 6.4 and 64 3 grid points.Kinetic energy spectra at t = 34t ref with the SF model and with the MSSF model on a log-log plot.

Figure 10 .
Figure 10.LES of decaying isotropic turbulence at high Reynolds number with k C /k i (t fin ) = 8 and 64 3 grid points.Kinetic energy spectra at t = 35t ref with the SF model and with the MSSF model on a log-log plot.

Figure 11 .Figure 12 .
Figure 11.LES of decaying isotropic turbulence at high Reynolds number with k i (0) = 10 and 64 3 grid points.Time evolution of the supergrid kinetic energy with the SF model, the SSF model and the MSSF model on a log-log plot.

Figure 13 .
Figure 13.LES of decaying isotropic turbulence at high Reynolds number with k i (0) = 10 and 64 3 grid points.Kinetic energy spectrum at t = 34t ref with the SF model, the SSF model and the MSSF model on a log-log plot.

Figure 14 .Figure 15 .
Figure 14.LES of decaying isotropic turbulence at high Reynolds number with k i (0) = 10 and 64 3 grid points.Time evolution of the supergrid kinetic energy with the SF model, the MSSF model and Smagorinsky's dynamic model on a log-log plot.

Figure 16 .
Figure 16.LES of decaying isotropic turbulence at high Reynolds number with k i (0) = 10 and 64 3 grid points.Kinetic energy spectra at t = 34t ref with the SF model, the SSF model and Smagorinsky's dynamic model on a log-log plot.

Figure 17 .
Figure 17.LES of decaying isotropic turbulence at low Reynolds number, comparison with the experiment by Comte-Bellot and Corrsin [3].Time evolution of the supergrid kinetic energy with the SSF model and the MSSF model compared with the experimental data on a linear-linear plot.

Figure 18 .
Figure 18.LES of decaying isotropic turbulence at low Reynolds number, comparison with the experiment by Comte-Bellot and Corrsin [3].Kinetic energy spectra at the three experimental locations: * , x/M = 42; •, x/M = 98; ×, x/M = 171.The results of the SSF model (dashed curve) and of the MSSF model (full curve) are compared with the experimental data on a linear-linear plot.

Figure 19 .
Figure 19.Turbulent channel flow.Mean velocity profile in wall units: comparison of the SSF and MSSF model predictions with Piomelli's [31] results.The theoretical linear and logarithmic laws are also plotted on a linear-log plot.

Figure 20 .
Figure 20.Turbulent channel flow.Profile of the ratio of the average eddy-viscosity coefficients given by the SSF model and the MSSF model (see equation (39)).

Figure 21 .
Figure 21.Turbulent channel flow.Rms of the three velocity components in wall units: comparison of the SSF and MSSF model predictions with Piomelli's [31] results on a linear-linear plot.

Table 1 .
Parameters of the various LES.Resolution k i (0) k C /k i (t fin ) Angle α c ( • )