Large Eddy Simulation with Energy-Conserving Schemes and the Smagorinsky Model: A Note on Accuracy and Computational Efficiency ^†

Abstract

Despite advances in turbulence modelling, the Smagorinsky model remains a popular choice for large eddy simulation (LES) due to its simplicity and ease of use. The dissipation in turbulence energy that the model introduces, is proportional to the Smagorinsky constant, of which many different values have been proposed. These values have been derived for certain simulated test-cases while using a specific set of numerical schemes, to obtain the correct dissipation in energy simply because an incorrect value of the Smagorinsky constant would lead to an incorrect dissipation. However, it is important to bear in mind that numerical codes may suffer from numerical or artificial dissipation, which occurs spuriously through a combination of spatio-temporal and iterative errors. The latter can be controlled through more iterations, the former however, depends on the grid resolution and the time step. Recent research suggests that a complete energy-conserving (EC) spatio-temporal discretisation guarantees zero numerical dissipation for any grid resolution and time step. Therefore, using an EC scheme would ensure that dissipation occurs primarily through the Smagorinsky model (and errors in its implementation) than through the discretisation of the Navier-Stokes (NS) equations. To evaluate the efficacy of these schemes for engineering applications, the article first discusses the use of an EC temporal discretisation as regards to accuracy and computational effort, to ascertain whether EC time advancement is advantageous or not. It was noticed that a simple non-EC explicit method with a smaller time step not only reduces the numerical dissipation to an acceptable level but is computationally cheaper than an implicit-EC scheme for wide range of time steps. Secondly, in terms of spatial discretisation on uniform grids (popular in LES), a simple central-difference scheme is as accurate as an EC spatial discretisation. Finally, following the removal of numerical dissipation with any of the methods mentioned above, one is able to choose a Smagorinsky constant that is nearly independent of the grid resolution (within realistic bounds, for OpenFOAM and an in-house code). This article provides impetus to the efficient use of the Smagorinsky model for LES in fields such as wind farm aerodynamics and atmospheric simulations, instead of more comprehensive and computationally demanding turbulence models.

1. Introduction

Large eddy simulation (LES) is a fine alternative for simulating high Reynolds number flows as compared to the computationally demanding direct numerical simulation (DNS). LES is tractable because it only resolves the large, energy-containing scales of a flow, while modelling the smaller and more numerous dissipative scales mathematically with a subgrid scale (SGS) model. Given the numerical complexity of this process, care must be taken in terms of detecting and handling errors stemming from spatial discretisation and the SGS model. A first insight into these errors was provided by [1], who demonstrated that the errors from spatial discretisation and turbulence modelling are equally consequential; which were later discovered to be capable of cancelling each other in certain cases, leading to a desirable result [2,3,4]. It was also noted that approaches designed to suppress the two sources of errors simultaneously, may not necessarily lead to a low overall error. References [4,5] illustrated how the same SGS model could produce different results under the influence of different spatial discretisation schemes. The former also commented on the dependence of the Smagorinsky constant on the grid resolution [6]. Further, the accuracy of LES can also be affected by time-integration, although this can be corrected easily by reducing the time step or by using higher-order time-integration schemes.

Most of the quoted studies used a finite volume method (FVM) for solving the Navier-Stokes (NS) equations. These were combinations of 2nd and 4th order discretisation of the convective and diffusive terms of the NS equations, advanced in time with a 2nd order explicit Runge-Kutta scheme. Further, they commented mostly on the effect of spatial discretisation and the resulting numerical dissipation on the overall accuracy of LES. As an adjunct to these studies, this article discusses the use of energy-conserving (EC) spatio-temporal discretisation schemes for FVMs. EC schemes are designed to conserve the kinetic energy (KE) of an inviscid flow in the absence of boundaries and body forces. As a result, for a viscous flow, EC schemes ensure that the loss of KE is brought about only through the action of molecular viscosity, which is an invariant property of incompressible flows [7]. An EC spatial discretisation precludes numerical dissipation that may increase (or even reduce) the dissipation of the turbulence kinetic energy (TKE). However, concomitantly, an implicit EC time integration scheme is necessary to ensure that the conservation of KE remains enforced with time-stepping because numerical errors through time integration could also bring about numerical dissipation.

In keeping with above, this article first addresses whether an implicit EC time integration scheme is crucial for averting numerical dissipation or could one use an explicit non-EC scheme and yet reduce the numerical dissipation to a value that the resulting solution is accurate? This is solely because implicit EC schemes are computationally demanding and numerically complex to implement, as a result of which, they may not be most suited for engineering applications that concern the research described in this article, namely, the simulation of the atmospheric boundary layer (ABL) and wind farm aerodynamics.

Secondly, the article compares an EC spatial discretisation with the simple central difference scheme. An EC spatial scheme is based on arithmetic means but on a staggered grid. On the other hand, a central difference scheme is based on arithmetic means on a collocated grid. As a result, both schemes do not have a dissipative error per se. This article considers only Cartesian grids on which the theoretical advantages of EC schemes have been established. but have a non-zero dispersive error. This article assesses whether there is a difference between the EC and central difference approaches (on their respective grids) in terms of how they affect the outcome in LES.

The tests described in this article are simulated using the simple central difference scheme on collocated grids available in the open source field operation and manipulation (OpenFOAM) toolbox [8] and the EC schemes within the Energy-Conserving Navier-Stokes (ECNS) code [9] owned by the Energy research Centre of the Netherlands (ECN, now TNO).

To avoid confusion, the term dissipation-free instead of energy-conserving will be used to indicate the property of conservation of energy obeyed by at least, the EC scheme and must be obeyed by any dissipation-free scheme in general. Therefore, the abbreviation EC, henceforth strictly refers to the numerical scheme offered within the ECNS and not the property of being dissipation-free. Further, the term EC-discretisation refers to discretisation of the convective operator in the NS equations.

Table 1. Salient features of the two dissipation-free spatial discretisation methods studied in this article.

Scheme	ECNS	OpenFOAM
Grid	Staggered	Collocated
Interpolation	Arithmetic mean	Linear
Differencing	Central	Central

