Evaluation of sub grid scale and local wall models in Large-eddy simulations of separated flow

The performance of the Sub Grid Scale models is studied by simulating a separated flow over a wavy channel. The first and second order statistical moments of the resolved velocities obtained by using Large-Eddy simulations at different mesh resolutions are compared with Direct Numerical Simulations data. The effectiveness of modeling the wall stresses by using local log-law is then tested on a relatively coarse grid. The results exhibit a good agreement between highly-resolved Large Eddy Simulations and Direct Numerical Simulations data regardless the Sub Grid Scale models. However, the agreement is less satisfactory with relatively coarse grid without using any wall models and the differences between Sub Grid Scale models are distinguishable. Using local wall model retuned the basic flow topology and reduced significantly the differences between the coarse meshed Large-Eddy Simulations and Direct Numerical Simulations data. The results show that the ability of local wall model to predict the separation zone depends strongly on its implementation way.


Introduction
Large Eddy Simulations (LES) are widely used to simulate flows that are dominated by large-scale structures such as the flows that are often encountered in atmospheric boundary layer and wind energy applications.In order to use the LES models in these applications, they should be able to deal with the anisotropic nature of the flow over terrains and the near surface turbulence.Resolving the buffer and viscous sub-layers requires a mesh resolution that scales with Reynolds number.This is often not achievable in such applications where the Reynolds numbers are very high and therefore the traditional no-slip boundary conditions cannot be used.Instead, approximate boundary conditions are devised to represent the effects of the unresolved flow at the surface neighborhood on the outer flow.The need of approximate boundary conditions for LES of high Reynolds number wall bounded flows was recognized in the early stage of the LES development (see for example [1,2]).A variant of a wall model [3][4][5][6] is derived based on the asymptotic behavior of the attached near wall flow at high Reynolds number where the averaged velocity shows a logarithmic profile.However, for an inhomogeneous configuration there is no obvious way to evaluate the mean wall stress.In the absence of any known rigorous formulation which will hold in this case, the law of the wall is adapted by enforcing it locally and instantaneously.The use of the law of the wall locally can be motivated by the short lifetime and length scale of the near wall vortices.If the cell dimensions are chosen to be much larger than the length scale of near wall vortices in conjugate with a time step that is as large as many vortices lifetime, these local values can be seen as the statistical average of the effect of many vortices.Masson and Callen [7] reported that the validation of this approach depends strongly on the grid size and the statistical base of this approach can fail if too fine grid is used.This approach is frequently used in LES for atmospheric flow over complex terrain [see for example 8-11].However, it has also been shown that the use of a wall model in conjunction with coarse grid can lead to wrong evaluation of wall stress because of poorly predicted flow in the vicinity of the wall [12,13].The use of coarse grid become seriously questionable in case of a separated boundary layer, since the details about the growth of the boundary layer and its separation depend to a large extent on the nearwall momentum, which may be distinctly inaccurate using wall stress models.
The goal of this study is to investigate the performance of a number of the sub grid scale models implemented in OpenFOAM [14] and the legitimacy of using a local wall model for separated flow of a wavy channel.The question to be answered here is whether the use of local wall model can reproduce correctly the separation, reattachment points and the statistically averaged quantities or not.

