Evaluating large eddy simulation results based on error analysis

The quantiﬁcation of the prediction accuracy in large eddy simulations (LES) is very challenging due to various interacting errors associated with this approach. When dealing with errors in LES using implicit ﬁltering, numerical and modeling errors have drawn the interest of many researchers. Little attention has been paid to other sources of discrepancies between LES results and reference data, namely sampling errors, inﬂuence of the initial conditions, improper boundary conditions or uncertainties issuing from reference data. A framework of metrics that includes all these issues is addressed in the present paper to study subgrid-scale (SGS) models for LES and to quantify their prediction accuracy and computational costs. The method is applied to a simple wall-bounded turbulent ﬂow at moderate Reynolds number. It turns out from the results obtained with six commonly used SGS models that wall-adapting models (WALE and SIGMA) and localized dynamic models reproduce the physics of the ﬂow ﬁeld more faithfully, reveal a superior prediction accuracy and have a similar computational cost than models using van Driest wall damping. Especially at the viscous wall region ( r + < 50), wall-adapting and localized dynamic models are more accurate, reﬂecting the proper near wall behavior of such models. Relying on the analysis of sources of various errors, uncertainties in LES are estimated and systematically assessed, and their inﬂuence on simulation results is quantiﬁed. Finally, engineering estimations of the required averaging time to obtain basic estimates of statistical quantities with a predetermined degree of accuracy are suggested. mean and rms velocity proﬁles predicted by the LES models in comparison with the generated DNS dataset. All velocity proﬁles are normalized by the friction velocity of the DNS. Results are exclusively shown for the Smagorinsky model with van Driest wall damping, WALE Smagorinsky model, as representatives for different near wall treatments. Similar results are obtained using the other models.

Several aspects were investigated, which play an essential role in the interpretation of the quality and accuracy of results obtained by large eddy simulations, e.g., [6,8,18,30,39]. Regarding numerical and modeling errors, one approach, the so-called error landscapes, was recently suggested to assess LES while facilitating the formulation of new insights into the behavior of subgrid-scale (SGS) models, numerical discretization and their interaction [37]. This consists of a systematic variation of simulation parameters resulting in a database analysis that allows to gain a general overview of the error behavior. Furthermore, the error landscapes approach was combined with error estimators, in order to assess the quality and reliability of LES [31]. In fact, it is worth mentioning that with LES and implicit filtering, the model depends inherently on the grid size, thus linking physical modeling characteristics to numerical characteristics. These two kinds of errors were estimated separately in [30] by means of a systematic grid and model variation method based on variant of Richardson extrapolations. Besides these quantitative estimations, the LES index of quality [6] or other qualitative measures [6,10,17] can be used to determine whether a given numerical grid might be suitable to produce accurate results.
While such numerical and modeling errors have drawn the interest of many researchers, very little attention has been paid to other sources of uncertainties, such as insufficiently large computational domains, improper inflow conditions, finite initial transient or sampling errors interfering with results of LES [7,9,28,46]. Furthermore, reference data reported in the literature of the same configuration can differ significantly [63] from each other. Additionally, since in theory each LES with implicit filtering becomes as accurate as desired with increasing spatial resolution, a comparison of the modeling error and prediction accuracy alone appears not to be suitable to assess the capability of LES methods. In this respect, it is also important to estimate the computational cost of a model [51], in particular for daily engineering and practical applications while assessing LES methods with implicit filtering. However, a rigorous characterization of LES methods in terms of both accuracy and computational cost is rarely reported in the literature.
In the present study, an attempt is made to consider all these issues and thus to fill this gap in the literature in order to further improve the assessment capability and accuracy of LES results. The evaluation procedure includes an extensive error analysis and a rational characterization of SGS models in terms of prediction accuracy and computational cost. In addition, several sources of discrepancies between LES results and reference data, namely sampling errors, influence of the initial conditions, discrepancies in the boundary conditions and uncertainties issuing from the reference data, are estimated and systematically assessed, and their effects on simulation results are analyzed.
To demonstrate the feasibility of the present evaluation procedure, a turbulent pipe flow test case is selected, since it features a wall-bounded turbulent flow as encountered in many engineering applications, e.g., [21,33,66]. A moderate Reynolds number flow (Re τ = 180 based on the friction velocity) is chosen similar to that found in several DNS and experimental studies [1,12,13,15,68]. Therefore, it enables an assessment of variations in the available reference data as well. Since classical assumptions in LES models are only justifiable for high-Reynolds number turbulence and apart from walls, considering extreme departures from such assumptions might be of particular interest when evaluating LES result at low Reynolds numbers in wall-bounded flows. However, since the moderate Reynolds number selected in the present study is not fully representative for a large number of industrial-like applications, the outcome of the comparative study of different SGS models regarding prediction accuracy and computational costs might be different for higher Reynolds number flows. Nevertheless, the presented framework of metrics is definitely largely applicable for flows with higher Reynolds number and more complex geometries.
To avoid any error contribution of the near wall modeling, the near wall flow structures are fully resolved in the present study. Six commonly used LES subgrid-scale models, that are thought to provide proper near wall behavior, are selected. The first two models, the Smagorinsky [58] and the one-equation SGS kinetic energy model [69] with van Driest wall function, integrate a damping function following the much earlier near wall treatment similar to low Reynolds number turbulence models in the RANS context [62]. The third and fourth SGS models belong to the wall-adapting LES models generation, namely the wall-adapting linear eddy-viscosity model (WALE) [44] and the SIGMA model [43], which do not require any ad hoc treatment to reproduce the correct flow behavior in near wall regions. The last two models, namely the localized dynamic Smagorinsky model [19] and the localized dynamic one-equation SGS kinetic energy model [29], use a dynamic procedure as proposed by Germano [16], which automatically adapts the model coefficients in order to provide a proper near wall behavior. This paper is organized as follows. In Sect. 2, the numerical methods, the adaption technique for filtered statistics and the evaluation method are introduced. Next, the test case chosen to demonstrate the feasibility of the method, the fully developed turbulent pipe flow, is described in Sect. 3. Discrepancies issuing from reference data, aspects related to finite averaging time, finite domain size and finite initial transient are reported and quantified in Sects. 4-6. The evaluation of the SGS models in terms of accuracy and computational cost is discussed in Sect. 8. Finally, concluding remarks are provided in the last section.