lists the features of the ECNS and OpenFOAM in terms of spatial discretisation. The only aspect in which these two FVM-based codes differ, is their formulation on non-uniform grids. The EC scheme remains dissipation-free as it uses the arithmetic mean to interpolate and average the velocity field across cell boundaries [7]. Whereas, the central difference scheme in OpenFOAM switches to linear interpolation on non-uniform grids. This article assesses if this difference makes one scheme more favourable than the other and under what circumstances. Finally, the article discusses whether the absence of numerical dissipation promotes the use of the simple Smagorinsky model for the LES of turbulent flows with FVMs.

It is important to emphasise that this article is aimed solely at assessing how effectively the Smagorinsky model can be adapted to simulate flows with the EC schemes, towards ascertaining if such an approach could indeed find its way into practical applications in engineering.

2. LES and the Smagorinsky Model

The first step in LES involves modifying the NS equations to represent only the large scales, specifically, the scales greater than a cut-off length, $\Delta$ [10,11]. This operation that is formally defined as filtering, when applied to the incompressible form of the NS equations leads to,

$$\frac{\partial \widetilde{u_i}}{\partial x_i}=0,$$

(1)

$$\frac{\partial {\tilde{u}}_i}{\partial t}+\frac{\partial \widetilde{u_i{u}_j}}{\partial x_j}=-\frac{\partial \tilde{p}}{\partial x_i}+\frac{\partial }{\partial x_j}\left(\nu \left(\frac{\partial {\tilde{u}}_i}{\partial x_j}+\frac{\partial {\tilde{u}}_j}{\partial x_i}\right)\right),$$

(2)

where the pressure p has been divided by the constant density of the fluid and $\nu$ is the kinematic viscosity.

As opposed to the standard representation of the component of velocity in the i-th direction as $u_i$, ${\tilde{u}}_i$ represents the velocity field of the scales greater than $\Delta$. In practice, $\Delta$ is generally equal to the grid size and hence, all scales that are smaller than $\Delta$ are called the subgrid scales. Next, the filtered convective term is further simplified using Leonard’s decomposition [12] or,

$${\tau }_{ij}=\widetilde{u_i{u}_j}-\widetilde{u_i}\widetilde{u_j},$$

(3)

which leads to,

$$\frac{\partial {\tilde{u}}_i}{\partial t}+\frac{\partial {\tilde{u}}_i{\tilde{u}}_j}{\partial x_j}=-\frac{\partial \tilde{p}}{\partial x_i}+\frac{\partial }{\partial x_j}\left(\nu \left(\frac{\partial {\tilde{u}}_i}{\partial x_j}+\frac{\partial {\tilde{u}}_j}{\partial x_i}\right)\right)-\frac{\partial {\tau }_{ij}}{\partial x_j},$$

(4)

where ${\tau }_{ij}$ is the SGS tensor that represents the interaction between the subgrid scales and the large resolved scales. The SGS tensor is generally represented by SGS models that are semi-empirical. For more details, the readers are referred to two important reviews on this subject by [13,14]. The Smagorinsky model that is used in this article is based on expressing ${\tau }_{ij}$ in terms of the resolved velocity field and eddy viscosity as [6],

$${\tau }_{ij}-\frac{1}{3}{\delta }_{ij}{\tau }_{kk}=-2{\nu }_T{\tilde{S}}_{ij},$$

(5)

where ${\tilde{S}}_{ij}$ is tensor of the resolved strain rate or the strain rate of the resolved velocity field. ${\nu }_T$ is the eddy viscosity, which in terms of the resolved velocity field is defined as,

$${\nu }_T={\left(C_S\tilde{\Delta }\right)}^2\left(2\right|{\tilde{S}}_{ij}{{|}^2)}^{\frac{1}{2}},$$

(6)

where $C_S$ is the flow-specific Smagorinsky constant. The values reported for $C_S$ range from 0.1 to 0.23, of which most have been deduced from experimental or DNS data (see [15]).

In practice, LES codes report different values of $C_S$ purely based on the study of decaying turbulence alone. Reference [16] proposed a value of 0.17. Reference [17] obtained a value of 0.144 with the EllipSys3D LES code [18,19], whereas [5] used 0.15. This article will also report different values of $C_S$ for the ECNS and OpenFOAM.

These different values stem from different discretisation schemes and their errors that are always present unless the grid is very fine and the time step very small, all of which affect how the Smagorinsky model must be tuned to obtain the correct overall dissipation. It is important to note that the Smagorinsky model does not involve explicit filtering of the velocity field. The grid itself acts as the filter that only resolves all scales bigger than the grid size. Further, $\tilde{\Delta }$ is the cut-off scale in space that in practice should be set equal to ${\Delta }_g$, the overall grid size, which is ${\Delta }_g=\sqrt[3]{{\displaystyle \prod _{i=1}^3}{\Delta }_i}$, where ${\Delta }_i$ is the average grid size along the i-th direction. One may choose a larger value of $\tilde{\Delta }$ in which case, the Smagorinsky model only parameterises the scales smaller than $\tilde{\Delta }$ that must then be referred to as the subfilter scales (SFS). In this article, ${\Delta }_i$ is fixed as the overall grid size. More details on the relationship between subgrid and subfilter scales can be found in [3,20].

3. Dissipation in a Numerical Scheme

Within the discrete NS-LES system, the errors that arise can be classified as those from the regular operators of the discrete NS equations or those from the SGS model. The former comprises errors from the discretisation of the convective and diffusive terms, while the latter comprises the modelling error due to the nature of the SGS model and the discretisation error through the implementation of the SGS model in the NS solver.

One can roughly express the overall dissipation in an LES scheme as,

$$\mathcal{P}={\epsilon }_C(\Delta x)+\nu (\mathcal{D}+{\epsilon }_D(\Delta x))+\mathcal{S}\left(\tau \right)+{\epsilon }_S(\Delta x)+{\epsilon }_T(\Delta t).$$

(7)