Sub Grid Scale Models
In LES, the computational cost of solving all the turbulence scales is reduced by only solving the most energetic part of the turbulent energy spectrum.The LES are based on the assumption that for high Reynolds number flow the dependent variables can be decomposed into large-or Grid Scale (GS) components and small-or Sub Grid Scale (SGS) components, which represent the unresolved fraction of turbulence.Classically, the separation between scales is done by applying a convolution filter of low pass type to the unsteady Navier-Stokes Equations (NSE).In the absence of body forces, the incompressible LES are governed by: The convolution process generates an additional term, the SGS stress tensor: The SGS stress tensor includes unfiltered quantities and therefore it cannot be computed directly and it needs to be modeled.A widely used approach to model the effect of the SGS stress on the filtered field, which is mainly considered as energy drain, is based on the eddy viscosity hypothesis that relates the eddy stress to mean flow gradient by assuming that the deviatoric part of the SGS stress tensor is linearly proportional to the filtered strain rate tensor: where ߜ is the kronecker delta, ܵ ̅ is the filtered strain rate tensor and ‫ݒ‬ ் is SGS eddy viscosity, which is unknown quantity that must be specified for the model to be closed.
The best-known model of eddy viscosity type is the one proposed by Smagorinsky [15]: where ‫ܥ‬ ௦ is the Smagorinsky constant and ∆ is the filter width.The other widely used approach also based on the eddy viscosity hypothesis is the one-equation model (see for example [16]).In the oneequation model, a balance equation of SGS kinetic energy ݇ ௌீௌ is derived by constructing the exact balance equation for SGS stress tensor: Where D is the diffusive effect and assumed to be equal to ݇ ௌீௌ ൰ where ‫ݒ‬ is the effective viscosity ( ‫ݒ‬ = ‫ݒ‬ ௌீௌ + ‫.)ݒ‬The last term in the Eq 5 represents the dissipative effect and it is modeled by ߝ = ‫ܥ‬ ݇ ଷ ଶ ⁄ /∆.In the above formulation the SGS stress tensor is modeled according to the eddy viscosity Eq 3 and the SGS eddy viscosity is evaluated as ‫ݒ‬ ் = ‫ܥ‬ ∆√݇ .The ‫ܥ‬ and ‫ܥ‬ are model constants.
The determination of model constants in the above formulation of SGS stress tensor requires knowledge of turbulence nature.Since the only available theoretical analysis for the turbulence is for homogeneous isotropic turbulence, most of these constants are derived within the framework of this theory.To deal with a more generic form of turbulence, the SGS model constants are modified often by an ad hoc manner.The error rising from tuning the SGS model can be reduced by introducing dynamic procedures that calculate the constants from the flow features.Germano [17] suggested a dynamic procedure providing a systematic way to calculate the value of the Smagorinsky coefficient at every time and position in the flow based on the dynamics of the smallest resolved scale.This procedure is based on the Germano identity, which links the SGS stress tensor to the equivalent tensor obtained at larger filtering width.The Germano dynamic procedure requires smoothing to guarantee numerical stability and to avoid excessive fluctuation in the model coefficient that might be result from the dynamic procedures.Typical averaging can be over direction of homogeneity or in case of complex terrain over flow pathlines using Lagrangian averaging [18,19].A similar dynamic idea was applied to evaluate the ‫ܥ‬ and ‫ܥ‬ in the one-equation model.Other modifications are also used to deal with the increase in anisotropy in both the resolved and SGS velocities due to strong mean shear near the surface.The most common one is the damping functions of van Driest type which ensure that the SGS viscosity vanishes as the wall is approached as the one used by Moin and Kim [20] and the wall damping correction proposed by Mason and Thomson [21].
An alternative approach to model the SGS stress is to construct their components from the filtered quantities by deriving their balance equation from the filtered NSE: where P is the production, M is the generalized triple correlation, Π is the pressure velocity gradient tensor, and E is the dissipation tensor.Advantage of such approach is more likely to be able to deal with the flow or grid anisotropity.Deardorff [22] proposed a first successful model of that type.He modeled the terms in Eq 6 by where Ρ = ൫‫߬ܮ‬ ் + ‫ܮ‬ ் ߬ ൯, L here is the Leonardo stress, and Much effort has been devoted to SGS modeling results in variant of models.We limited the above description to only the models have been used in this study.For further information about other SGS models, the reader is referred to Sagaut [23].