Governing equations and modeling
In the case of a Newtonian constant density fluid flow, the applied governing equations of isothermal LES are the continuity equation and the momentum equation where P is the modified kinematic pressure, U i the resolved velocity field, ν the laminar viscosity and ν sgs the eddy viscosity of the residual content. Closure of Eq. (2) is thereby obtained by modeling the subgrid-scale viscosity ν sgs . In the present study, the employed SGS viscosity models can be generically represented as where C m is a model coefficient, Δ grid = Δ x Δ y Δ z 1/3 the grid filter width, and D m an operator which depends on the specific corresponding model. Six commonly used SGS viscosity models, namely the Smagorinsky model [58] with van Driest wall function [62], the one-equation SGS kinetic energy model by [69] with van Driest wall function, the WALE model [44], the SIGMA model [43], the localized dynamic Smagorinsky model [19] and the localized dynamic one-equation SGS kinetic energy model proposed by [29], are applied. In the present study, classical model coefficients are employed as proposed in the literature [43,44,58,69]. It follows for the constant in the Smagorinsky model that C S = 0.18, in the one-equation SGS kinetic energy model C k = 0.1, C e = 1, in the WALE model C W = 0.5 and in the SIGMA model C σ = 1.5. Notice, that in the case of wall-bounded turbulent flows, it is a common practice to use a reduced Smagorinsky constant of C S = 0.1 [11], which is not employed in the present study for the sake of consistency. In the case of localized dynamic models, the model coefficients are automatically computed using the least mean squares approach as proposed by [36]. Thereby, the least mean square formula is applied over the surrounding grid points of the actual cell, and remaining negative values are clipped to reduce the variability of the model coefficients over space.
In the case of SGS viscosity models using van Driest wall damping function, the grid filter width Δ grid is replaced by the expression where r + is the dimensionless wall distance, κ = 0.41, C Δ = 0.158 and A + = 26. The operators D m and the coefficient C m of the models applied in the present study are summarized in Table 1. Since the operators for the localized dynamic Smagorinsky and one-equation SGS kinetic energy models equal those of the corresponding models without dynamic procedure, they are simply omitted in Table 1.

Adaption technique for statistical moments in LES
In contrast to a DNS, the velocity U i and any statistic obtained by a LES are filtered quantities and thus depend on the grid resolution. To compare statistics obtained by a LES with reference data from DNS, the residual content of the statistics has to be determined and added to the filtered statistics [49]. Following the procedure described in [49,54], the mean value of the resolved velocity is approximately the temporal averaged velocity Table 1 Operators D m and model coefficients C m of the LES models

Model
Operator S i j denotes the resolved strain rate tensor, k sgs is the subgrid-scale kinetic energy, S d i j represents the traceless symmetric part of the squared velocity gradient and σ 1,2,3 are the three singular values of the velocity gradient predicted by a DNS. In contrast, the covariance matrix can be approximated by adding the temporal averaged residual stress tensor, thus where < . > t denotes temporal averaging. In Eq. 5, an expression for the subgrid-scale kinetic energy k sgs is required, which is usually unknown. Several models for estimating k sgs exist, for example [2,35,42,57,67].
In the case of wall-bounded flows, care must be taken such that the models are not too dissipative in the near wall region and, consequently, that k sgs vanishes as close to the wall. Therefore, this is accounted for in the present study by using the algebraic expression of [35] k sgs = ν sgs in which the eddy viscosity of the WALE model is utilized for all models under investigation. Thereby, the coefficient 1 is selected, consistent with the one-equation SGS kinetic energy [69] and with the other model coefficients in Table 1. Notice that it is also possible to use the eddy viscosity of the SIGMA model, since both models ensure that ν t , and consequently k sgs , vanish close to the wall.