$\mathcal{D}$ and $\mathcal{S}$ represent the diffusive (physical diffusion) and the SGS operator, respectively. ${\epsilon }_C$ and ${\epsilon }_D$ are the discretisation errors (dissipative) accompanying the convective and diffusive operators, respectively. ${\epsilon }_D$ (dissipative) arises from the approximation used to estimate the physical diffusion through molecular viscosity. Dissipation through the SGS model is represented by a function of the model $\tau$, which can be tuned (see Section 2) and the spatial discretisation error in the model’s implementation ${\epsilon }_S$, which is a function of the grid resolution $\Delta x$. The pressure term also leads to numerical dissipation but not with an EC scheme, this dissipation will be ignored for now.

${\epsilon }_T$ is the dissipation through time integration, which is a function of the time step and must approach zero with reducing time step. It is also known from [21] that numerical dissipation through time integration can alter how the flow evolves and hence, must also be accounted for in the above equation.

Physical dissipation or $\mathcal{P}$, pertains to the dissipation that one would try to achieve with a numerical method, which should equal what the flow would experience in nature. In a high Reynolds number flow, the molecular diffusion is very small as compared to the turbulent diffusion and so is the error from the discrete diffusive operator, both of which can be neglected. In LES of wind farms and the atmospheric boundary layer, the diffusive term is set to zero. However, the diffusive term has not been neglected for the low Reynolds number tests reported in this article. Assuming that the error from time-integration is reduced with a small time step or an EC time integration scheme, Equation (7) then reads,

$$\mathcal{P}={\epsilon }_C(\Delta x)+\mathcal{S}\left(\tau \right)+{\epsilon }_S(\Delta x).$$

(8)

In a dissipation-free spatial scheme, the error from the convective operator is predominantly dispersive (see [21]). For a high Reynolds number flow, the convective term is large as compared to the term representing molecular diffusion. It is therefore likely that numerical dissipation from the discrete convective term could be comparable the term representing the discrete molecular diffusion. In doing so, the numerical dissipation through ${\epsilon }_C$ can influence the flow to a great extent and should be removed. As a consequence with an LES code based on a dissipation-free scheme, one would expect Equation (8) to read,

$${\mathcal{P}}_{EC}\approx \mathcal{S}\left(\tau \right)+{\epsilon }_S(\Delta x).$$

(9)

Whereas, with a non-EC scheme one obtains,

$${\mathcal{P}}_{Non-EC}\approx {\epsilon }_C(\Delta x)+\mathcal{S}\left(\tau \right)+{\epsilon }_S(\Delta x).$$

(10)

Next, the dissipative errors from the implementation of the SGS model and ${\epsilon }_S(\Delta x)$ are combined with the dissipation arising from the model itself i.e., $\mathcal{S}\left(\tau \right)$, to rewrite the above equations as,

$${\mathcal{P}}_{EC}\approx \mathcal{S}(\tau ,\Delta x),$$

(11)

$${\mathcal{P}}_{Non-EC}\approx {\epsilon }_C(\Delta x)+\mathcal{S}(\tau ,\Delta x).$$

(12)

It is hard to assess the contribution of $\Delta x$ within $\mathcal{S}(\tau ,\Delta x)$ but it is nonetheless apparent that the contribution of $\Delta x$ as a whole to $\mathcal{P}$, reduces with the removal of ${\epsilon }_C(\Delta x)$. Therefore, it is fair to hypothesise that with a non-dissipative scheme (EC or pseudo-spectral), the dissipation is unaffected by the grid resolution to a greater extent that what would be the case with a scheme that has inherent numerical dissipation.

This property may help tune the SGS model more easily as the grid resolution has a relatively minor role to play, when ${\epsilon }_C(\Delta x)$ is absent or very small compared with $\mathcal{S}\left(\tau \right)$ (if some numerical dissipation is present). In fact, the property of inherent numerical dissipation in terms of ${\epsilon }_C(\Delta x)$ is put to good use by Implicit LES methods wherein, the dissipation introduced by the discretisation replaces the SGS model (see [22]).

Although one cannot completely rule out the dependence of the SGS model on grid resolution (or ${\epsilon }_S(\Delta x)$), evidence from literature suggests that within a reasonable range of grid resolutions, the dependence is moderate enough for practical applications such as the simulation of high Reynolds number flows within the ABL (confirmed for the Smagorinsky model in [23]). It has been observed that at poor grid resolutions, the optimal Smagorinsky constant is quite dependent on the grid resolution. However, once the grid resolution is sufficient to promote adequate resolution of the large scales relevant for high Reynolds number flows, the dependence on increasing grid resolution is relatively weaker.

4. ECNS and OpenFOAM

This section provides information on the two codes considered in this article, ECNS and OpenFOAM; and the relevant spatial discretisation and time integration schemes.

4.1. ECNS

The ECNS offers 2nd order spatial discretisation on a staggered and Cartesian grid. For time integration, one can choose from the classical explicit 4th order but non-EC Runge-Kutta scheme (Ex4) and the implicit 2nd and 4th order EC Runge-Kutta schemes based on the Gauss quadrature, henceforth referred to as the Im2 and Im4. These time integration schemes are summarised in Appendix A; one can find the details in [21].

4.2. OpenFOAM