Wall models:
The law of the wall is derived by neglecting the acceleration, pressure gradient and the viscous effects in the streamwise momentum equation at the first grid point above the wall.Further, by assuming constant shear stress between the wall and first grid point, the integration of the momentum equation leads to either the form: Where ‫ݑ‬ * is the friction velocity ‫ݑ(‬ * = (߬ ௪ ߩ ⁄ ) ଵ ଶ ⁄ ), here ߬ ௪ , ߩ are the wall shear stress and the fluid density respectively, ݇ is the von Karman constant, ܼ ା is the normal distance from the wall in wall units (ܼ ା = ‫ݑܼ‬ * ‫,)ݒ/‬ B is a constant and ܼ is the aerodynamic roughness.This velocity profile is the same as the one that is widely used in atmospheric studies that resulting from the Monin-Obukhov [24] similarity theory in case of natural stratified flow.
In this study an approach similar to the one suggested by Schumann [4] and Grötzbach [5] but in local manner is used where the instantaneous velocity between the first off-wall grid point and the wall itself assumed to have a logarithmic profile Eq 8.The instantaneous wall shear stresses for a given velocity at first off-wall grid points, which then serve as a wall boundary condition for outer LES domain, are zero except for ߬ ௫௭ and ߬ ௬௭ components and are estimated as following: where x, y, and z are streamwise, spanwise and normal directions respectively, ‫ݖ‬ ଵ denotes the height of the center of the cell adjacent to the wall and the angle brackets denote a horizontal average.The friction velocity is calculated from Eq 8.The localization of the model is done by replacing the plane average velocity in equation 8 and 9 by local velocity magnitude.In complex terrain the local wall stresses are computed by combining Eq 8 and Eq 9 [8] : where ܷ ഥ is the magnitude of the tangential velocity calculated at the first off-wall grid points.ߠ ௫ , ߠ ௬ are the local angles of inclination of the topography in x and y directions.
The calculated wall shear stress can be added to the SGS stress tensor term in the momentum equation.Computing the eddy viscosity at the center of the first computational cells along the wall would be a problem in this case, because it requires the values of stress and strain tensors which are not totally correct at the wall since it is not no-slip condition and hence the horizontal velocities are not actually specified.This problem can be remedied by averaging the no-slip values of the one-sided velocity differences at the first off-wall point from interior [12].Another way to estimate the SGS eddy viscosity is to evaluate it at midway between the first off-wall cells and the cell directly above it.In this study we will use the later method (hereafter we call it LW1).
Another simple way to implement the calculated wall shear stress is to modify the surface SGS viscosity and the multiplication of it with the normal gradient of velocity at the wall equals the wall shear stress, we will call this model as LW2.The wave amplitude to wavelength ratio is 0.05.The NSE are solved numerically using the finite volume method with second-order central difference discretization scheme for space and second-order backward implicit scheme for time discretization.The pressure-momentum system is decoupled using the PISO method.Highly-resolved LES are performed first without wall model by using three different SGS models: Smagorinsky (Smag) [15], One Equation model (One) [16] and Deardorff differential model (Deard) [17], where the λ are discretized by using 64 grid points for both x-and y-directions and 96 grid points are used in z-direction.The mesh is stretched slightly toward the channel center with (∆z i+1 =1.022*∆z i ) ratio.To isolate the effects of SGS model from the wall model effects, the above mentioned SGS models are then used to simulate the flow on relatively coarse grid with no-slip boundary condition (no wall model is used).The mesh resolution is reduced by factor three in all directions where (20 20 32) equidistantly distributed grid points are used for each λ in respectively directions.To examine the effects of using dynamic models on the results, the dynamic Lagrangian Smagorinsky (dynLag) [18] and the dynamic one equation model (dynOne) [15] are also used.Two additional cases are simulated on the same mesh by using Smagorinsky with van Driest damping function (SmagDamp) [25] and Spalart Almaras (SA) Detached Eddy Simulation (DES) models [26].

Figure1. The flow domain
Finally the effectiveness of using local wall models (LW1 and LW2) with both Smag and dynLag are studied.SA model combined with Nut-wall function (SA1) [15] and Spalding wall function (SA2) [27] are also studied.In this study the SGS models are used in their original implementation form in OpenFOAM 2. 1.x [15] with their defaults constants.The flow is initialized from well resolved LES and the statistical quantities are taken under the last 20 seconds (40s < t < 60s).  ) , are less than 7% (Table 1).The three SGS models are nearly indiscernible where the modeled SGS energy in this case is a small fraction of the total turbulent energy and the modeling approach has no significant effect.The three SGS models predict the separation and reattachment points at almost the same (x/λ) positions, within the mesh streamwise resolution marginal which is 1/64 in this case, (Table 1).Table 1.Modeled parameters of various LES cases. in all cases the computational domain extends over (6λ 2λ 1λ) in x,y and z directions.N denotes the number of points for each λ, (x/λ)Sep and (x/λ)Reat stand for the stream wise position of separation and reattach points, respectively, normalized by wavelength.(u*/U) is the normalized friction velocity at the wavy wall.Ξ is the relative L2 norm of the differences between the LES and DNS data.