Accuracy and computational cost
In order to evaluate numerical simulation results, the prediction accuracy has to be quantified [39,65]. For this purpose, the error measure of the simulation is estimated. It is defined as the deviation of the predicted statistic Q sim x and the corresponding statistic observed by the reference Q ref x interpolated at position x. In contrast to other evaluation studies of numerical simulations [22,24,40,64], the present error measure is normalized by the difference d(.) between the maximal and minimal value of the reference data Q ref , corresponding to the interval of interest. This leads to a normalized error as: which is a non-dimensional metric and therefore well suited to compare the prediction accuracy for different statistical quantities. Contrary to the relative error often used for evaluating simulation results, Eq. 7 never gives infinite values except in the irrelevant case where all reference data are equal. Alternatively, the error can be normalized using another characteristic quantity that has to be nonzero, e.g., mean or median value, which might be advantageous in cases of wide spreading reference data.
In the case of LES, discrepancies between the reference data and calculated flow properties arise from inaccuracies of the physical modeling ε mod x , the numerical error ε num x , the sampling error ε rand x , the influence of the initial conditions ε init x , the discrepancies in the boundary conditions ε bc x and the error of the reference data ε ref In order to determine the accuracy of a model, the error contributions ε num x , ε rand x , ε init x , ε bc x , ε ref x have to be much smaller than the inaccuracies of the physical modeling.
By means of explicit filtering, when the filter width and the grid spacing are independent, the modeling and the numerical error can be determined separately, such that a grid-independent solution can be reached [5]. In the context of implicit filtering as used in the present study, these two errors interact and are difficult to estimate separately as addressed in [6,30,64]. Therefore, we combine the modeling and the numerical errors to build ε LES x , which allows to characterize the prediction accuracy of the entire LES method. In order to obtain a global error metric of a predicted statistic Q sim , the normalized mean absolute error (nMAE) is introduced as where i is the location and n the number of spatial points. Following the criteria of appropriate global metrics [45], the nMAE is non-dimensional, nonnegative, symmetric and fulfills the triangle inequality. Furthermore, it can be used for datasets which have mean or median of zero (e.g., periodic data). The framework proposed here is intended to allow control of the trade-offs between accuracy, uncertainties and computational costs as function of the degree of fidelity expected. In this respect, it is important to address the required computational cost of a SGS model and also its scaling with regard to the spatial resolution. The relative computational cost of a subgrid-scale model C PU h * can be defined as the ratio of the CPU time spent for the calculation of the SGS model C PU h sgs and the total computation time of the simulation C PU h tot . It follows that Notice that the relative computational cost depends generally not only on the SGS model, but also on the selected test case, the particular code implementation and the solution procedure applied.

Test case of fully developed turbulent pipe flow
The methods and error metrics established in the previous section will be utilized to quantify the agreement between LES results and DNS data of a turbulent pipe flow at Re τ = 180 (based on the friction velocity u τ ). For this purpose, various LES simulations and a separate DNS (for comparison and evaluation) are conducted. The computational setup and a summary of the test cases are provided in this section.

Computational domain and numerical grids
A representation of the flow domain and numerical grid is shown in Fig. 1. Characteristic quantities of the numerical grids used in the current study are summarized in Table 2. Thereby, the numerical grid with the finest resolution (grid 5), which consists of more than 21 million control volumes, serves as a reference DNS.

Initial and boundary conditions
An isotropic turbulent velocity field is used to initialize the pipe flow simulations. Therefore, a random velocity field Ω is generated in a channel with equally spaced grid, which has a mean velocity value of zero by construction. In order to fulfill the continuity equation in Fourier space, the random field is cross-multiplied with the wavenumber vector κ and rescaled. In the next step, the autocorrelation spectrum similar to that found in [59,71], is imposed on each component of the random velocity field. Subsequently, the resulting random flow field U * is multiplied with a random phase n and inverse Fourier transformation F (κ) −1 is carried out leading to an initial fluctuation field where i is the imaginary unit. Afterward, u is linearly interpolated to the pipe geometry, and finally, a turbulent mean velocity profile [70] is superposed, where the velocity at the pipe wall is set to zero. Notice that the resulting initial velocity field generated by means of this procedure is only divergence-free away from solid bodies but not in the near-wall region. At the inflow and outflow, periodic boundary conditions are applied in stream-wise direction for the pressure and velocity. In the case of pipe wall, no-slip condition is set for the velocity and Neumann condition for the pressure. The pressure gradient, that drives the flow, is adjusted dynamically to maintain a constant mass flux. Therefore, the pressure is split into a periodic pressure and a nonperiodic pressure part [32]. A source term for the nonperiodic pressure gradient is added to the momentum equation. Regarding pimpleFoam, the velocity and pressure are coupled by means of a merged PISO [25] and SIMPLE [48] algorithm. Thereby, the governing equations are numerically solved using a momentum predictor, pressure solver and momentum corrector. In the present study, second-order central differencing scheme is applied for the convection terms and a second-order, conservative scheme for the Laplacian terms. In the case of convective transport of k sgs , a limited second-order scheme is used to avoid non-physical values. A second-order backward integration method is utilized for the time derivative terms. Thereby, the time step is chosen small enough to ensure that the CFL-number remains smaller than one. Convergence optimization and acceleration techniques are incorporated in order to speed up the calculations. In particular, geometric agglomerated algebraic multigrid solver is considered for the resolution of the pressure Poisson equation and for the momentum predictor. Convergence of the iterative procedure is obtained if all normalized residuals are smaller than 10 −4 . No model (DNS) 14 5 Grid no. see Table 2 (a) 3 10 50 200 Toonder et al [29] Ahn et al. [25] Eggels et al. [26] El Khoury et al. [59] Fukagata et al. [28] Wu et al. [27] Unger et al. [26] r + U +  [1,[13][14][15]68] and experiment [12] 3.4 Summary of the case studies Important features of the various numerical simulations are listed in Table 3. Cases 1-8 are used for a systematic variation study of the computational length. Thereby, the cell number in stream-wise direction is adjusted such that the grid width remains constant and equals grid no. 3 of Table 2.

