Scale-resolving simulations, such as large eddy simulations, have become affordable tools to investigate the flow in turbomachinery components. The resulting time-resolved flow field is typically analyzed using first- and second-order statistical moments. However, two sources of uncertainty are present when statistical moments from scale-resolving simulations are computed: the influence of initial transients and statistical errors due to the finite number of samples. In this paper, both are systematically analyzed for several quantities of engineering interest using time series from a long-time large eddy simulation of the low-pressure turbine cascade T106C. A set of statistical tools to either remove or quantify these sources of uncertainty is assessed. First, the Marginal Standard Error Rule is used to detect the end of the initial transient. The method is validated for integral and local quantities and guidelines on how to handle spatially varying initial transients are formulated. With the initial transient reliably removed, the statistical error is estimated based on standard error relations considering correlations in the time series. The resulting confidence intervals are carefully verified for quantities of engineering interest utilizing cumulative and simple moving averages. Furthermore, the influence of periodic content from large scale vortex shedding on the error estimation is studied. Based on the confidence intervals, the required averaging interval to reduce the statistical uncertainty to a specific level is indicated for each considered quantity.
Computational fluid dynamics (CFD) simulations are a well-integrated part in the design process of modern turbomachinery components. The current state-of-the-art technique relies on solving the steady Reynolds-averaged Navier–Stokes (RANS) equations with second-order accurate finite volume schemes in combination with appropriate RANS models, analyzing the steady-state flow field. Unsteady effects are assessed using unsteady RANS approaches or, to save computation time, frequency domain methods such as harmonic balance [1,2]. Yet, all the mentioned steady and unsteady approaches include RANS models, which are generally well tuned and mature at and close to the aerodynamic design points. However, compressors and turbines in off-design feature a high level of unsteadiness, flow separation, and laminar-to-turbulent transition, and these effects are typically difficult to reliably model, cf. Ref. . Methodologies of higher fidelity are needed in these cases to gain a deeper insight into the flow physics, to obtain high-quality reference data for the improvement of RANS models, or even to use them in the design-process itself. Naturally, scale-resolving simulations (SRS), such as direct numerical simulation (DNS) or large eddy simulation (LES), which spatially and temporally resolve all or most of the energetic content of the turbulent spectrum, are the methods of choice to close the knowledge gap for applications of engineering interest [4–6].
While scale-resolving simulations are still challenging in terms of computational requirements for multi-stage turbomachinery components, they are already commonly applied to simulate linear turbine or compressor cascades [7–10]. Due to the statistically stationary nature of the flow in linear cascades, main aspects of engineering interest in the analysis of SRS results are first- and second-order statistical moments. Quantities such as the averaged isentropic Mach number around the blade profile, pressure losses in a wake cut or Reynolds stresses and their budget terms are typically of interest and used to validate the results against experimental data or compare them with RANS simulations.
Unfortunately, two fundamental uncertainties are present when recording such statistical moments from SRS. First of all, the flow has to have progressed from arbitrary initial conditions to a fully-developed, statistically stationary state to allow to start the sampling of the statistical quantities. If the averaging process starts too early, transient phenomena can significantly influence the final statistical quantities. If the averaging process starts too late, high-quality samples are ignored and cost-intensive CPU-time will be wasted. Therefore, the reliable detection of the end of the initial transient is crucial. The other uncertainty relates to the averaging process of a chaotic process itself. Due to the limited computational resources and time, the statistical quantities can only be averaged over a finite number of samples and are, therefore, subject to statistical errors. These errors are seldom reported in the literature, but especially for higher-moments or low-frequency dominated quantities, they can be of a considerable magnitude depending on the averaging time. To add even more complexity, the initial transient time and the required number of samples for reliable statistical moments highly depend on local timescales and, hence, on the location in the flow field, e.g., wake or separation bubble, and on the sampled quantity itself. In order to fairly evaluate LES results and draw meaningful conclusions, these uncertainties would optimally be removed by taking an extremely long sample of data. Usually, this cannot be realized in practice. Hence, the uncertainties need to be at least quantified to be able to identify statistically significant differences. Looking at it from another direction, knowing and quantifying uncertainties of the averaged quantities of interest can also be used to reduce the simulation time by accepting a certain error level. Thus, reliable statistical error analysis can be one key aspect to make SRS of flows with very high Reynolds numbers feasible and to pave the way of SRS into the design process [11–13].
Indeed, the aforementioned uncertainties are not exclusive to scale-resolving simulations. Statistical moments obtained from experiments suffer from similar error sources, which are, in contrast to simulations often not as crucial. First, the recording of statistics in experiments is usually performed over durations much larger than the turbulent timescales. Thus, the measurement interval can be planned in advance to be sufficient to reduce the statistical errors without significantly increasing the cost if the dominant timescales are faster than minutes or seconds. On the contrary, the main cost factor in simulations is the simulation time itself, making it very attractive to have tools available which help decide on the appropriate length of the data acquisition. Second, other sources of uncertainties and systematic measurement errors such as boundary conditions, invasive measurement techniques, leakage flows, manufacturing tolerances, etc. are often dominant in experiments. These can easily conceal the small remaining statistical errors. In contrast to that, setups of different SRS can be aligned with only little effort. Hence, the comparability of different simulations is generally enhanced, reducing systematical discrepancies and statistical errors emerge.
To the author’s knowledge, limited work has been published in the LES community and especially for turbomachinery-related test cases focusing on the systematic detection of the initial transient and the quantification of the sampling error. In most of the publications, the reader is left with the information that the quantities have been averaged over a certain amount of time units after reaching a statistically steady-state. While this is positive for the reproducibility, the reader has to believe that the time estimates have been appropriate to remove or significantly reduce the statistical error and initialization bias. Besides being only a few, some studies have been published addressing these issues, but focusing mainly on academic test cases. For example, Oliver et al.  present a systematic and unified approach to estimate errors from two uncertainty sources, namely the finite statistical sampling and the discretization of the Navier-Stokes equations, tested on model problems and the DNS of a channel flow. The sampling error is estimated a posteriori from the correlated time series using classical standard error relations . Ries et al.  apply classical standard error relations to estimate and reduce the averaging errors and a windowed test for statistical stationarity to remove the initial transient for an LES of a pipe flow. Based on the outcome, the authors choose simulation parameters with minimized uncertainties to compare different LES models. Talnikar et al.  include the sampling error in an optimization to reduce the time-averaged drag from LES of a channel flow using traveling wave boundary conditions at the walls. Mockett et al.  utilize a window-based estimation of confidence intervals and the initial transient time. The presented methods are validated on a time signal of the integral lift coefficient obtained from a detached eddy simulation of the NACA0021.
In the current work, we focus on the detailed investigation of the initial transient time and the statistical convergence behavior of various quantities of engineering interest at different flow field locations of an LES. As a test vehicle, we use a popular test case for SRS code validation; the T106C cascade at an exit Reynolds number of 80,000 and exit Mach number of 0.65, which is an advanced test case in the High-Order CFD Methods workshop series. In Ref. , we investigated the impact of several properties of the LES setup on the mean quantities and second-order moments, i.e., the spanwise domain size, the grid resolution and the choice of boundary conditions.
Here, we concentrate on the uncertainties concerning the time-averaging, presenting and extensively analyzing a heuristic method, well-known in the field of statistical simulations, to automatically determine the end of the initial transient of LES data, namely the marginal standard error rule (MSER) . Moreover, we evaluate an established method to estimate the averaging error for correlated time series based on the standard error and an integral timescale, cf. Refs. [15,20], for the applicability to time signals obtained from the LES of the T106C. We put emphasis on the plausibility of the estimations and validate them through varying averaging window positions and sizes.
We remark that, apart from integral quantities such as blade forces or boundary values, all presented analyses are purely based on time signals probed at specific locations without exploiting spatial averaging. For statistically two-dimensional problems, such as the LES of the T106C cascade, an additional averaging in the spanwise direction could indeed reduce the statistical error if the domain size was sufficiently large such that the two-point correlation decays below zero. The latter however, cannot be guaranteed for either quantities or areas in the flow which are affected by the strong 2D vortex shedding . In the wake measurement plane, only the spanwise velocity component decorrelates completely within the computed domain, allowing the statistical error or number of required samples to be reduced by a factor of 2 at most. However, for configurations with end walls featuring pronounced three-dimensional effects, spatial averaging cannot be exploited at all. Therefore, we limit the presented analysis to the general case, knowing that for this particular test case class a further reduction of the statistical errors could be achieved using spatial averaging.
2.1 Identification of the Initial Transient.
Finding and removing the initial transient from the analyzed time signal is an essential part in the process of quantifying the statistical convergence of a simulation. Often, the duration of the initial transient is judged by intuitive inspection of time signals of representative quantities such as total pressure on the inflow/outflow planes or acting forces on the blade surface. This can however produce inconsistent results between publications, which was studied for example by Ref. . In this study, scientists were asked to visually identify the equilibrium time in molecular dynamics simulations and it was found that there is no mutual consent between the participants and, moreover, the appearance of the plot has a significant impact on the predictions. Several non-visual and automatable techniques can be found in the literature, cf. Ref.  for an overview. Following Ref. , the techniques can be sorted into five classes, namely graphical methods, heuristic approaches, statistical methods, initialization bias tests and hybrid methods. For the purpose of this paper, heuristic methods, which provide relatively simple rules to determine the truncation point, seem to be most suitable. We choose the MSER proposed by Ref. , which offers a good trade-off between computational efficiency and accuracy, cf. Refs. [23,24]. The method is designed to select a truncation point that minimizes the width of the confidence interval about the truncated sample mean.
2.2 Sampling Error.
Assuming we have successfully removed the initial transient of the time signal and the signal is statistically stationary, the sampling error could be easily estimated if the data points were independent and identically distributed (IID) using standard error relations. However, the single samples in signals obtained from scale-resolving simulations are generally not independent because they represent the temporal evolution determined by the Navier–Stokes equations. One possible approach is to sample the time signal with enough separation between the samples so that they can be treated as independent and apply standard relations for IID random variables, cf. Ref. . This approach can, however, lead to either an underestimation of the sampling error if the sample separation is too small. If the sample separation is too large, the sampling error will be overestimated as valid samples are ignored.
Note that different approaches to estimate the autocorrelation are possible. For example, Oliver et al.  fit an autoregressive model and directly use Eq. (5) to compute the decorrelated sample size. This can be beneficial for modest samples sizes, i.e., DN ∼ N. However, the number of samples obtained from SRS of turbine cascades is usually several orders of magnitude larger than the estimated number of correlated samples and, thus, no significant difference in the outcome has been found in the preliminary studies between fitting an autoregressive model and using Eq. (7).
For the sake of completeness, in modern statistics and analysis of measurement data, bootstrapping is a prominent alternative to the classical standard error relations to estimate statistical error and confidence intervals. Especially, since it does not rely on the central limit theorem, it allows for simple estimations of confidence intervals of complex quantities. Benedict and Gould  compare, for example, the estimated uncertainties from bootstrapping and from the classical approach for higher-order velocity fluctuations of turbulence measurements made from a flow over a backward-facing step using a laser Doppler anemometer. They found that the results of both error estimation approaches agree well. Thus, we stick to the classical approach in the remainder of this paper, but recognize that bootstrapping is an interesting alternative for the desired use case, which could be investigated in future work.
3 LES of the T106C Low-Pressure Turbine Cascade
The T106C low-pressure turbine (LPT) cascade has been used as a test vehicle for the validation and comparison of numerical simulations from RANS to DNS in various forms of abstraction. The experiment at the VKI was arranged as a linear cascade of six prismatic blades with an aspect ratio of 2.4. The blades with a chord length of c = 93 mm are staggered under an angle of , which leads to an axial chord length of cax = 79.9 mm. This highly loaded case has a pitch to chord ratio of 0.95. At a fixed isentropic exit Mach number of 0.65, a test matrix of Reynolds numbers ranging from 80,000 to 250,000 with free-stream turbulence intensities of 0.9%, 1.8%, and 3.2% was published by Michalek et al. . Here, we consider the case with Reynolds number of 80,000 and a laminar inflow, which were also the boundary conditions of the case at the High-Order CFD Methods workshops. The selected operating conditions produce a laminar separation in the aft section of the suction side. In the separated shear layer, two-dimensional Kelvin-Helmholtz rollers develop and eventually break down into three-dimensional turbulent flow. In the averaged flow field, reattachment can be seen very close to the trailing edge (TE). The wake is dominated by a strong vortex shedding such that two-dimensional structures can be observed far downstream of the blade.
This laminar inflow case can be compared with the measurement data with a turbulence intensity of . The agreement of results from scale-resolving simulations with the experiment, however, is rather poor as was shown in several studies [7,33,34]. Deviations are present in terms of the blade loading near the leading edge (LE), a delayed separation bubble and a clear underestimation of the wake depth. There have been attempts to mitigate the differences in blade loading by adapting the inflow angle, but discrepancies in the appearance of the separation bubble and the wake depth remained . We will see that, by using the statistical error analysis performed in Sec. 4.2, the uncertainty due to the averaging process can be excluded as a reason for these deviations. Hence, it is relatively safe to assume that systematic issues in the reproduction of the experiment are the dominant factor, e.g., the assumption of spanwise periodicity is not justified and the growing endwall boundary layers affect the measurements at midspan due the relatively small aspect ratio of the cascade [7,18,33]. Despite the known deficiencies, the test case has been heavily used as a validation case in the High-Order CFD Methods workshops for different CFD solvers and the observed differences between numerical results are far smaller [10,18,34], indicating the relevance of thorough statistical error analysis to compare numerical approaches. The experimental results are, nevertheless, shown in this paper to put the numerical results into perspective.
DLR’s flow solver for turbomachinery applications, TRACE, was used to perform an LES of the above configuration. TRACE’s finite volume method solves the filtered compressible Navier–Stokes equations using a second-order accurate, density-based scheme applying MUSCL reconstruction with κ = 1/3 . A fraction of 10−3 of Roe’s numerical flux  is added to a central flux to avoid odd-even decoupling. Time integration is performed using a third-order accurate explicit Runge-Kutta method. The subgrid stresses are computed by the wall-adapting local eddy-viscosity (WALE) model . A one-dimensional non-reflecting boundary condition, which drives the time and surface averaged boundary state towards prescribed values (stagnation pressure, stagnation temperature and flow angles at inflow, static pressure at outflow) by means of incoming characteristics, is used for both inflow and outflow .
The computational domain is sketched alongside with turbulent structures visualized by spanwise vorticity in Fig. 1. The inflow plane is located 0.1 m = 1.075c upstream of the blade leading edge which was shown to be sufficient for the non-reflecting boundary conditions used here . To be able to analyze a longer portion of the wake, the outflow plane is 0.126 m = 1.351c downstream of the blade TE. A slice of one blade is computed with periodic boundary conditions in the pitchwise and spanwise directions. We choose a spanwise extent of 0.3c in contrast to 0.1c in the workshop case after carefully analyzing the setup of the computation, which showed a strong dependency of the dynamics in the wake on the spanwise extent for values lower than 0.3c. Even with a spanwise extent of 0.4c, the two-point correlations of the in-plane velocity components do not go to zero in the measurement plane at x = xTE + 0.4cax .
The domain was meshed with DLR’s in-house block structured mesh generation tool PyMesh . Typically, turbine blades are meshed with an OCGH-topology consisting of an O-type block around the blade surface, wrapped by a C-type block to ensure a smoother transition to the H-type blocks downstream of the blade. H-type blocks are also used upstream and in the blade passage. The G-type block connects the passage, wake and exit H-blocks introducing a singularity at which only 3 grid lines come together in 1 point. For the LES, additional H-type blocks were inserted downstream of the passage and in the near wake to ensure a high-quality mesh aligned with the mean flow direction for a longer section than normally possible in multi-stage configurations. On the suction side, the number of points is roughly 2.5 times larger than on the pressure side. Table 1 summarizes the properties of the mesh.
4 Analysis of Time Series
In this section, we provide an in-depth analysis of the time-dependent behavior of the test case shedding light on the initial transient time and the statistical convergence. Moreover, the applicability of the methodology to determine the end of the initial transient and to evaluate the sampling error for several quantities at different locations will be reviewed. Time signals at various positions of several quantities of interest have been probed over 355 convective times with a sampling frequency of . To set this in relation, Tyacke and Tucker  proposed an averaging time of 2 to 10tc as a best practice guideline for cascade LES. Due to storage constraints, it is prohibitive to sample well resolved time series of this length and resolution for the complete 3D domain. Hence, subsets of the domain are written to HDF5 files in order to be later processed using standard python libraries. Figure 1 shows the positions in the flow field, which are probed for the complete simulation duration. All probes shown here are located at a single spanwise position, i.e., at the midspan.
4.1 Time of the Initial Transient.
When planning the data acquisition for an LES, one of the first questions that arises is, at what point in time the averaging procedure should be started. If the averaging starts too early, transient effects can significantly influence the mean results and, therefore, bias the conclusions made from the LES. If the averaging starts too late, valuable samples are not considered in the process and, thus, reaching a statistically converged state requires more computational time. However, we will show in the following section, that the initial transient time varies highly for different quantities and different positions in the flow field. So, if one wants to or has to average multiple quantities at multiple locations with the same interval, one has to accept that at least one of the two options is going to occur.
The actual duration of the transient phase depends highly on the quality of the initialization. In this study, we initialize the simulation with a converged steady RANS solution. Because no resolved fluctuations are present in the RANS solution, the transient phase consists of 2D vortex shedding which eventually breaks down into 3D turbulence, and consequently a potential shift in mean quantities. If one uses a better initialization, e.g., the instantaneous flow field of LES on a coarser mesh, the initial transient time will probably decrease, but nonetheless be present.
4.1.1 Integral Forces on the Blade.
A typical starting point to quantify the time of the initial transient found in the literature is to analyze the temporal evolution of the integrated blade forces . The time signal of the normalized force in pitchwise direction Fy(t) up to is given in Fig. 2(a) and shows a significant initial range, which differs qualitatively from the remaining time signal. After a certain point in time, the changes of the force have an almost periodic character, oscillating around a mean value. Using MSER to quantify this point, results in an estimated end of the initial transient of , which seems reasonable by visual judgement. To further analyze the initial transient, Fig. 2(b) shows the truncated mean of the integrated force, which is the average of the signal starting from the end of the simulation following Eq. (2). One can observe that differences of the truncated mean to the full mean are very small even if only few convective times are considered. After the estimated end of the initial transient marked by the vertical line, we note a significant increase of the truncated mean despite having already averaged over more than 340 convective times. The decrease of the truncated mean right before the vertical line is within the order of magnitude of the other variations and, thus, not considered as a part of the initial transient. Overall, MSER yields a reliable estimate for this integrated quantity.
4.1.2 Boundary Conditions.
The non-reflecting boundary conditions at the inflow and outflow are weakly enforced. Therefore, the instantaneous surface averaged boundary values can differ from the boundary values prescribed by the user. Naturally, these signals also feature an initial transient before reaching a statistically stationary state as the boundary condition controller reacts to the changing inner flow field. The initial transient times computed with Eq. (1) for different quantities at the inflow and outflow boundaries are presented in Table 2. Compared to the blade forces, longer initial transient times can be observed for quantities at the inflow as well as at the outflow. The defining feature, which needs to develop from the 2D RANS initialization, is the separation bubble. Until , the axial velocity in the separation region is strongly correlated in the spanwise direction because large 2D structures separate from the blade, which decay into three-dimensional vortices afterwards. After , the laminar roll-ups begin to become turbulent directly inside the separation bubble, i.e., separation-induced transition starts to develop. Thus, the initial transient of the integral blade force ends shortly after the characteristics of the separation are in a reasonably steady-state. Then, the information needs to progress from the aft region on the suction side of the blade to the boundary planes and, as a consequence, the observed initial transient times at the in- and outflow are longer. The initial transient time at the inflow plane is affected by acoustic waves, which move upstream from the developing separation bubble and are reflected between the blades. The relevant characteristic at the outflow is the turbulent wake. At first glance, it is counterintuitive that the initial transient of the turbulent phenomena ends earlier than the laminar one, but we will make and discuss similar observations when considering signals inside the domain in the following section. Generally speaking, due to the low variance present in signals of laminar regions, the MSER is more sensitive to events which are triggered by the initial transient progress. On the other hand, developing large fluctuations in turbulent regions tend to cover events from the transient phase leading to shorter transient time predictions.
4.1.3 Field Values.
In addition to the analysis of integral values, it is instructive to investigate how the estimated end of the initial transient varies in the flow field. As time series data are required, only a limited subset of the domain can be used for this analysis. Figure 3(a) shows the analysis for the static pressure in the recorded suction side boundary layer probes, a relatively coarsely sampled region around the TE of the blade and two cuts through the wake. An instantaneous vorticity field is shown in grey scale for orientation. Several observations can be made: The quantity shows a mostly smooth, yet significant variation throughout the domain from very short transients up to more than . The largest initial transient times are detected in the laminar region near the LE. Interestingly, these values are significantly larger compared to the estimation from the static pressure of the inflow boundary, cf. Table 2. The reason for this could be that these are point probes, which are probed at the midspan, whereas the static pressure of the inflow boundary is an integral quantity. Thus, spanwise and pitchwise variations are averaged out and the fluctuations over time can be smaller. Looking at the raw time signal A in Fig. 3(b), a visible change in the amplitude can be identified after the estimated end of the initial transient. The level of fluctuations is significantly higher before due to upstream moving acoustic waves from the still not fully developed separation bubble. These acoustic disturbances, which occur in the initialization phase, are especially visible in signal A due to the relatively low variance at this position. The MSER processes each time signal individually and therefore disturbances, which are relevant for the laminar region but irrelevant when compared to the variances of turbulent regions, are treated equally. In combination with a desired error bound of statistical quantities, the MSER could be improved by ignoring variations below a certain limit.
In contrast to the location near the LE, short initial transient times ranging from 0 to are estimated in the region of the separation and around the TE. One can see that the signals B and C indeed exhibit only small portions of transient phenomena. After less than , the fluctuations in signal B, probed in the separation region, are already in a statistical steady-state. The point C, which is located near the TE, features a nearly constant signal at the beginning of the simulation. After around , the variations are directly in the same magnitude as for the remaining signal. Therefore, the truncated MSE does not increase and no initial transient is detected. While this is initially unintuitive, it makes sense when thinking in terms of averaging the signal. It means that if the values at the beginning of the signal are already very similar to the long-time mean and the fluctuations are below the fluctuations of the remaining signal, the statistical error of the mean is lower when the beginning of the signal is included. Thus, the long-time mean is more reliable in a statistical sense by averaging the full signal. Consequently, instead of determining the overall initial transient of the signal, the MSER rather detects the initial transient relevant for the average of the signal. To make the estimation more general, one could also include errors of higher statistical moments, such as the variance, to account for changes in these terms. However, for the scope of analyzing time-averaged results of LES, the MSER definition is sufficient. Another interesting point can be identified in the wake probes, marked by D. Its estimated initial transient differs significantly from the estimations of neighbor points, i.e. to . The time plot of signal D indeed shows an increase in the truncated MSE before . The increase is, however, only very small and, thus, the global minimum of the truncated MSE is at including more samples. If the truncated MSE exhibits a very flat minimum as in this case, small differences in the time signals can result in large changes in .
Summarizing the above observations, no clear trend can be identified as the estimated initial transient depends highly on the given initial field. Laminar regions tend to have a longer initial transient due to the fact that variations in the beginning of the simulation are more significant compared to the final steady-state fluctuations. However, if no variations are present, the estimated initial transient is also low as observed in the region on the suction side in front of the separation, cf. Fig. 3(a). Overall, the detected initial transient times from integral values cannot be simply applied to the flow field. Consequently, the initial transient has to be detected separately for each probe and, concluding from the investigated time signals, the MSER is appropriate for that purpose, although potential improvements have been addressed.
When the simulation is started from a given initial solution, we have observed in the previous sections that the initial transient depends on the quantity of interest and the location in the domain. Thus, finding a generally appropriate beginning of the averaging for multiple quantities and locations can be quite difficult. Furthermore, we showed that only considering integral quantities, such as the blade forces and boundary values, can lead to a false estimation of the initial transient of the inner flow field. Based on that, we formulate the following guidelines with regard to the starting point of the averaging for the remaining paper:
A time signal is available for the quantity and location: Use the initial transient of the time signal.
A time signal is not available for the quantity and location: Use the maximum transient time of all globally available time series as the beginning of the averaging.
The averaged quantity of interest is a function of multiple instantaneous quantities for which time series are available: Average each quantity individually with its own initial transient. Thus, the end of the transient of the derived quantity is the maximum of the initial transient of the involved signals.
Time series are available for the quantity on a subdomain, which should be averaged, e.g., blade cuts: Use the maximal detected initial transient time in the subdomain. One could also think of using each individual initial transient in this case. However, we found that exploiting programming benefits by using homogeneous numpy nd-arrays is more appropriate.
4.2 Statistical Convergence and Error.
4.2.1 Isentropic Mach Number.
We observe a good agreement of the LES results with the measurements on the pressure side. In contrast, even for the largest averaging time, significant discrepancies are present in the leading edge region on the suction side and in the region of the separation between x/cax = [0.58, 1]. These discrepancies are neither within the confidence interval of the LES nor within the uncertainties specified in Ref. , which range from to ( is shown in Fig. 4, but is only visible in the zoomed view due to the marker size). Hence, the differences cannot be linked to statistical errors due to insufficient averaging time. As discussed in Sec. 3, similar deviations have been observed in other publications.
Focusing on the LES results, we cannot distinguish between the results with different averaging lengths in the full view Fig. 4(a). Zoomed closer into the separation region at the suction side TE, see Fig. 4(b), small differences can be spotted between the short averaging times and the long-time average. Most of the long-time average lies well within the estimated confidence interval of the short-time averages, which is at maximum for and for averaging window. Yet, in the aft region between x/cax = 0.97 and x/cax = 1.00, the long-time average is biased towards the upper border of the CI of the average over . This could be an indication of a slow transient drift, which is more pronounced in the total pressure losses inside the wake and will be discussed in detail in Sec. 4.2.3. Nevertheless, the deviations are small and the relative CI is already below , such that the isentropic Mach number can be regarded as well converged after only .
To further address the statistical convergence and survey the estimated CI, the time signals of two locations with a more significant estimated error, marked with vertical lines in Fig. 4(b), are separately analyzed in Fig. 5. The first point is located in the separation region at x/cax = 0.93 and the second near the TE at x/cax = 0.997. In the upper subplot, the time signal of an instantaneous isentropic Mach number using the instantaneous reference total pressure p0,inlet(t) and static pressure p(t) is shown. We note a significant initial transient in the time signal at both locations, whose duration corresponds with the estimated transient duration of the inflow boundary condition. The subplot in the middle shows the CMA of the isentropic Mach number (blue line) following the guidelines from Sec. 4.1.4. Additionally, the CI of the mean is plotted for each averaging interval. To visually judge the evolution of the mean, the long-time-averaged isentropic Mach number is also shown in the figure as a constant over the normalized simulation time (orange line). Differences between the long-time mean and the short-time means, i.e., averaging window sizes of up to , are clearly visible in Fig. 5(a). However, the long-time average lies well within the predicted CI of the short-time averages. This is further substantiated by the lower subplot of Fig. 5(a), where estimated error is plotted against the observed difference between the CMA and the long-time average. The observed error is always below the estimated error, which indicates the plausibility of the estimate. We emphasize here that one realization of an averaging window has to be interpreted as one experiment and, following the definition of a confidence interval, there is a chance that the confidence interval encompasses the true mean. Hence, it is not guaranteed. As the mean values with longer averaging intervals also contain the same samples of the short-time mean, they are to some extend linked together and follow a similar bias. Looking at Fig. 5(b), deviations of the long-time mean and the CI of the CMA between to can be noted. Though, after until the shown , the CI encompasses the long-time average as already observed for the probe in the separation. Considering the estimated and observed errors, one can see that the observed errors match the estimated ones very well when averaging the signal until to . On the other hand, there are intervals in which the observed error is comparatively small. If only a small series of selected averaging intervals is chosen to judge statistical convergence without considering the estimated averaging error, misleading conclusions on convergence may be drawn.
Summarizing the observations, the isentropic Mach number around the blade converges statistically in a fairly short time. The maximal predicted CI is already below after averaging the flow field for only . Analyses of time signals from two locations indicate that the used estimation is plausible even if it rather overestimates the observed error for this signal.
4.2.2 Skin Friction Coefficient.
The MSER estimates the maximum initial transient of the distribution to , which is longer than the involved signals from the boundary conditions and, thus, determined by one of the shear stress signals. All averaging lengths result in similar distributions of the mean skin friction coefficient featuring a long separation bubble and a secondary recirculation region in which the backflow separates from the blade. The separation point of first bubble is predicted to be at in the long-time reference, which was determined through the zero crossing of the mean shear stress wall in the axial direction. The shorter averaging intervals differ by a maximum of (). Similar observations can be made for the location of the secondary recirculation region differing by a maximum of with respect to the long-time average. Hence, the key features are already well converged after averaging for . Uncertainties are only present for the magnitude of the friction coefficient in the secondary bubble and near the TE on the suction side.
A way to study the uncertainty and the plausibility of the CI is to keep the averaging duration constant and move the averaging window along the time signal, i.e., using SMAs. This is carried out using non-overlapping averaging intervals of in Fig. 6(b) for the position marked in Fig. 6(a). First, a significant variation of the mean values depending only on the position of the averaging window can be noticed. The measured mean skin friction coefficient, visualized by the dashed red lines, ranges from to , which is a difference of relative to the long-time mean. The high uncertainty is also reflected in the width of the CIs, which is on average of the long-time mean. Comparing both, 30 of the 33 CIs from independent averaging windows contain the reference mean and only 3 realizations do not. This seems to be low, as in ideal stochastic conditions, 10 to 11 of the 33 CIs should not encompass the true mean. Though, to approach these conditions, vastly more realizations and, therefore, longer simulation times are required, which is not realizable at the moment. Consequently, keeping the low number of realizations and possible correlations in mind, the plausibility of the estimated CI is at least indicated.
Concluding from an engineering perspective, the large relative uncertainty of the friction coefficient in that specific location is probably not important enough to justify longer averaging times as key aspects are already well converged. Hence, with the knowledge and estimation of the uncertainty, an averaging window of or even less is enough for the friction coefficient in this case.
4.2.3 Total Pressure Loss 0.465cax Behind the Trailing Edge.
Interestingly, the estimated confidence interval does not vanish in the region between the wakes where almost no turbulent structures are present and the mean pressure loss is equal for all averaging times. This can be related to the phenomenon where the periodic content of the signal results in an offset of the estimated error, which will be discussed in Sec. 4.3. Overall, in comparison to quantities on the blade, e.g., the isentropic Mach number and shear stress, the predicted relative uncertainties are about an order of magnitude larger and visually noticeable variations occur between the different averaging windows, especially in the middle of the wake.
In Fig. 7(b), the observed drift in the wake profile is investigated using the CMA (top) and two SMAs (bottom) of the total pressure loss coefficient at the position marked in Fig. 7(a). While the CMA of isentropic Mach number shown in Fig. 5 converged relatively quickly toward the full average of the signal, a different behavior can be observed here. It takes roughly for the confidence intervals of the cumulative and the full averages to touch. As can be seen in the SMA with a window of , there is a transient behavior on a very long timescale which could not be detected by MSER. With this knowledge, the full average should be taken over the interval 100 → 355 leading to a loss, which differs by . This is just outside of CI of the average over 8.1 → 355. The slow drift is extremely hard if not impossible to detect without having access to the full time sample. This is illustrated by the green lines in Fig. 7(b) which show the analysis which could have been performed if only of time signal had been available. The CMA (top) still shows a rather noisy signal within this range but the value moves within the confidence interval of the average over 8.1 → 100 from as early as about . The fluctuations in the SMA (bottom) appear to be random and lie well within the CI, which is, of course, considerably bigger due to the reduced window size. All things considered, it is very likely that one could have called the short simulation statistically converged with a relative estimated error on the total pressure loss coefficient at this position of , when, in fact, this result deviates more than from the average over 100 → 355. In terms of integral loss coefficient, the deviation is not as big as detailed above but still because the effect can be observed over large parts of the pressure loss profile.
4.2.4 Shear Stress 0.465cax Behind the Trailing Edge.
Reynolds stresses extracted from LES are often used in RANS modeling. It is essential to be certain of their statistical convergence to draw the right conclusions. The error analysis for the shear stress component is shown in Fig. 8. Because the Reynolds stresses are second moments, they are expected to converge slower than the mean values. This translates into larger statistical errors and confidence intervals, which are present in Fig. 8(a). We can observe a variation of the shear stress distribution with increasing averaging time, while the confidence intervals intersect each other. At the position of peak stress marked by the dashed vertical line, the one-sided confidence interval after averaging over is around and after still around .
Again, the plausibility of the estimator is validated in Fig. 8(b) for the marked location. The estimated confidence interval covers the observed differences for most averaging times. Note that the long-time average is only an estimate for the true mean and by definition of the CMA the observed error must eventually go to zero. A slow transient phase similar to that in the pressure loss can be observed as well. However, in this case, the variation of the SMA occurs mostly within the CI, so the drift is statistically less significant. Note that especially from on, the CI appears to be overestimated. It has to be considered, though, that the SMA with a window size of 100 yields just about 3 statistically independent samples for the mean. To reliably assess the quality of the estimator for such a broad window size, the simulation would have to run at least 10 times longer. Averaging with a window size of gives a hint that very low frequency oscillations are present in the shear stress. This motivates longer averaging times such as which seems to be more appropriate in this case judging by the rather stochastic oscillations of the corresponding SMA.
4.2.5 Boundary Layer Analysis.
The incoming laminar boundary layer () has no visual statistical error and is, not surprisingly, equally well predicted for all averaging times. After , the appearance of the separation is nearly identical for all averaging times with only little uncertainties. The largest uncertainties are present at the boundary of the developing separation bubble, i.e. cut , and in the boundary layer cut near the TE, where strong turbulent phenomena are present. Compared to the quantities in the wake, the mean tangential velocity profiles along boundary layer cuts do not require such a long averaging time to decrease the uncertainties to reasonable values, i.e., maximal or standard error relative to the mean outlet velocity in case of averaging until or , respectively. If one is only interested in the overall structure of the separation bubble, even averaging windows below are sufficient.
To get an indication for laminar to turbulent transition, we investigate the development of the maximum resolved turbulent kinetic energy in the boundary layer over the axial coordinate , see Fig. 9(b). The shaded area represents the CI of the respective location of the maximum k, which is computed following Sec. 2.2. A clear, continuous, but slow increase of the maximal turbulent kinetic energy starting from can be observed. Interestingly, between and , the level of the resolved turbulent kinetic energy is significantly higher compared to the other shown intervals at the blade cuts . The CI barely intersects the other CIs, which shows the high uncertainty in that region. Despite the uncertain magnitude, the location of the transition onset is very similar for all averaging windows at with maximum difference of . It has been determined by the crossing of two tangential lines approximating the two slopes. The first one is the best fit of the region between and , and the second one is the best fit between and .
The largest uncertainties of the maximum resolved turbulent kinetic energy are present in the fully turbulent region after . The maximum differs from , when averaging from to , to using an averaging window from to . Nevertheless, the CIs intersect each other and the long-time reference lies well within the CIs. Hence, longer simulation times might be required if the exact value is of interest. On the contrary, the overall appearance of the maximum turbulent kinetic energy converges faster and averaging over is already suitable to draw conclusions about the basic transition phenomena.
4.2.6 Integral Total Pressure Loss.
We choose a short section of of the signal starting at the estimated end of the initial transient of and use it to estimate the number of correlated samples as well as the variance. For a relative error on of , Eq. (12) yields a required averaging time of . Indeed, the statistical error obtained from the simulation over (excluding the initial transient) is estimated very consistently as . In light of the measurement uncertainty , a specified relative error of yields a required averaging time of , showing that a first, rough evaluation of integral quantities is already possible after averaging only over few convective timescales. This demonstrates another use of the error estimation approach assisting in the planning of scale-resolving simulations.
4.3 Error Estimation in the Presence of Periodic Signals.
In the previous section, we showed the plausibility of the confidence interval based on the sample variance and integral timescale motivated decorrelation times for different quantities of interest. The applied error estimation procedure is usually suited for random processes including correlated samples. However, especially in the wake, significant low-frequency periodic phenomena are present due to huge size of the laminar separation bubble and its roll-up, which can be seen in Fig. 1 as wavy structure in the vorticity contours. When we look at the power spectral density (PSD) of the axial velocity sampled at the location marked in Fig. 7(a) shown in Fig. 10, dominant frequencies exist at f = 1.7/tc, 2.3/tc, and 3.0/tc. Especially in laminar regions such as between two wakes, the estimated error can be dominated by periodic oscillations.
As discussed before, the estimator for the MSE with a given confidence interval of can be interpreted in the following way: Suppose a series of independent random experiments for which a sample mean is computed. Then this sample mean will lie within the confidence interval of the true mean in of the cases. The statement can be reformulated such that the actual standard deviation of the sample means of all experiments is equal to the average standard deviation estimated from the single experiments. This is easily confirmed using normally distributed numbers from any random number generator for sufficiently large samples of the order of 100 samples per experiment and 100 experiments. Using a lower number of samples per experiment or experiments will yield more spread in the number of means within the given confidence interval.
The results are summarized in Table 3. We can observe how the standard deviation of the sample means grows only very slowly with increasing amplitude of the synthetic signal. The estimated error, however, grows much faster, leading to an increasing ratio of sample means within the estimated confidence interval of the true mean. While the error is already slightly overestimated for the base signal (which can be due to the low number of samples), the overestimation increases with increasing periodic content. This has to be considered when assessing flows in which stochastic turbulence is superimposed with periodic content such as vortex shedding or passing wakes from previous blade rows.
|0.00||0||0.76||9.91 · 103||7.45 · 103||1.33|
|0.02||5||0.80||10.32 · 103||7.45 · 103||1.38|
|0.04||5||0.86||11.50 · 103||7.52 · 103||1.53|
|0.08||5||0.94||15.66 · 103||7.83 · 103||2.00|
|0.00||0||0.76||9.91 · 103||7.45 · 103||1.33|
|0.02||5||0.80||10.32 · 103||7.45 · 103||1.38|
|0.04||5||0.86||11.50 · 103||7.52 · 103||1.53|
|0.08||5||0.94||15.66 · 103||7.83 · 103||2.00|
The above analysis is also applied to time signals from the wake of the T106C at x = 0.465 cax. To this end, the pitchwise probes are separated into two regions, i.e., inside the wake defined by probes with and outside the wake containing the remaining probes. An example spectrum for each region is shown in Fig. 10, displaying the significant difference in their ratio between periodic and turbulent content. For the considered signals, the true mean of quantities of engineering interest such as the average total pressure is generally not known. We can, however, use the long-time average as an estimator for the true mean. With a window size of 6tc, 42 windows are obtained for this interval with the effective sample sizes ranging from 32 outside the wake to 81 in its core. Outside the wake, the ratio of the estimated standard deviation to the computed standard deviation of the sample means σestimated/σsample means increases to 3.8, showing an overestimation of the error as expected due to the limited strength or total absence of turbulent fluctuations. Inside the wake, on the other hand, an average ratio of 1.08, including portions with ratios below 1, indicates an accurate estimation of the error. This further substantiates the findings of the previous sections, in which the estimated error at several locations of interest has been analyzed with CMAs and SMAs and found to be plausible. Finally, because the estimated error inside the wakes is still larger by a factor of up to 3.7, the overprediction of the error between the wakes would not lead to false overall conclusions regarding the required averaging time.
5 Summary and Conclusion
A set of known tools to systematically assess the statistical convergence of SRS results has been introduced covering two main ingredients: the detection of the end of the initial transient and an estimator for the confidence interval of sample means based on correlated time signals. We have applied the methods to a long sample of data obtained from an LES of the T106C LPT at an exit Reynolds number of 80,000 and exit Mach number of 0.65. In general, the marginal standard error rule (MSER) method has been shown to reliably detect the initial transient except for very slow drifts. We formulated guidelines on how to deal with spatially varying initial transients to obtain the largest of statistically independent samples possible. With the initial transient reliably removed from the time signals, we assessed the statistical convergence of various quantities of engineering and turbulence modeling interest. The short averaging times in the order of 10 convective times found in the literature can be justified for the blade pressure distribution and skin friction.
A different conclusion needs to be drawn for the wake flow. Significantly longer averaging times in the order of 100 convective times are required to obtain confidence intervals of below for the total pressure loss. Not surprisingly, higher moments such as the Reynolds stress tensor show even larger relative confidence intervals. We observed that signals with a significant periodic content compared to the turbulent fluctuations need to be treated with caution because there, the error can be consistently overestimated.
An interesting observation was made in terms of an extremely slow drift seen especially in the wake. This slow drift could not be observed in the boundary values and integral quantities. We discussed that even with the information of confidence intervals this would have been extremely hard to spot with only short data samples available. This might be an artifact of this concrete case and initialization, but the possibility of effects like this has to be considered when LES results are evaluated.
Considering possible slow drifts and periodic content, further work is required to fully automate estimates of the mean and the confidence interval. This can involve the specification of a tolerable error when determining the end of the initial transient as well as a reliable estimation of error in the presence of periodic signals. With robust error estimates incorporated in the online analysis of CFD solvers, statistical convergence criteria could be formulated such that the simulation terminates when the specified error is reached. However, even with the known limitations, the discussed methods can be used to base the duration of the computation on a tolerable confidence interval in terms of statistics. As LES enters engineering applications more and more, this can potentially save a lot of resources and avoid misleading conclusions.
Conflict of Interest
There are no conflicts of interest.
- CFD =
computational fluid dynamics
- CI =
- CMA =
cumulative moving average
- DNS =
direct numerical simulation
- IID =
independent and identically distributed
- LE =
- LES =
large eddy simulation
- LPT =
- MSE =
- MSER =
marginal standard error rule
- PSD =
power spectral density
- RANS =
- SMA =
simple moving average
- SRS =
- TE =
- WALE =
wall-adapting local eddy-viscosity
- c =
- f =
- g =
example signal from an LES
- k =
resolved turbulent kinetic energy
- p =
- t =
- E =
expected value of a random variable
- N =
number of samples
- T =
- cf =
skin friction coefficient
- eN =
estimated error based on N samples
- tc =
characteristic convective time ()
- ttransient =
estimated end time of the initial transient
tangential velocity component
- xtrans =
axial coordinate of the transition onset
- DN =
number of samples between two independent observations
- Mis =
isentropic Mach number
vector of the integral forces on the blade surface
- VarN =
vector of the Cartesian velocity components
tangential vector on the blade surface pointing in the positive mean flow direction
vector of the Cartesian coordinates