OpenFOAM (http://www.openfoam.com) is an open-source simulation package that permits the user to customise an NS solver using a combination of spatial-temporal schemes and a turbulence model. It uses a finite volume approach like the ECNS.

Amongst the host of schemes offered, the 2nd order central difference-based Gauss Linear scheme has inherently zero numerical dissipation but a non-zero dispersive error [24]. In addition, the inherently dissipative quadratic upwind interpolation scheme for convective kinematics (3rd order, QUICK) and the linear upwind (2nd order, LU) schemes are also considered. As regards to Table 1, QUICK uses a collocated uniform grid but with upwind differencing and interpolations defined by the QUICK scheme itself. For time integration, the Crank-Nicolson 2nd order implicit forward time stepping is used with a very small time step to avert numerical dissipation through time integration [24].

It is important to mention that improper coupling of pressure and velocity within a solver can affect the KE. Appendix A provides the relevant details. For simplicity, the spatial discretisation approach with the EC scheme on staggered grids is directly referred to as ECNS and the central difference scheme on collocated grids as OpenFOAM.

5. Tests

The first objective is to determine what time integration schemes require the least amount of computational effort while guaranteeing that the numerical dissipation through time integration does not affect the flow.

Secondly, this article establishes qualitatively whether Equations (11) and (12) hold for the spatial discretisation approaches mentioned in Table 1. Next, the consequences that result should Equations (11) and (12) not hold, would be discussed followed by an assessment of whether the differences between a central difference scheme on a collocated grid and an EC scheme on a staggered grid, lead to marked differences the results so-obtained.

5.1. Grid-Generated Turbulence

The test case simulated in this article is the recognised experimental campaign by [25] (CBC), who used a grid to generate an isotropic homogeneous turbulent field in a wind tunnel. The experiments provided details on the energy spectrum and the statistical properties of the turbulent flow at three instances of time, as the flow moves through a wind tunnel at a near constant average velocity, $U_o=10$ ${\mathrm{ms}}^{-1}$. This case was chosen as it has been widely used for the validation of LES codes and also for the tuning of the Smagorinsky model [17].

One must bear in mind that the flow’s Reynolds number is low, around 34,000 based on the average velocity and the size of the grid used to create turbulence ($M_g=5.08$ cm) and that the flow is not truly isotropic with an isotropy ratio of 0.95 (the ratio of the standard deviation in streamwise velocity to that in transverse velocity). Hence, the TKE is higher in the transverse direction than in the streamwise direction.

Although the slight anisotropy can be overlooked, the low Reynolds number requires attention given that the fundamentals of LES, by definition, consider a high Reynolds number flow with an extended inertial range [10]. Further, the experimental energy spectra in the inertial range deviate a little from the $-5/3$ $Kolmogorovspectrum$ that is observed in turbulent flows [25,26,27] (see Figure 1). These cannot be verified even through DNS because of the wide range of length and time scales in the flow, for which the current computational resources are insufficient.

Nevertheless, this article provides an assessment of how closely the simulations can replicate the spectra from the largest to the smallest scales on the grid, which indirectly should also estimate the correct decay of TKE with time. It is fair to mention that the intention is not to validate the ECNS or perform a follow-up study of the experiments by [25]. Other relevant test cases from literature were also considered; Appendix B provides the details. Apart the cases mentioned in Appendix B, the ECNS was also tested against the roll-up of a shear layer, the simulation of the ABL and the experimental analysis of a single wind turbine in turbulent flow, while being compared against a range of LES codes available in literature. More details on these simulations can be found in [28,29].

5.2. Effect of Time Step and Time Integration

One notices three experimentally obtained spectra in Figure 1 (CBC). These correspond to three different locations in a wind tunnel at which the flow was analysed. The codes are used to simulate the flow between Stages I and III, referred to as $U_{o}t/M_g$ = 42 and 171, respectively, in [25] with an intermediate stage, $U_{o}t/M_g$ = 98 or Stage II, where t is the time. These numbers refer to distances downstream of the grid that generates the turbulence. The generation of the initial turbulent velocity field using the experimental data is explained in Appendix C.

The CBC test case is simulated on ${32}^3$, ${48}^3$ and ${64}^3$ uniform 3D grids with average grid sizes (${\Delta }_g=\sqrt[3]{{\displaystyle \prod _{i=1}^3}{\Delta }_i}$) equalling 1.72 cm, 1.14 cm and 0.86 cm, respectively. These sizes are small as compared to the longitudinal length scale of 2.4 cm but larger than the Taylor micro scale of 0.48 cm at Stage I, which places the cut-off filter (grid-based) somewhere in the inertial range, as the case should be. Therefore, a ${16}^3$ grid would be quite coarse. In certain cases, ${80}^3$ and ${96}^3$ grids are also used. With regards to the time-stepping, the arguments are centred on a time step of $\Delta t=$ 0.0001 s. This time step was chosen because it was used as reference by [17] to ensure that the Courant number was below 0.1 for a second order implicit time integration scheme, while tuning an LES code with the CBC case for two different grid resolutions.

Further, the time integration schemes are operated with a range of pre-determined time steps. When a time integration scheme is used with $\Delta t=$ 0.001 s, it will be quoted by its name alone, i.e., as Ex4 or Im2. For other time steps, an extra letter is added at the end of the scheme’s name, for example, Ex4Q. The time steps and their corresponding letters are summarised in

Table 2. A summary of the time steps used.

Notation	Time Step
Ex4Q	0.004 s
Ex4D	0.002 s
Ex4	0.001 s
Ex4h	0.0005 s
Ex4q	0.00025 s
Ex4te	0.0001 s

with Ex4 to serve as an example.

For most simulations with the Smagorinsky model, the $C_S$ is set to 0.1 [30], which will later be determined more accurately for the ECNS. As an output, the energy spectrum, the rate of decay of TKE and the computational time are monitored. Based on how accurately the spectra and decay of TKE are predicted and with what computational effort, this article will outline which time integration scheme is most appropriate for the ECNS. Using the ${32}^3$ grid, the CBC experiment is simulated with Ex4, Im2 and Im4 (the three schemes available within the ECNS), in combination with different time steps, the acronyms for which are defined in Table 2. The reason behind choosing the ${32}^3$ grid apart from being computationally efficient, is that the energy spectra from CBC suggest that the inertial range is not that wide; contrary to what is expected of a high Reynolds number flow. Although LES can resolve the large energy-containing range and a part of the inertial range but one must avoid resolving the lower inertial range that is comparable with the Kolmogorov scales, most of which are represented by the SGS model. Ideally, the cut-off scale should be in the inertial range such that it is comparably smaller (say a factor of 10 to 100) than the large scales but larger by a similar factor with respect to the Kolmogorov scales (see [10] for a discussion).

5.2.1. Stability

The stability of the Ex4 is conditionally determined by the Courant Friedrichs Lewy (CFL) criterion for explicit time integration schemes. Therefore, at time steps larger than 0.004 s, the simulation was unstable with a Courant number greater than ∼2.4. Reducing the time step to obtain a Courant number below 1, for example $t=$ 0.001 s, guarantees stability. On the other hand, owing to its unconditional stability, the Im2 is stable even at a time step of 0.01 s that is ten times the reference time step, as opposed to the Ex4 that is divergent [21].

Now to compare the three methods, one must ensure that they all are stable at the time step considered. Therefore, the Ex4 would be used only with time steps smaller than 0.002 s.

5.2.2. Error Analysis

First, a reference solution is calculated using Ex4 with a very small time step of 0.00005 s, which was sufficiently low to ensure an error below $\left({10}^{-9}\right)$ through time integration, for all the time integration schemes considered here. Next, the root-mean-square of the difference between the TKE of the field obtained during a simulation and the TKE of the reference solution (interpolated in time for the relevant time step), is analysed over the duration of the simulation. The TKE is defined as,

$$\mathrm{TKE}=\frac{1}{2}\left(\sum _{i=1}^3\overline{u_i^{{}^{\prime }}{u}_i^{{}^{\prime }}}\right),$$

(13)

where $\overline{u_i^{{}^{\prime }}{u}_i^{{}^{\prime }}}$ is the standard deviation in the component of velocity along the i-th direction. This value is in fact, the average KE associated with the fluctuations in velocity across the range of eddies within the flow (presented as energy per unit mass). The TKE is chosen instead of the KE as the changes in KE in the CBC case are hard to detect due to the high average streamwise velocity of 10 ${\mathrm{ms}}^{-1}$ as compared to velocity fluctuations that are only 2.22% of the streamwise value (see [25] for experimental values). The error is calculated using the following definitions of $e_2$ and $e_{\infty }$,

$$e_2={\left(\frac{1}{n}\sum _{i=1}^n{\left|{\epsilon }_i\right|}^2\right)}^{\frac{1}{2}},$$

(14)

$$e_{\infty }=\max \left|{\epsilon }_i\right|,$$

(15)

where $\epsilon$ is the difference between the TKE of the reference solution and the simulated TKE at a given instance of time and n is the number of such observations. Figure 2 shows the change in TKE with time as obtained with various time steps and time integration schemes. The errors are shown in Figure 3 with the computational time of the simulation. T is the time needed to simulate the test case on a ${32}^3$ grid using Ex4 and $\Delta t=$ 0.001 s.

5.2.3. Order

Figure 3 shows the 2nd order accuracy of the Im2 and the 4th order accuracy of the Ex4 and Im4 (design orders).

For a time step larger than 0.001 s, the Ex4 incurs more error as compared to the Im2. This is best explained in terms of the methods’ orders of accuracy. The Im2 converges with second order accuracy, hence the reduction (or increase) in error with the refinement of time step should be less than what is noticed with the fourth order Ex4. Therefore, Im2 is more accurate at larger time steps with a lower magnitude of error. But as the time step is reduced, it is the non-EC Ex4 method that is more accurate.

Im4 leads to the lowest global error at nearly all time steps. However, at very small time steps, one observes that Im4 incurs a larger error than Ex4 and loses its 4th order accuracy. One of the reasons behind this is that the iterative method being used to solve the Im4 system reaches its limit of accuracy as the tolerance cannot be reduced further due to computational restrictions. Further, at very large time steps (0.004 s), the accuracy is compromised because the time step is not sufficient to resolve the turbulent phenomena that may have a time scale that is smaller than the time step. Figure 4 presents a similar analysis on ${48}^3$ and ${64}^3$ grids. One notices slopes corresponding to the accuracy of the relevant method, independent of the grid resolution. Further, the error is larger with refined grids at a given time step, as one must ideally reduce the time step to maintain the same Courant number.

5.2.4. Accuracy

To evaluate the accuracy of the time integration schemes, one must have the correct value of $C_S$. Therefore, this article evaluates the effect of changing the time step and the time-integration scheme on a reference solution calculated with a fixed value of $C_S$. These changes are evaluated in terms of the energy spectrum, which is the distribution of the TKE across the various scales of turbulence. This will be helpful in the next section that outlines the tuning of the Smagorinsky model. If the model is tuned on the basis of the TKE, one would ensure that the total TKE is correctly predicted. But tuning with respect to the spectra will not only ensure a correct value of TKE but also its distribution within the various scales of turbulence.

In Figure 5 (compare the energy spectra at Stages II and III), one notices that a time step of 0.002 s with the Ex4 results in a large amount of numerical dissipation as compared to 0.001 s, causing a marked change in the energy spectra; and the decay of TKE, as shown previously in Figure 2. The plot also presents the ratio of the energy at a given wavenumber for a certain time step, to the energy at the same wavenumber obtained with the reference solution of Ex4 with a time step of 0.00005 s; mentioned as $A_{Tw}$ in the bottom-right plot within Figure 5. This helps visualise the distribution of the error in predicting the energy across wavenumbers. A value close to 1 indicates an accurate prediction, whereas values less than or greater than 1 indicate an under prediction or an over prediction, respectively. One notices that reducing the time step below 0.001 s does not lead to marked changes in the spectra ($A_{Tw}$ is nearly 1).

Comparing the error analysis (see Figure 4) with the change in TKE with time shown in Figure 2, one can establish that an error below ${10}^{-4}$ pertains to a solution obtained without significant influence of numerical dissipation through time integration. With the Ex4, this time step is below 0.005 s or a Courant number of less than 0.15, indicating that even with a higher order method may need a small time step to avert numerical dissipation.

On the other hand, Im2 requires a time step of 0.0025 s due to its lower order of accuracy, to obtain a final error below ${10}^{-4}$. But the Im4 being 4th order accurate and EC by design, is accurate enough in this case at a time step of 0.002 s itself (in fact, even at a time step of 0.002 s), as shown in Figure 6. As a result, one also notices that the spectra obtained with Im4 do not show much apparent variation in comparison with the spectra obtained with the reference solution of Ex4 with 0.00005 s. The $A_{Tw}$ remains nearly 1 at almost all the wavenumbers (see bottom-right plot in Figure 6).

One could conclude that Im4 is highly accurate as compared to the other two methods. However, it is important to consider the computational costs that one must reckon with while using higher-order implicit EC methods.

5.2.5. Efficiency

From the preceding discussion on accuracy, it is clear that Ex4 with a time step of 0.0005 s on a ${32}^3$ grid, leads to an accurate solution that is independent of time step, with a computational time of $1.3T$, where T is the time required by Ex4 (see Figure 3). This corresponds to an error below ${10}^{-4}$. The computational times on finer grids show a similar trend with the various schemes and time steps. These are not discussed for succinctness.

Upon reducing the time step, the error is driven below ${10}^{-8}$ with Ex4 using a time step of 0.0001 s, costing about $6.7T$ in computational time due to the method’s high order of accuracy. On the other hand, the Im4 permits using a larger time step of 0.002 s but with a very high computational time of $8.8T$. This computational time is very high to the extent that Ex4 with a time step of 0.0001 s is as fast as Im4 with a time step of 0.004 s.

Hence, the advantage of an implicit method permitting a larger time step, is annulled by the increased computational time required per time step. Although the ECNS is not parallelised, it is safe to believe that the ratios of computational times with the implicit methods presented in Figure 3, would at best remain the same upon parallelisation, as compared to the Ex4. An improved computational efficiency with the implicit methods would require a more optimised approach to solving the coupled system of equations (see Appendix D for a detailed discussion).

5.3. Effect of Spatial Discretisation Schemes

The first step to a fair comparison of the ECNS and OpenFOAM with regards to LES, is to tune the Smagorinsky model for the two codes. As per the hypothesis and Equations (11) and (12), the absence of numerical dissipation may reduce the influence of grid resolution on the overall dissipation, which in the case of the Smagorinsky model can be controlled with the Smagorinsky constant or $C_S$ (this represents ${\epsilon }_S$ and $\mathcal{S}$ in Equations (11) and (12)). Therefore, the absence of numerical dissipation may ensure that $C_S$ is grid-independent to a greater extent with non-dissipative schemes as compared to a code that uses dissipative spatial discretisation.

Sensitivity of the Spectra to Changes in Grid Resolution

For the time being, assume $C_S=$ 0.12. It shall be shown later that $C_S=$ 0.12 or 0.125 is appropriate for the ECNS. This value is merely used to help the reader understand the discussion that follows. for the ECNS. Figure 7 shows the energy spectra at different grid resolutions. Further, the bottom-right plot in the Figure 7 shows the ratio of the simulated and the experimentally determined value of the energy as a function of the wavenumber; represented by the symbol $A_{CBC}$.

One notices that the predicted spectra follow the experimental curve till a certain wavenumber beyond which, the predicted spectra show a slight deviation (dampening in magnitude) followed by an abrupt and large reduction in energy just before the cut-off wavenumbers (these abrupt changes and large deviations would be called damped tails). It is observed that upon increasing the resolution beyond ${48}^3$, the abrupt dampening continues to occur but shifts to higher wavenumbers. However, the lower wavenumbers are accurately predicted and are affected to a minor extent by the increasing grid resolution. This can be seen through the plot of $A_{CBC}$ wherein the lines are clustered till a wavenumber of 100.

One may conclude that beyond ${96}^3$, changes in grid resolution will not alter the spectra for the large scales (up to a wavenumber of 100) because the contribution of $\Delta x$ to $\mathcal{S}(\tau ,\Delta x)$ in Equation (11) (or the discretisation error in the SGS model) will be reasonably small. This leads to a more or less constant overall dissipation, governed solely by Smagorinsky model and its constant $C_S=$ 0.12, in this case. Therefore, it is reasonable to assume that for refined grids, one could probably resort to a single value of $C_S$, to correctly predict the spectra up the wavenumber at which the spectra seem to dampen.

Additional tests were conducted using different SGS models and even without one, to find out what causes the excessive dampening at high wavenumbers near the limit of the grid’s resolution. Figure 8 shows the results obtained with a pseudo-spectral code with the Germano model. The Germano model was also used in combination with the ECNS but the results obtained were poor, especially near the higher wavenumbers without any changes in the tails. This supported a hypothesis that the tails may arise from the use of a FVM. Therefore, it was deemed sensible to exclude regions that show these damped tails from further analyses on the spectra. The reason behind their occurrence has little relevance to this article but could be worth investigating.

Next, an analysis similar to the one followed in [31] is preformed to tune the Smagorinsky model for the ECNS.

5.4. Tuning the Smagorinsky Model

To tune the Smagorinsky model in the ECNS, one must first calculate what part of the energy spectra (Stage II and Stage III) can be predicted accurately through the simulations. For example, one may check how closely can the solution on a ${64}^3$ report the spectra up to a wavenumber equal to the cut-off wavenumber on a ${32}^3$ grid and likewise for different resolutions.

To begin with, the Smagorinsky constant is set to $C_S=$ 0.12, on a ${32}^3$ grid. The results are compared with LES data from [26] in Figure 8, which appear to be quite similar, apart from the sudden dampening of energy beyond a wavenumber of $\kappa =$ 100. Next, the same value of $C_S=$ 0.12 is used (and other values close to 0.12, for example, 0.11, 0.115, 0.125 and 0.13) to obtain the energy spectra on grids up to ${96}^3$. In either case, it is ensured that the time integration introduces no numerical dissipation (see Section 5.2).

Next, one sets three limits $\mathcal{A},\mathcal{B}$ and $\mathcal{C}$, where $\mathcal{A}$ is only a part of the simulated energy spectrum on a ${32}^3$ grid, $\mathcal{B}$ a part of the simulated spectrum on ${48}^3$ and $\mathcal{C}$ for a ${64}^3$ grid. Figure 9 shows these three limits with grey lines. One notices that these limits are not equal to the cut-off wavenumber of the pertinent grid. This however ensures that the damped tail noticed so far with the spectra obtained using the ECNS, are removed from the analysis due to the reasons mentioned earlier.

To estimate the error incurred while predicting the part of the spectrum up to the cut-off wavenumber obtained on a ${48}^3$ grid or $\mathcal{B}$ but using a ${80}^3$ grid, the spectra obtained with the ${80}^3$ grid are integrated till the corresponding grey line that indicates $\mathcal{B}$ in Figure 9. Then the difference in energy with respect to the experimental solution integrated till the wavenumber corresponding to $\mathcal{B}$ is calculated. This is repeated for all logical combinations (an illogical combination would be trying to assess the prediction on $\mathcal{C}$ with a ${48}^3$ grid) and the resulting errors are summarised in percentage within

Table 3. A summary of errors (reported in percentage) using the energy spectra obtained by CBC [25].

		II (%)					III (%)
${\mathit{C}}_{\mathit{S}}$	Case	32	48	64	80	96	32	48	64	80	96
0.115	$\mathcal{A}$	15	19	20	13.4	7.8	10.9	18.8	20.3	15.6	9.4
	$\mathcal{B}$		12.2	14.9	12.6	8.1		7.5	10	7.5	3.9
	$\mathcal{C}$			10.8	10.1	6.5			3	1.1	1.1
0.12	$\mathcal{A}$	11.2	15.5	18.1	12.9	7.8	6.3	14.1	18.8	14	9.3
	$\mathcal{B}$		8.1	12.6	10.1	6.8		2.5	7.5	6.3	3.8
	$\mathcal{C}$			7.1	7.7	4.8			≈0	≈0	2.2
0.125	$\mathcal{A}$	6.9	12.9	16.3	12	6.9	≈0	10.9	15.6	14.1	9.4
	$\mathcal{B}$		4.1	9.5	8.8	6.1		1.3	5	3.8	2.5
	$\mathcal{C}$			3.6	5.4	3			3.3	2.2	3.3
0.13	$\mathcal{A}$	12.1	14.3	17.6	13	7.6	≈0	14.2	18.4	15.7	9.6
	$\mathcal{B}$		7.7	10.1	9.9	6.3		2.4	7.3	6.2	3.5
	$\mathcal{C}$			6.8	7.9	4.6			3.3	2.6	3.5

(see Figure 9 to interpret results).

It is important to note that the error is calculated as the difference between the experimentally determined (from [25]) and the simulated TKE and not as the root-mean-square error of the differences observed in the spectra at different wavenumbers. This is because the simulated spectra show two regions of marked deviations at smaller wavenumbers (indicated in Figure 9 as R1 and R2), which as it appears, do not change with grid resolution and could lead to a constant offset in the calculated errors. Further, these region are also noticed in Figure 8, when simulated with a pseudo-spectral code and the Germano dynamic model (indicated in the plot with Men).

In general, the errors are lower for Stage III than Stage II. With $C_S=$ 0.12 one notices a negligible error in predicting $\mathcal{C}$ for Stage III, especially with finer grids. But in this case, the corresponding errors for Stage II are higher. On the other hand, $C_S=$ 0.125 leads to a better result for Stage II than Stage III. Increasing the $C_S$ beyond 0.125 or reducing it to 0.115, leads to a higher error. Therefore, one may conclude that an optimal $C_S$ belongs to the range [0.12, 0.125], with regards to the test case. Also, as hypothesised through Equation (11), this range of values is quite independent of the grid resolution, provided it is enough to resolve the large scales correctly; possibly due the absence of numerical dissipation guaranteed by the ECNS approach (and the low time step).

Another point worth mentioning is the shape of the experimental spectrum at Stage II (see R1 in Figure 9). One notices a sharply placed data point at a wavenumber of ∼25 (indicated with a cyan line), as opposed to the smooth spectra at the other stages and the ones obtained through simulations. It is probable that the deviation of the simulated spectra from this data point, contributes to the error at Stage II to a great extent, given the relatively high energy content at that wavenumber.

OpenFOAM and the Influence of Numerical Dissipation

The procedure in Section 5.4 is repeated with OpenFOAM using the spatial discretisation mentioned in Table 1, to obtain a value of $C_S=$ 0.135 on a ${32}^3$ grid, as shown in Figure 10. It was ensured that the time integration does not introduce numerical dissipation by using a very small time step of 0.0001 s.

Upon increasing the grid resolution, one notices that the OpenFOAM approach of a central difference scheme on a collocated grid behaves similar to the ECNS approach of an EC scheme on a staggered grid. One also obtains the damped tails that were noticed with the ECNS, indicating that these are possibly inherent to an FVM using the Smagorinsky model. Further, the tails dampen above a certain wavenumber, below which, the scales (the large energy-containing ones) remain nearly unaffected with changes in grid resolution (see Figure 10). Although on a ${96}^3$ grid one does see some dampening of the large scales at Stage III, Stage II of the process remains nearly unaffected.

OpenFOAM is also used with the upwind schemes, QUICK and LU mentioned in Section 4. An immediate consequence is the strong effect of numerical dissipation leading to a damping of the spectra, as shown in Figure 11. As the grid resolution in increased (only shown with QUICK), it is observed that there is a reduction in dissipation as the large scales move towards their experimentally established values. However, the spectra as whole does not recover its correct shape indicating incorrect decay rates at multiple wave numbers. However, the 3rd order formulation of QUICK introduces relatively lesser dissipation than its 2nd order counterpart, LU (compare top-left and top-right in Figure 11).

In short, one is able to propose a suitable value of $C_S$ that allows finer grids to predict a fair regime of the large scales correctly with the ECNS and OpenFOAM. But on the other, one must not forget that the overall error, especially on the coarser grids, is a combination of the contributions from the diffusive and the SGS operators, which can only be studied separately and practically if the Reynolds number is high enough to neglect the diffusive operator.

Upon conducting relevant tests on the LES of the ABL and the wake of a model wind turbine in a wind tunnel; and comparing the results so-obtained against those generated with other LES codes, it was noticed that the values of $C_S$ deduced in this article are rather accurate for the simulation of high Reynolds number flows. Details on these tests and the results so obtained can be found in Mehta [28].

Out of the tests conducted, the simulation of the ABL is of much relevance to the wind energy community. LES has evolved as an important tool for the simulation of wind farm aerodynamics to study phenomena such as the effect of large wind farms on local weather, the meandering of a turbine’s wake, wind power generation, to name a few [32]. In either case, one requires an accurate simulation of the ABL to create a suitable velocity field which encompasses the turbulent conditions that are present in the ABL. Therefore, after obtaining a possibly appropriate value of $C_S$ for the ECNS, it was first assessed if the value was suitable in terms of simulating the correct velocity profile and turbulence spectra within a neutral ABL. The results from this test will be described briefly in the next section.

6. Simulation of the ABL with the ECNS

To standardise the ABL simulations, the domain is chosen to match the one described by Porté-Agel et al. [33] in the study of neutral ABLs, the dimensions of which are listed in Appendix E. Further, Appendix E contains important details on the velocity profile expected inside the simulated neutral ABL (log-law profile) and the simulation set-up. The results discussed in this section are obtained with the ECNS in combination with the Smagorinsky model to perform an LES.

It is important to note that the values for $C_S$ derived earlier, were done so for isotropic flows. However, the flow within an ABL is anisotropic near the surface, which in turn will affect the dissipation rate. Mason [34] proposed the following ad hoc function to reduce the value of the Smagorinsky constant near the surface in keeping with increased anisotropy and loss of scale invariance [33].

$$C_{S_y}={[C_{S_o}^{-n}+{\{\kappa +(\frac{y}{\Delta }+\frac{y_o}{\Delta }\left)\right\}}^{-n}]}^{-\frac{1}{n}},$$

(16)

where $C_{S_y}$ is the Smagorinsky constant at height y above the ground, $C_{S_o}$ is the Smagorinsky constant for homogeneous isotropic turbulence, $\Delta$ is the overall grid resolution and n is an integer of choice.

To tune the ECNS for ABL flows, the values of $C_S$ found earlier were used in combination with a value of 1 or 2 for n. The results obtained were analysed in terms of the average streamwise velocity profile that should be close to the log-law for the ABL (see Appendix E) and the energy spectra as a function of height. For clarity, only the results relevant to this article are discussed here; more details can be found in Mehta [28]. The combinations of $C_S$ and n discussed in this section are $C_S=$ 0.125, n = 1 and $C_S=$ 0.12, n = 2, which rendered the most accurate velocity profiles.

The tuning was performed on ${64}^3$ and ${96}^3$ grids. Figure 12 depicts how the velocity prediction improves with resolution, especially in the lower part of the ABL close to the ground, regardless of the tuning parameters. Outside of this region, above $y/L_y=$ 0.1, the tuning parameters do not affect the profile. On a ${96}^3$ grid, the combination $C_S$ = 0.125, n = 1, approximates the log-law closely till around 20 m above the ground, followed by a under estimated velocity regime. On the other hand, the combination $C_S=$ 0.120, n = 2, underestimates the velocity below $y/L_y=$ 0.1. Although increasing the grid resolution improves the prediction with both $C_S=$ 0.120, n = 2 and $C_S=$ 0.125, n = 1, the latter renders a more accurate average velocity profile.

For better insight into the nature of turbulence within the simulated ABL, one must look at the spectra at various heights. One expects the spectra to segregate into two regimes. The first with a slope of nearly $-1$ in the region near the ground and the second with a slope of $-5/3$, located away from the ground, as noticed by various studies (see Porté-Agel et al. [33] and references therein). Figure 13 shows that around 400 m above the ground, the spectra collapse onto the standard $-5/3$ slope for isotropic homogeneous turbulence, more or less, independently of the tuning parameters.

Further, the inertial range clearly exhibits a slope of $-1$ followed by the sharp dampening of the spectra at higher wavenumbers that was also noticed in the CBC simulations. This can once again be attributed to the FVM used in the ECNS. Therefore, in this case too, one must only take the spectra up to this sharp tail into account.

Figure 13 depicts that the spectra near the ground do not depend much on the combination of the $C_S$ and n. In fact, they collapse well onto the $-1$ slope, apart from the sharp dampening at higher wavenumbers. Thus, with nearly similar spectra for either combination, $C_S$ = 0.12, n = 2 and $C_S=$ 0.125, n = 1, the latter is chosen as it renders a more accurate velocity profile.

7. Conclusions

With regards to an implicit EC time integration scheme, one may conclude that it is capable of averting numerical dissipation at large time steps, while remaining stable and ultimately leading to an accurate solution. But in terms of the test case considered in this article, it appears that a non-EC explicit scheme is computationally efficient with small time steps, as long as it is of the same order as the implicit EC scheme (4th order in this case). Using a 2nd order explicit non-EC scheme may not be favourable as even an implicit EC 2nd order scheme requires a reasonably small time step for accuracy.

Further, it is fair to mention that the implicit methods may retain all their advantages in simulations that involve local grid refinement and the presence of walls, which may require a non-uniform grid to obtain accurate results. In such cases, the stability criteria may be very strict for non-EC explicit methods to resolve the refined regions, possibly rendering them unworkable. On the other hand, most LES studies on wind farm aerodynamics mentioned in literature, prefer to model the shear stress at the surface instead of refining the flow near it to reduce computational requirements, thus favouring the non-EC explicit methods. To harness the accuracy and stability of implicit-EC methods, a detailed study of the procedure used to solve their system of equations must be performed to assess if the computational time for implicit-EC methods could be reduced. With regards to spatial discretisation, it is noticed that the OpenFOAM approach leads to large scale behaviour that is similar to what one notices with the ECNS approach. No pronounced advantages of latter were found, which would make it more favourable for applications such as wind farm aerodynamics or atmospheric simulations. Further, the damped tail in the energy spectra, noticed with either approach can be attributed with confidence to the use of an FVM with the Smagorinsky model. This simple model does not account for changes in dissipation with scale, leading to excessive damping at smaller scales (as confirmed with two different codes). With either code, one obtains a value of the Smagorinsky constant that is nearly independent of the grid resolution (within reasonable bounds), as long as the numerical schemes are dissipation-free.

Following a series of tests, of which the simulation of a neutral ABL is presented here, one may conclude that an EC scheme on a staggered grid and the central difference scheme on a collocated grid could be similar in the simulation of high Reynolds number flows (on uniform grids) with the simple Smagorinsky model, which is the motivation behind developing the ECNS.