Variations of the reference data (case 33)
Generally, in the available DNS, there is no complete information about the uncertainties with respect to the reported data. Furthermore, DNS results reported in the literature for the same configuration can differ from each other [63] and from experimental data. For the turbulent pipe flow at Re τ = 180 relevant to the present study, Fig. 2 shows a comparison of mean and rms velocity profiles of various DNS databases [1,[13][14][15]68] and of the experimental dataset provided in [12]. It should be mentioned that the profiles were extracted from graphs, and only the dataset of [15] was listed in an online database.
As it can be seen in Fig. 2a, mean velocity profiles of the DNS databases agree very well in comparison with each other and to the experimental data, while rms velocity profiles in Fig. 2b vary slightly, especially at the buffer layer and near the wall. Regarding the experiment, it is mentioned in [12], that reflections of the laser light from the pipe wall may cause these discrepancies. In the case of DNS databases, these variations are apparently smaller and may be caused by different numerical methods, spatial resolution, sampling errors and/or domain size, etc.
It is worth mentioning that an error for each DNS dataset is difficult to estimate since the true value is unknown, which impedes to assign the most appropriate data for an evaluation study. However, in the case several DNS of the same application are available, the variability of the datasets in comparison with each other can be quantified. For this purpose, Fig. 3 depicts the normalized standard deviation of mean and rms velocity datasets with respect to the non-dimensional wall distance. Thereby, the standard deviation is normalized by the minimal and maximal value of the respective statistic over the entire radius, corresponding to the region of interest.
In line with the visual observation in Fig. 2, it appears that the dispersion of predicted mean values is small (< 1%), while rms velocity data show slightly higher variations (≈ 1-4%). Furthermore, it can be seen in Fig. 3, that rms velocity data become more disperse toward the wall, with highest deviations at the buffer  4 Mean (a) and rms (b) velocities with respect to non-dimensional wall distance. Comparison of results from the present DNS with reference data of [1] layer (≈ 3-4%). Even though the overall variations are small in the case of DNS (≈ 2%), it can influence the outcome when evaluating LES results, especially in the case of highly resolved LES.
To avoid error contributions from the reference data caused by different numerical approaches, averaging time etc., it might be advantageous to use the same approach for generating the reference data as it is employed in the LES evaluation study. For the sake of consistency with the numerical approach applied for the present LES, an additional DNS is conducted for comparison purpose (see Table 3 case 33). To establish the validity of the present DNS, computed mean and rms velocity profiles are compared with the DNS dataset of [1] as shown in Fig. 4.
Apparently in Fig. 4, mean and rms velocities show excellent agreement with the DNS dataset of [1]. Using the normalized global error metric as introduced in Eq. 9, the discrepancies between the present DNS results and the reference data are ≈ 0.3% for the mean and ≈ 0.7% for the rms velocity components. Furthermore, considering the commonly used DNS spatial resolution criterion [23], the ratio of mean grid width Δ = (Δ r Δ ω Δ x ) 1/3 and the Kolmogorov length scale η K = ν 3 /ε 1/4 is below 1.5 in the entire domain, which ensures sufficient spatial resolution.

Sampling error (case 20)
In most experimental and numerical studies of turbulent processes, statistical moments are used to describe the properties of the flow field. These moments are usually estimated by a finite number of sample averages over a finite time period or over a finite spatial distance. They consequently contain stochastic errors. The central limit theorem can be applied to estimate these errors [4,34,61]. However, in contrast to experimental studies, this is rarely addressed in the case of LES or DNS.

Method
In the case of temporal averaging of stationary turbulent flows, statistically independent events occur after approximately two integral time scales, corresponding to the so-called temporal decorrelation distance. Fol-lowing the procedure described in [4,61], the number of statistically independent samples N t at position x can be estimated as N t = t av /(2τ i ), where t av is the averaging time and τ i the integral time scale of the velocity component i. Thereby, τ i can be calculated using the temporal autocorrelation function ρ i (s) of the time series (see, e.g., [52]) as where u i is the velocity fluctuation of U i , σ t U i the temporal standard deviation of the resolved velocity, t the time, s the lag-time, and < . > t denotes the temporal averaging. In this work, ρ i (s) is integrated only up to the first zero-crossing, representing the smallest time shift for which the velocity becomes decorrelated. However, other integration domains like integration over the entire available domain, up to the value where ρ i (s) is a minimum or up to the value where the autocorrelation function falls to 1/e, are also possible. For further discussion, please refer to, e.g., [47]. In addition, it should be noted here, that in the case of modest sample size, it is advantageous to estimate ρ i (s) by fitting an autoregressive time series models to reduce spurious oscillations as described in [46]. By means of the central limit theorem, the standard deviation of the estimator of the mean velocity U i,mean can be determined as σ t U i,mean = σ t U i / √ N t . Similarly, the standard deviation of the estimator of the rms velocity U i,rms is given for normal distribution as σ t U i,rms = σ t U i / √ 2N t . Normalizing σ t U i,mean and σ t U i,rms by the expectation of the resolved and the rms velocity at the same position x, respectively, leads to the following sampling errors of the mean and rms velocities Note that in the case of turbulent flows that are intrinsically not statistically stationary, e.g., in-cylinder flows in internal combustion engines or start-up processes, Eqs. (14) and (15) are not appropriate to estimate the sampling error. In such cases, ensemble-averaging has to be applied, whereby N t is the number of independent ensembles. Regarding statistically homogeneous flows, sample averages can be obtained by means of spatial averaging along a homogeneous direction. Thereby, a statistically independent event occurs after approximately two integral length scales, corresponding to the so-called spatial decorrelation distance, which can be estimated by means of two-point autocorrelation function R ik (see, e.g., [52]) as where r is the spatial lag-distance and e k is the unit vector in the x k -coordinate direction. Just as for ρ i (s), R ik (r ) is integrated only up to the first zero-crossing, representing the smallest length shift for which the velocity becomes decorrelated. In analogy to temporal averaging, the number of statistically independent samples N s along a homogeneous direction k can be estimated as N s = L/(2δ ik ), where L is the length of the homogeneous direction and δ ik the corresponding integral length scale. Using the central limit theorem leads to the following sampling errors for the spatial mean and rms velocities ε rand where < . > s denotes spatial averaging. Finally, regarding LES and DNS of turbulent pipe or channel flows, it is a common practice to reduce the stochastic error by means of spatial averaging of temporal averaged statistics, which is also done in the present study. In such cases, both temporal sampling errors ε rand t U mean and ε rand t U rms decrease by a factor of 1/ √ N s = √ 2δ ik /L.  Table 3)  . 6 Non-dimensional integral length scales δ + ik (a) and observational errors ε rand for a non-dimensional averaging time of t + av = 1.5e4 and additional spatial averaging in axial and azimuthal direction (b) as function of the wall distance. (case 20 in Table 3)