E3S Web of Conferences
The turbulent kinetic energy (k) and the vertical shear stress ‫‪ത‬ݓ́ݑ(‬ തതത ) are also in good agreement with DNS data.The phase averaged profiles of (k) and ‫‪ത‬ݓ́ݑ(‬ തതത ) at different (x/λ) positions are shown in figure [4].The LES results show the same trend as DNS data and have the profiles maxima at almost the same positions as DNS which indicates that the LES correctly determine the position of shear layer.However, the LES underestimate the (k) at all (x/λ) positions.The largest Ξ(k) is found in the separated zone and it reduces to its minimum value near the reattachment point and increases again within the developed boundary layer at the upslope of the wave crest.This special variations can be caused by the intensive mixing and the additional turbulence production across the separation zone due to the strain rate induced by stream lines curvature, while at the reattachment point the interaction between the shear layer and the wall causes a strong modification in pressure fluctuation which in turn causes intensive turbulent energy redistribution among the stress component and reduces the stress anisotropty.The ‫‪ത‬ݓ́ݑ(‬ തതത ) predicted by LES is smaller than the DNS near the wavy surface in the leeward side of the wave crest where the flow decelerate and the velocity gradient is small, while LES results overestimate it due to the high velocity gradient in the thin boundary layer downstream the reattachment points.The Deard model predicts superior results to those obtained with SGS models based on the eddy viscosity hypothesis (Table 1).

Low-resolution LES
The differences between SGS models are more pronounced in the coarse mesh where SGS models are supposed to model a considerable part of turbulent kinetic energy.Without using a wall model, all the tested SGS models failed to predict the separation zones.The results show that the Smag model overestimates the velocity near the wall in all (x/λ) positions (Fig 5 left).This can be due to the use of the default value of Smag constant (Cs=0.17The theoretical value of Smag constant was determined by Lilly [29] to be (Cs=0.23)for homogeneous isotropic turbulence with cutoff in the inertial subrange and ∆ equal to the grid size.The constant is often reduced to improve the model performance at high shear flow near the wall.
Deardorff [22] used the value of (Cs=0.1)with filter width equal to grid size, while Mason and Callen [7] used 0.2 for sufficiently fine mesh and found that this value should be reduced if the mesh is coarse.The same behavior is noticed by using the default values of the One model constants.The present results show that the error can be reduced dramatically by using the dynamic models (Fig 5 Right and Table 1).The use of van Driest damping functions reduces the error to a comparable level to that obtained by using dynamic models.Solving Reynolds Average NSE near the wall and LES in the bulk of the flow by using SA model predict superior results to those obtained with Smag and One models but just like the other tested models in this study it could not predict the separation of the flow.The Deard model shows less sensitivity to mesh resolutions than the Smag and One models and it predicts the same friction velocity as in the highly-resolved case.The error analyses (not shown here) reveal that both Smag and One models produce almost the same high Ξ (U) regardless the streamwise position, while the other models exhibit an enormous spatial variation with largest value of Ξ (U) at (x/λ = 0.2) and minimum between (x/λ = 0.5 and 0.8).2,…,0.9,Left: by using by using (Smag, Smag+LW1, Smag+LW2) SGS models, right: by using by using (dynLag, dynLag+LW1, dynSmag+LW2 SGS models.

Wall model influence
Adding the local wall model in its two implementation variants (LW1 and LW2) to Smag model enhance the prediction of velocity profiles, see Fig 5. Vague differences between the smag+LW1 and smag+LW2 are noticed.However the use of smag+LW1 produces a small and weak separation zone while smag+LW2 model fails to predict any separation zone.The second order statistical moments of the resolved velocities are also improved by using the wall models (Table 1).The effects of localization of the wall models can be seen in the spatial variation of (k) and ‫‪ത‬ݓ́ݑ(‬ തതത ) in  ) However, the two wall models exhibit distinguishable differences when they are used with dynLag model.The dynLag+LW1 model predicts a stronger separation zone than the Smag+LW1 and just like Smag+LW2 the dynLag+LW2 does not predict any separation zone.Both wall model enhance the results significantly and reduce the differences in velocity profiles between LES and DNS data to the highly-resolved LES level (Table 1).Large differences between the wall models are found in (k) and ‫‪ത‬ݓ́ݑ(‬ തതത ) in Fig 7 .Modifying the SGS eddy viscosity to take into account the wall model in case of dynLag+LW2 produces a very high (k) near the wall and results in overestimation of (k) everywhere in the whole domain.dynLag+LW2 overestimates the vertical shear stress at all (x/λ) positions and shifted the profile peaks upwards while dynLag+LW1 underestimate a bit the ‫‪ത‬ݓ́ݑ(‬ തതത ) but produces a profiles that has the same shape as DNS.
The two other wall models that are tested with SA model show superior results to those obtained with SA with no-slip boundary condition but they fail to predict any separation zone.

Conclusions
Several SGS model are used to simulate a separated flow in a wavy channel.The highly-resolved LES predict correctly the flow features and compared well with DNS data of Maass and Schumann [28].The differences between SGS models are indiscernible at highly-resolved LES.However, by using coarser mesh the differences between SGS model become distinguishable.At coarse mesh without wall models, no model could predict the separation zones.The Deardorff stress SGS model exhibits superior results to those obtained by using Smagorinsky or One-Equation SGS models.The present results revel that adding dynamic procedures to the models reduces the model differences between the LES and the DNS results to almost half.The study has also demonstrated that the use of wall model based on logarithmic profile that is held locally in space and instantaneously in time with course meshed LES improves the results and returns the closest fidelity to the highly-resolved LES.However, these models show sensitivity to the implementation and depending on the way they are implemented they can predict the separation zone and the gross flow features with different accuracies.

Figure 2 :
Figure 2: contours of LES phase averaged streamwise velocity obtained by using Smagorinsky SGS model.

2nd
Symposium on OpenFOAM in Wind Energy 03002-p.5 The phase average (the averaging over y-direction and time, with an additional averaging over all identical 6 waves) of streamwise velocities are shown in Fig 2. The flow is characterized by the Reynolds number and wave steepness.Waves with such relatively large steepness add a degree of complexity to the flow by causing streamline expansion in the wave leeward side.The flow decelerates and the adverse pressure is sufficiently high in this case to cause flow separation.A thin shear layer develops between the outer flow and the separated flow.Further downstream, the streamlines are curved inward and the flow reattaches to the surface forming a bubble of separated flow.A thin boundary layer develops downstream the reattachment point on the windward side of the next wave crest.Fig 2 displays that the flow perturbation caused by the undulated surface, which its amplitude is about 10 % of the boundary layer height (the half of the channel height), influences the overall flow behavior.The flow exhibits observable spatial variation with flow separation, reattachment and shear layer.These flow features are often encountered in wind energy applications.

Figure 3 .
Figure 3.Comparison of well resolved LES by using (Smag, One, Deard) SGS models and DNS by means of vertical flow profiles of phase averaged streamwise velocity normalized by the bulk velocity at positions (x/λ)=0.1,0.2,…,0.9.The results of highly-resolved LES exhibit a good agreement with DNS data of Maass and Schumann [28].The phase averaged velocity profiles at different (x/λ) positions obtained by three SGS models (Smag, One, and Deard) are compared to DNS results in Fig 3.The relative differences between the LES and DNS data of phase averaged streamwise velocity, calculated by using L 2 norm (Ξ = ඥ∑ |ாௌିேௌ| మ ඥ∑ |ேௌ| మ
).This high constant value causes excessive damping of large-scale fluctuations because of the high shear near the wall, manifested in the low values of turbulent kinetic energy and vertical shear stress predicted by standard Smag in the left part of Fig 5.
Fig 6 comparing with the use of Smag model without any wall treatment that has almost the same profiles and values regardless the near wall flow topology.Near the wavy surface the local wall models have the similar trends as DNS kinetic energy profile but with lower values and with profile peak shifted upward due to the vertical resolution and the absence of the large separation zone.