Results
Along the wall direction, Fig. 5 shows the non-dimensional integral time scales τ + i of the radial, azimuthal and axial velocity components predicted by the LES (case 20) and the corresponding observational errors ε rand t for a non-dimensional averaging time of t + av = 1.5e4, which is used in the present study. All time scales are non-dimensionalized in terms of viscous wall units (τ + i = τ i u 2 τ /ν). As it can be seen in Fig. 5a, integral time scales increase toward the wall and reach a maximum at the buffer layer (5 < r + < 30). Thereby, axial integral time scales are significantly larger than radial and azimuthal ones. In line with τ + i , the corresponding observational errors increase toward the wall and remain approximately constant at r + < 10 (see Fig. 5b). Highest values of ε rand t occur for the axial rms velocity (≈ 3.5%), while ε rand t is small in the case of mean velocity (< 2%). Obviously longer averaging time is required to obtain a desired temporal sampling error, when dealing with flow field statistics close to the wall.
In order to further reduce the sampling error, spatial averaging of the temporal mean and rms velocity components is applied along the axial and azimuthal directions. The integral length scales along these directions and the corresponding overall sampling errors after the spatial averaging are shown in Fig. 6a, b, respectively. All length scales are non-dimensionalized in terms of viscous wall units (δ + ik = δ ik u τ /ν). Note, that integral length scales in azimuthal direction for r + > 150 are omitted, since two-point autocorrelation functions do not reach zero in this region, leading to no additional independent samples.
Similar to the integral time scales, it is visible in Fig. 6a, that axial integral length scales (δ + r x , δ + ωx , δ + x x ) are small in the outer region (r + > 50), increase toward the wall and reach a maximum at the buffer layer. In contrast, integral length scales in azimuthal direction (δ + r ω , δ + ωω , δ + xω ) are large at the outer region and decrease toward the wall. Obviously, the number of independent samples depends inherently on the distance to the wall, the sampling direction and the velocity component when dealing with spatial averaging. It appears by comparing Figs. 5b and 6b that the stochastic errors can be reduced by a factor of approximately ten by means of spatial averaging. Especially radial and azimuthal velocity components in the viscous wall region (r + < 50) benefit from spatial averaging, due to the large number of independent samples in this region. It can be deduced from the present averaging time of t + av ≈ 1.5e4 and by means of additional spatial averaging in axial and azimuthal directions, that a sampling error less than 0.2% for mean and less than 0.4% for rms velocities is likely expected in the present study.

Influence of the domain size (cases 1-8)
Computed flow properties of LES and DNS are directly affected by the boundary conditions and/or the size of the computational domain, e.g., improper inflow boundary conditions [28], too short development lengths in the case of synthetic inflow turbulence [50] or uncertainties due to an insufficiently large computational domain [9]. In a large number of cases, an evaluation of how much such constrains are contributing to the uncertainty of the computed flow properties can be addressed by means of sensitivity tests or systematic variation studies. This is the focus of the present section.

Method
Regarding periodic boundary conditions as employed in the present study, care must be taken such that (1) the extent of the computational domain is sufficiently large to avoid any artificial constraints on the formation of fluid structures and (2) the computed turbulence statistics are no more affected by the size of the domain. For this purpose, a systemic variation study is applied in order to investigate the effects of the domain length on fluid structures and turbulence statistics. Therefore, the length of the domain L is gradually increased from L = 2D up to L = 16D with increments of 2D. Variations of the turbulence statistics are judged by means of the nMAE as introduced in Eq. 9, where Q ref corresponds to the statistics computed with the largest domain size (L = 16D). Furthermore, following the procedure described in [9], two-point correlations of the axial velocity in axial direction R x x are used to analyze the effects of the domain size on the formation of fluid structures (see Eq. 16). Thereby, a sufficiently large computational domain is needed that profiles of R x x converge to values close to zero. This is assessed by means of 95%-confidence levels estimated as C I 95 = 0 ± 1.96/ √ N , where N is the number of independent samples.

Results
To analyze the convergence of turbulence statistics, Fig. 7 shows the nMAE of the mean (a) and rms velocities (b) as a function of the pipe length.
As displayed in Fig. 7a, the nMAE of the mean velocity is fairly small (< 1%) and decreases rapidly with increasing domain size. Convergence is reached at a length of ∼ 6D, which corresponds to a non-dimensional distance of Δx + ≈ O (2000). Regarding the rms velocity components, the nMAE is more affected by the domain size associated with significantly higher values, especially for short pipe lengths. Convergence is reached at a domain size of approximately ten diameters, corresponding to a non-dimensional distance of Δx + ≈ O(3400) for all rms velocity components. Next, the spatial two-point correlations of the axial velocity component in stream-wise direction R x x are analyzed for different domain sizes at r + ≈ 15 (Fig. 8a) and at the centerline of the pipe (Fig. 8b). For the calculation, 4000 independent samples are utilized.
At the near wall region in Fig. 8a, it is apparent that correlations decrease with increasing pipe length. Significant correlations occur for L < 12D, where values of R x x do not reach zero, which indicates periodicity effects in this region. Convergence is reached at L ≥ 12D, corresponding to a non-dimensional distance of Δx + ≥ O(4000). Regarding spatial correlations at the centerline (Fig. 8b), profiles fall more rapidly and tend to zero at smaller lag distances compared to the near wall region. Convergence of the correlation curves is already reached at L ≥ 6D, which corresponds to a non-dimensional distance of Δx + ≈ O (2000).
In summary, to avoid any error contribution due to an insufficiently large pipe domain, a length of L ≥ 12D is required to achieve that mean and rms velocities and the formation of fluid structures are no more affected.

Influence of the initial transient (case 20)
To avoid uncertainties caused by the initial solution, the start-up phase of the simulation has to be long enough to ensure fully developed flow. Stationarity might be assumed, when first-and second-order moments of the signal do not vary with respect to time. In general, this can be tested for instance by visual inspection, statistical tests (see e.g., [3,53,55]) or by an analysis of statistical quantities for several segments of the recorded time series [27].

Method
In the present study, a segmentation technique is used to determine at which time the mean and the rms velocities are statistically stationary. For this purpose, the recorded velocity signal has to be much longer than the longest characteristic time scale of the process (t rec 2τ i ). First, estimations of the integral time scale, the mean and the rms velocities are calculated using the second half of the recorded signal. Thereby, it is assumed that the second half of the signal is stationary. Next, sample mean and rms velocities are calculated for consecutive segments of the dataset, each with a window size of t win = 30τ i . This window size is chosen to ensure a systematic error much smaller than the stochastic error. For the given window size, confidence intervals with a 95%-confidence level C I 95 are determined using the estimations of the second half of the recorded signal. This leads to for the mean and rms velocities [3], respectively. In order to calculate ε rand t , Eqs. (14) and (15) are applied with t av = t win . Stationarity of a signal is assumed when 95% of the sample mean and rms velocities averaged over the window size t win lie inside the corresponding estimated confidence intervals.  Figure 9 displays the start-up phase of the axial velocity signal (gray line) at the centerline of the pipe and its sample mean and rms values averaged over the moving window (black lines). Blue shaded areas represent the 95%-confidence region, and red lines highlight the time at which the statistics become statistically stationary.

Results
Apparently in Fig. 9, non-stationarity effects appear in the velocity signal during the start-up phase. Thereby, the sample mean velocity increases before it reaches stationary state at ≈ 80τ i . Afterward, the predicted values lie very well within the 95%-confidence limits. In contrast, in the case of sample rms velocity, stationarity is already reached at ≈ 30τ i . Therefore, it can be concluded that for the present initialization method, a cycle time of at least 80 integral time scales is required to ensure stationarity of first-and second-order moments. However, to avoid any uncertainties caused by the initial transient, a cycle time of 500 integral time scales is used in the present study, which ensures statistical stationarity of mean and rms velocities.

LES models evaluation (cases 9-33)
After minimizing the uncertainties of the simulations in Sects. 4-7, the evaluation of the physical behavior, the prediction accuracy and the computational cost of the selected LES models are addressed in this section. This is achieved by a systematic comparison of the obtained LES results with the generated DNS dataset using the error analysis (see Sect. 2) and by providing an estimation of the computational cost. Thereby, results of the Smagorinsky and one-equation models with van Driest damping using the two coarsest grids (cases 9, 10, 13, 14 in Table 3) are excluded from the analysis, since both models overdiffuse the flow so that it was not possible to get a fully developed flow field. Figure 10 shows the mean and rms velocity profiles predicted by the LES models in comparison with the generated DNS dataset. All velocity profiles are normalized by the friction velocity of the DNS. Results are exclusively shown for the Smagorinsky model with van Driest wall damping, the WALE model and the dynamic Smagorinsky model, as representatives for different near wall treatments. Similar results are obtained using the other models.

Near-wall physical consistency and comparison of statistical moments
As expected, results from the LES become more accurate with increasing spatial resolution and collapse onto the DNS data in the case of highest resolution. Regarding the axial mean velocity profiles (Fig. 10a), values are underestimated at the buffer layer region and near the wall. This tendency holds more or less for all models and spatial resolutions, but it is most notable for low mesh resolutions and for models using van Driest damping function. In the case of rms velocity profiles as depicted in Fig. 10b, results are more affected by the spatial resolution. Peak values of the rms velocities are shifted toward the log-law region, especially in the case of the Smagorinsky model with van Driest wall damping and low spatial resolution. This non-physical behavior is less pronounced for the wall-adapting and dynamic models. Regarding the rms velocity components without  Table 1). Dotted transparent lines denote the corresponding statistical quantities without adding the residual contribution (see Sect. 2.2) (color figure online) adding the residual contribution (dashed lines in Fig. 10b), values are inherently smaller, in particular for radial and azimuthal components. The benefit of adding the residual contribution will be discussed in detail later (see Sect. 8.2). Nevertheless, it can be observed here that results of radial and azimuthal components including the residual contribution are more faithful, while axial rms components are not affected by adding the residual contribution. It appears especially that models using a van Driest wall damping function underestimate the axial mean velocity at the buffer layer and feature a non-physical shift in the peak values of the rms velocity components. One reason for such inconsistencies might be an inappropriate near wall scaling of the eddy viscosity, which is analyzed in Fig. 11 for all LES models under consideration. Thereby, dashed lines represent the correct near wall scaling of the eddy viscosity (ν sgs = O(r 3 )).
It turns out that the wall-adapting models (WALE and SIGMA) well retrieve the theoretical asymptotic behavior near solid boundaries. Regarding the models using van Driest wall damping functions (Smagorinsky and one-equation model), the near wall scaling is well reproduced as well. However, a non-physical amount of eddy viscosity is produced at the buffer layer and near the wall. It seems likely that this apparent physical inconsistency is responsible for the shift in the rms velocity profiles, since the ratio ν sgs /ν does not differ significantly from the wall-adapting models in the remaining regions. Focusing on the localized dynamic Smagorinsky and one-equation models, both LES models are unable to reproduce the proper asymptotic behavior near the wall. As mentioned in [43], the proper asymptotic behavior can be only obtained in such models when the dynamic procedure is applied over homogeneous planes parallel to the walls, which is not feasible in complex geometries.
Summing up, it is found that first-and second-order statistics of the flow field are predicted more faithfully by the wall-adapting models and by models using the localized dynamic procedure than by models using van Driest wall damping function. Furthermore, only the wall-adapting models exhibit a proper near wall behavior of the eddy viscosity.

Error characteristics and prediction accuracy
The apparent performance of the LES models established by comparison with the generated DNS dataset in the previous section is now quantified by means of the error metrics as introduced in Sect. 2. For this purpose,  Table 1). Dashed lines represent the proper near wall scaling of the eddy viscosity (ν sgs = O(r 3 )) (color figure online) Fig. 12 Normalized error of predicted mean and rms velocities (grid no. 3, see Table 1). Smagorinsky with van Direst damping (blackline), one-equation with van Driest damping (redline), WALE (blueline), SIGMA (greenline), Smagorinsky with Germano (purpleline), one-equation with Germano (orangeline) (color figure online) Fig. 12 shows the normalized relative error of the mean and rms velocities with respect to the non-dimensional wall distance r + . Results are exemplarily shown for grid 3 (see Table 2).
As revealed in Fig. 12, deviations are small in the viscous sublayer (r + < 5), increase rapidly with increasing wall distance, reach a maximum at the buffer layer (5 < r + < 30) and finally decrease in the outer region (r + > 50). This trend in error characteristics is similar for all models and also for both statistics (mean and rms values). Furthermore, it can be seen that deviations are considerably higher in the case of models using a van Driest wall damping function, especially at the buffer layer and near the wall. In contrast, discrepancies in statistics predicted by wall-adapting and dynamic models are significantly smaller and comparable. Thereby, results of the dynamic models are slightly more accurate. After analyzing the error characteristics of the LES models for a given spatial resolution and for different flow regimes (viscous wall region, buffer layer, outer region), the overall prediction accuracy of the models with respect to the spatial resolution is examined next. For this purpose, the nMAEs of mean and rms velocities are calculated (see Eq. 9). Thereby, the locations, at which the nMAEs are computed, are logarithmically distributed along the pipe radius in order to obtain approximately an equal number of sample points in each flow regime (viscous sublayer, buffer layer, log-law region). Results with respect to the spatial averaged ratio of Kolmogorov length scale η and grid width Δ grid are depicted in Fig. 13. Thereby, η is computed for each location using the local turbulent kinetic energy dissipation rate obtained by the DNS. Solid lines represent the nMAEs of the rms velocities including the residual contribution and dashed lines the nMAEs of the rms velocities without adding the residual contribution.
As expected, the nMAEs decrease with increasing spatial resolution (η/Δ grid ↑). Especially for the walladapting and dynamic models, a lower spatial resolution is required to achieve an acceptable accuracy, reflecting a smaller modeling error of such models in comparison with LES models using van Driest wall damping. Differences in the prediction accuracy between wall-adapting and dynamic models are small. However, as shown in Fig. 13, the localized dynamic one-equation SGS kinetic energy model is the most accurate one for the selected test case. This holds for both mean and rms velocities over almost the entire resolution range. Furthermore, radial and azimuthal rms velocity components including the residual contribution have essentially lower nMAEs than those without using the adaption techniques described in Sect. 2.2. Axial rms velocity components do not benefit from this technique.
It turned out that wall-adapting and dynamic models feature a lower normalized relative error than models using a van Driest wall damping function, especially at the viscous wall region (r + < 50). In particular, a lower spatial resolution is required to achieve an acceptable accuracy of mean and rms velocities, especially for the localized dynamic one-equation SGS kinetic energy model, which is the most accurate LES model for the selected test case.

Estimation of computational cost
One of the key objectives for developing LES is to establish it as a model-based design method. This is possible if it is able to provide predictions using economically computational cost. Therefore, it is of practical interest to address the required computational cost of a LES model to achieve an acceptable accuracy. For this purpose, Fig. 14 shows the required relative computational cost for each subgrid-scale model as a function of normalized spatial resolution. Thereby, Δ grid represents the spatially averaged grid size and D the diameter of the pipe. For the estimation of the relative computational cost C PU h * (see Eq. 10), only one CPU-core was used and the maximal CFL-number was set to one.
As it might be expected, the Smagorinsky model with van Driest wall damping and the WALE model have the lowest relative computational cost of C PU h * ∼ 3%, whereas the one-equation model with Germano procedure and the SIGMA model are the most expensive ones with C PU h * ∼ 8%. The high computational cost of the one-equation model with Germano procedure is quite obvious due to the additional effort spent for solving the transport equation of k sgs along with the effort to automatically adapt the model coefficients. Less obvious is the high relative computational cost of the SIGMA model, which is purely algebraic and does not use any dynamic procedure. It appears that the calculation of the SIGMA model operator D m (see Table 1) is computationally intensive since the three singular values of the velocity gradient tensor have to be computed for each control volume and each time step. Nevertheless, the CPU time spent for the calculation of all subgrid-scale models under consideration is quite small compared to the total computation time of the simulation (< 10%). It can be therefore concluded that the computational cost plays a minor role in the case of such relatively simple subgrid-scale models.

Concluding remarks
The current work suggests a framework for comparative analysis of LES subgrid-scale models by implicating various sources of errors and uncertainties. This is demonstrated for a simple wall-bounded flow configuration, a fully developed turbulent pipe flow at Re τ = 180, where physical consistency and computational cost of six models, namely the Smagorinsky model with van Driest wall damping, the one-equation subgrid-scale kinetic energy model with van Driest wall damping, the wall-adapting local eddy-viscosity (WALE) model, the SIGMA model, the localized dynamic Smagorinsky model and the localized dynamic one-equation subgrid-scale kinetic energy model, have been additionally included. In particular, the sources of uncertainties in LES are estimated, systematically assessed and their effects on simulation results analyzed. Furthermore, various methods and measures to quantify and minimize these discrepancies are presented and discussed. Some highlights from the present study can be outlined as follows: 1. Examining the variability of available DNS datasets of turbulent pipe flow at Re τ = 180, it appeared that the dispersion in mean values is small (< 1%), while rms velocity data show higher variations (1 − 4%), in particular at the buffer layer. Such discrepancies can influence the outcome when evaluating LES results in comparison with DNS, especially in the case of highly resolved LES with prediction uncertainties of approximately 5%. Therefore, it might be advantageous to use the same approach for generating the reference data as it is employed in the LES evaluation in the present study (see Sect. 4). 2. When dealing with flow field statistics in the viscous wall region (r + < 50), it turned out that longer averaging time is required to obtain accurate results than away from the wall, due to large integral time scales in this region. Spatial averaging can help to reduce the sampling error in the case of statistically homogeneous flows (see Sect. 5).
Uc·ε 2 I is the turbulent intensity, τ the integral length scale, ε the desired maximal sampling error, l c a characteristic length scale of large turbulent scales and U c a characteristic velocity associated to the convection 4. Analyzing the influence of the initial transient, this study reveals that for an initialization method with synthetic isotropic turbulence, a cyclic time of at least 80 integral time scales is required to ensure stationarity of first-and second-order moments. This corresponds to an entrance length of L ≈ 20D (see Sect. 7). 5. In order to obtain accurate results in turbulent wall-bounded flows at moderate Reynolds numbers, an estimate of the required grid width in terms of non-dimensional grid spacing amounts + r ≤ 1, + ω ≤ 10, + x ≤ 20 (see Sect. 8). 6. Finally, it is found that wall-adapting subgrid-scale models (WALE and SIGMA) as well as localized dynamic models reproduce the physics of the flow field more faithfully while revealing a superior prediction accuracy than models using van Driest wall damping. Especially at the viscous wall region (r + < 50), wall-adapting and dynamic models are more accurate, reflecting the proper near wall behavior integrated into such models (see Sect. 8).
It is worth to note that some methods to quantify and minimize the uncertainties addressed in this work, especially estimations of the sampling error, can be very time-consuming and costly in particular for daily engineering and practical application with complex geometries. Engineering estimations for the practical use are then helpful. Therefore, a practical method to predict a priori the minimum record lengths to obtain basic estimates of statistical quantities, which is applicable for a large number of generic and engineering flow applications, is suggested in the appendix.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix: Practical approach for the required averaging time in LES
The temporal sampling error expressions derived in Sect. 5 provide ways to predict a priori the minimum record time to obtain basic estimates of statistical quantities with a predetermined degree of accuracy. By assuming stationary turbulent flow, normal distribution and an even autocorrelation function that follows the expression ρ i (s) ∼ e −τ/t av , the minimum total record lengths needed to achieve a desired accuracy are summarized in Table 4. Note that in the case of engineering estimations, it is further assumed that the hypothesis of "frozen turbulence" [60] is fulfilled and that ρ i (s) ≈ 0 for s ≥ 2τ .
In order to verify the engineering estimations of Table 4, estimates of the required averaging times at the centerline are compared with the present LES results. Therefore, l c = 0.25D, U c = 1.25 * U bulk , I = 0.04 and ε = 0.02 are chosen. This leads to a required averaging time for the mean velocity of t + av = 39 and for the rms velocity of t av = 1.21e4. These estimates are very close to t + av = 48 and t + av = 1.32e4, which are the required averaging times obtained directly by the velocity signal of the LES. In the case of complex wall-bounded flow geometries, in which the characteristic length scales and turbulent intensities are difficult to guess, RANS results can help to address this issue. In this context, the characteristic length equals the integral length scale l c = 0.43 · k 3/2 ε and the characteristic velocity equals the magnitude velocity. Thereby, k is the turbulence kinetic energy, and ε is the dissipation rate of k.