Abstract

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.

1 Introduction

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. [3]. 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 [46].

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 [710]. 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 [1113].

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. [14] 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 [15]. Ries et al. [16] 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. [13] 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. [17] 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. [18], 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) [19]. 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 [18]. 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 Methodology

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. [21]. 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. [22] for an overview. Following Ref. [22], 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. [19], 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.

Given a finite number of samples N of the time series obtained from an LES {gi : i = 0, …, N − 1}, the final sample of the initial transient $dN,M*(g)$ is determined by solving the following minimization problem:
$dN,M*(g)=argmin0
(1)
where $g¯N,d$ is the truncated mean defined as
$g¯N,d=1N−d∑i=dN−1gi$
(2)
Here, M is the number of realizations of the function to be minimized, which has been set to the total number of samples N in the original formulation. In the current work, we assume that the initial transient is only present in a small portion of the signal and set M = N/2 to speed-up the computation [24]. In the remainder of the paper, the term to be minimized in Eq. (1) is referred to as the truncated mean-square error (MSE). We remark that newer versions of the MSER, namely the MSER-5 [23,25] and the MSER-5Y [26], have been developed in the meantime. Instead of the raw time signal, the MSER-5 uses non-overlapping batch means with size 5 as the input signal in Eq. (1). This leads to a reduction of the sample size and, therefore, improves the computational speed of the algorithm. As the authors stated, the batch size of 5 is “somewhat arbitrary” [24], but showed the best results for the considered problems. However, with increasing batch sizes, the estimation of the truncated MSE becomes less accurate. Especially for signals which have a flat region around the minimum of the truncated MSE, these inaccuracies affect the determined initial transient times. Regarding the signals analyzed in Sec. 4.1.3, $97%$ of the initial transient times determined by MSER and MSER-5 deviate by less than $0.1tc$. Extreme deviations are found up to $5tc$. In the scope of the current paper, we aim for the most accurate solution and use the original MSER formulation applied to the raw time signal with a carefully chosen sampling frequency, which is also used to compute the sampling error. The improvements of the MSER-5Y focus on problems that occur if the initial transient is detected in the second half of the signal, i.e. $dN,M*>N/2$. This was not the case in the current study.

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. [27]. 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.

In this paper, we follow a more sophisticated approach based on the direct estimation of the sampling error from correlated data, cf. Refs. [15,20]. Given the signal obtained from LES {gi : i = 0, …, N − 1}, the sampling error of the population mean $μ=E[g]$ is defined as
$eN(g¯)=g¯−μ$
(3)
where $g¯=1N∑i=0N−1gi$ is the sample mean. Due to the finite simulation time, the population mean μ is generally not known. However, for IID random variables, the MSE of the sample mean $g¯$ can be estimated as
$MSE(g¯)=E[eN2(g¯)]=VarN(g)N$
(4)
where $VarN(g)=1N∑i=0N−1(gi−g¯)2$ is the biased estimate of the population variance [15]. In order to account for the correlation of the samples from SRS, a minimal number of samples between two independent observations can be defined [20,28] as
$DN(g)=1+2∑τ=1N−1(1−τN)ρN,τ(g)$
(5)
where
$ρN,τ(g)=1N−τ∑i=τN−1(gi−g¯)(gi−τ−g¯)VarN(g)$
(6)
is a straight-forward estimate of the autocorrelation with lag τ. Using this estimation, the autocorrelation typically obtained from SRS of blade profile flows tends to be noisy and exhibits oscillations. Thus, we follow the thoughts in Refs. [20,29] and integrate the autocorrelation only until the first zero crossing to estimate the number of samples between independent observations DN. Moreover, the sample size of our data is typically significantly larger than the lag of the first zero crossing. Thus, we assume that τ/N is negligible until the first zero crossing and Eq. (5) finally reduces to
$DN(g)≈1+2∑τ=1untilρN,τ<0ρN,τ(g)$
(7)
This leads to the estimation of the MSE for correlated time series used in the remainder of the paper
$eN2(g¯)≈DN(g)VarN(g)N$
(8)
where Neff = N/DN can be interpreted as the effective number of samples.

Note that different approaches to estimate the autocorrelation are possible. For example, Oliver et al. [14] 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., DNN. 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).

In the following analysis, we will study the statistical errors of quantities of engineering interest, always clearly stating the respective averaging procedure. To keep the overhead in the CFD solver as low as possible, only time signals of basic quantities are probed, i.e. primitive variables and their spatial gradients. Time signals of derived quantities such as the total pressure can be computed in the post processing step. The errors of the quantities involving averages from multiple sources such as the isentropic Mach number on the blade are computed using error propagation neglecting the correlations between the source quantities. We employ the python package gvar to automatically perform the error propagation for derived quantities [30]. The error of second-order moments such as Reynolds stresses is estimated using the variance of the combined signal and the maximum of the number of decorrelated samples obtained from the respective single signals, see Eq. (7), yielding
$eN2(g1′g2′¯)≈max{DN(g1′),DN(g2′)}VarN(g1′g2′)N$
(9)
where the total number of samples N is equal for both signals. The error of derived quantities, e.g., the turbulent kinetic energy, are again propagated using gvar. In the further analysis, the estimated error is used to define a confidence interval (CI) for the statistical quantity based on the outcome of the central limit theorem. In the remaining paper, we use the 68 % and the 95 % CI assuming a normal distribution defined as:
$CI68%≈[g¯−1eN;g¯+1eN]$
(10)
$CI95%≈[g¯−1.96eN;g¯+1.96eN]$
(11)

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 [31] 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.

The method described above can not only be used to estimate the error for statistical moments of a given signal but can also be employed to estimate the required number of samples N to obtain a specified error by rearranging Eq. (8) to
$N=DN1(g)VarN1(g)eN2(g¯)$
(12)
Note that in this case, $DN1(g)$ and $VarN1(g)$ are estimated using only a smaller number of samples N1N.

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 $30.7∘$, 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. [32]. 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 $0.9%$. 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 [7]. 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 [35]. A fraction of 10−3 of Roe’s numerical flux [36] 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 [37]. 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 [38].

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 [18]. 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 [18].

Fig. 1
Fig. 1
Close modal

The domain was meshed with DLR’s in-house block structured mesh generation tool PyMesh [39]. 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.

Table 1

Overview of mesh and time-step properties

 Number of cells spanwise 90 Total number of cells/106 14.82 $Δx+¯$ 12.05 $maxΔxSS+$ 17.12 $y+¯$ 0.21 $maxySS+$ 0.56 $Δz+¯$ 11.51 $maxΔzSS+$ 30.40 Number of time-steps per tc 28,569
 Number of cells spanwise 90 Total number of cells/106 14.82 $Δx+¯$ 12.05 $maxΔxSS+$ 17.12 $y+¯$ 0.21 $maxySS+$ 0.56 $Δz+¯$ 11.51 $maxΔzSS+$ 30.40 Number of time-steps per tc 28,569

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 $tc$ with a sampling frequency of $fs=460.8/tc$. To set this in relation, Tyacke and Tucker [11] 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 $F(t)$. The time signal of the normalized force in pitchwise direction Fy(t) up to $100tc$ 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 $ttransient=5.7$, 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.

Fig. 2
Fig. 2
Close modal

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 $4tc$, 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 $4tc$, 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.

Table 2

End of initial transient ttransient/tc for various quantities at the inflow and outflow boundary

Quantityttransient/tc
Inflow boundary
Stagnation pressurep0,inlet8.1
Stagnation temperatureT0,inlet8.1
Flow angleαinlet8.1
Static pressurepinlet8.2
Outflow boundary
Static pressurepoutlet7.0
Flow angleαoutlet7.8
Mass flow$m˙outlet$6.0
Quantityttransient/tc
Inflow boundary
Stagnation pressurep0,inlet8.1
Stagnation temperatureT0,inlet8.1
Flow angleαinlet8.1
Static pressurepinlet8.2
Outflow boundary
Static pressurepoutlet7.0
Flow angleαoutlet7.8
Mass flow$m˙outlet$6.0

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 $ttransient/tc$ shows a mostly smooth, yet significant variation throughout the domain from very short transients up to more than $15tc$. 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 $16.5tc$ 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.

Fig. 3
Fig. 3
Close modal

In contrast to the location near the LE, short initial transient times ranging from 0 to $5tc$ 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 $3tc$, 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 $1tc$, 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. $5tc$ to $7.5tc$. The time plot of signal D indeed shows an increase in the truncated MSE before $7.5tc$. The increase is, however, only very small and, thus, the global minimum of the truncated MSE is at $5tc$ 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 $ttransient/tc$.

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.

4.1.4 Guidelines.

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:

1. A time signal is available for the quantity and location: Use the initial transient of the time signal.

2. 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.

3. 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.

4. 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.

In this section, we address the statistical convergence of the LES of the T106C by analyzing time signals of relevant quantities at different locations in the flow field. We start with the blade loading and the pressure distribution on the blade represented by the isentropic Mach number. Following that, uncertainties in the position of the separation are investigated through the wall shear stress. The statistical convergence of the transition behavior and general appearance of the separation bubble is judged using boundary layer probes. Finally, we study the total pressure loss and turbulent shear stress in the wake. For all these quantities, we analyze the validity of the relations of the statistical error given in Sec. 2.2 using a simple moving average (SMA) and a cumulative moving average (CMA). The different averaging windows starting from tstart and ending at tend are denoted by the superscript tstarttend, i.e.,
$g¯tstart→tend=1tend−tstart∫tstarttendg(t)dt$
(13)
Unless stated otherwise, the shaded areas shown in the following figures represent the $68%$ CI of the respective mean quantity.

4.2.1 Isentropic Mach Number.

The first quantity of interest is the mean isentropic Mach number, which is defined as
$Mis¯(x)=2γ−1((p0,inlet(t)¯p(x,t)¯)γ−1γ−1)$
(14)
where p0,inlet(t) is the area averaged stagnation pressure taken from the inlet and p(x, t) is the static pressure on the blade surface. In Fig. 4, the isentropic Mach number averaged for 10, 20 and $337.9tc$ is shown and compared with the experimental results of Ref. [32]. Following the guidelines given in Sec. 4.1.4, we start the averaging after an initial transient of $17.1tc$, which is restricted by the transient time of the probes near the LE, cf. Fig. 3(a). Additionally, the $68%$ CI based on the estimated error is plotted as a shaded area for each averaging duration. As stated in Sec. 2.2, the estimated error of the Mach number is propagated from the individual averages of p0,inlet(t) and p(x, t).
Fig. 4
Fig. 4
Close modal

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. [32], which range from $0.4%$ to $1.2%$ ($0.8%$ 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 $68%$ confidence interval of the short-time averages, which is at maximum $0.8%$ for $10tc$ and $0.64%$ for $20tc$ 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 $20tc$. 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 $68%$ CI is already below $1%$, such that the isentropic Mach number can be regarded as well converged after only $10tc$.

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 $68%$ 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 $12tc$, are clearly visible in Fig. 5(a). However, the long-time average lies well within the predicted $68%$ 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 $68%$ 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 $68%$ CI of the CMA between $10tc$ to $15tc$ can be noted. Though, after $15tc$ until the shown $100tc$, 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 $75tc$ to $93tc$. 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.

Fig. 5
Fig. 5
Close modal

Summarizing the observations, the isentropic Mach number around the blade converges statistically in a fairly short time. The maximal predicted CI is already below $1%$ after averaging the flow field for only $10tc$. 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.

To investigate the separation behavior of the T106C, we are interested in the mean skin friction coefficient around the blade, which is given by
$cf¯(x)=sgn(t(x)⋅τ(x,t)¯)||τ(x,t)¯||p0,outlet(t)¯−poutlet(t)¯$
(15)
where p0,outlet(t) and poutlet(t) are the area averaged total and static pressures at the outlet and $τ(x,t)$ is the wall shear stress vector. The sign of cf is derived using the scalar product of the wall shear stress vector and the tangential vector of blade surface t(x) defined in the span plane, i.e., tx > 0 and tz = 0. Again, a separate error is estimated for each signal and error propagation is applied to determine the error of the derived quantity. The mean skin friction coefficient distribution including the $68%$ CI is presented in Fig. 6(a).
Fig. 6
Fig. 6
Close modal

The MSER estimates the maximum initial transient of the distribution to $16.6tc$, 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 $x/cax=0.676$ 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 $±0.002cax$ ($0.3%$). Similar observations can be made for the location of the secondary recirculation region differing by a maximum of $±0.003cax$ with respect to the long-time average. Hence, the key features are already well converged after averaging for $10tc$. 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 $10tc$ 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 $cf¯min≈0.13⋅10−3$ to $cf¯max≈1.02⋅10−3$, which is a difference of $158%$ relative to the long-time mean. The high uncertainty is also reflected in the width of the $68%$ CIs, which is on average $94%$ 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 $10tc$ or even less is enough for the friction coefficient in this case.

4.2.3 Total Pressure Loss 0.465cax Behind the Trailing Edge.

The error analysis of the wake velocity deficit at 0.465cax downstream of the TE is plotted in Fig. 7 in terms of the total pressure loss coefficient including the $68%$ CI over the relative pitchwise coordinate s. The mean total pressure loss is given by
$ζ¯(s)=1−p0(s,t)¯p0,inlet(t)¯$
(16)
where p0(s, t) is the total pressure at the pitchwise position s in the wake and p0,inlet(t) is the total pressure area averaged over the inflow panel. Significant differences between the pressure losses determined by the LES and the experiment of Ref. [32] are present. The peak pressure loss from the experiment nearly doubles the value of the simulation. On the contrary, the width of the wake matches fairly well in all cases. As discussed in Sec. 3, the assumptions made in the numerical simulation, such as the use of laminar inflow conditions rather than a turbulence intensity of $0.9%$, as reported in the experiment, can have an effect on the size of the separation and the mixing in the wake. The reported measurement uncertainty of 10–$20%$ ($15%$ is shown as shaded band in the figure) certainly cannot explain the deviation alone. The observed differences to the experiments are also larger than the estimated confidence intervals for all plotted averaging times. The total pressure loss in the experiment was measured with 5-hole probes and, therefore, the instantaneous total pressure has been averaged directly to obtain the shown total pressure loss. However, this procedure is not the only possible one to evaluate the total pressure loss from scale-resolving simulations. Another choice could be to average the primitive variables directly, which is in the context of a CFD solver more natural and could seem arbitrary at a first glance. Yet, due to the presence of turbulent fluctuations inside the wake, e.g., quantifiable by a maximal normalized turbulent Mach number of $VarN(M)/M¯=0.106$, the averaging strategies result in deviating predictions of the total pressure loss as shown in Fig. 7. This fact has to be considered and clearly stated when showing and comparing averaged quantities. Seemingly, the wrong averaging strategy results in a closer agreement with the experimental data, which is most likely due to a cancellation of errors.
Fig. 7
Fig. 7
Close modal

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 $170tc$ for the confidence intervals of the cumulative and the full averages to touch. As can be seen in the SMA with a window of $100tc$, 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 $3%$. This is just outside of $95%$ 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 $100tc$ 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 $30tc$. The fluctuations in the SMA (bottom) appear to be random and lie well within the $95%$ 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 $2.9%$, when, in fact, this result deviates more than $11%$ from the average over 100 → 355. In terms of integral loss coefficient, the deviation is not as big as detailed above but still $6.3%$ 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 $u′v′¯$ 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 $≈92tc$ is around $7%$ and after $347tc$ still around $4%$.

Fig. 8
Fig. 8
Close modal

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 $95%$ CI, so the drift is statistically less significant. Note that especially from $200tc$ 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 $30tc$ gives a hint that very low frequency oscillations are present in the shear stress. This motivates longer averaging times such as $100tc$ 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.

To evaluate the separation and transition behavior of the LES, boundary layer cuts normal to the blade surface have been placed along the suction side in the midspan plane (cf. Fig. 1). The tangential velocity at these boundary layer cuts gives insight into the appearance of the separation bubble, which is visualized for different averaging times in Fig. 9(a). The average tangential velocity is computed via
$uζ¯(η)=u(η,t)¯tx+v(η,t)¯ty$
(17)
where (tx, ty) is the tangential vector pointing along the blade surface in the positive mean flow direction and η is the distance from the blade surface.
Fig. 9
Fig. 9
Close modal

The incoming laminar boundary layer ($0.65cax$) has no visual statistical error and is, not surprisingly, equally well predicted for all averaging times. After $0.8cax$, 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 $0.8cax$, 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 $5.0%$ or $3.4%$ standard error relative to the mean outlet velocity $||uoutlet¯||$ in case of averaging until $30tc$ or $50tc$, respectively. If one is only interested in the overall structure of the separation bubble, even averaging windows below $30tc$ 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 $k(η)=12(u′(η,t)2¯+v′(η,t)2¯+w′(η,t)2¯)$ over the axial coordinate $x/cax$, see Fig. 9(b). The shaded area represents the $68%$ 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 $0.65cax$ can be observed. Interestingly, between $100tc$ and $200tc$, the level of the resolved turbulent kinetic energy is significantly higher compared to the other shown intervals at the blade cuts $0.7. The $68%$ 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 $0.89cax$ with maximum difference of $±0.002cax$. 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 $0.7cax$ and $0.85cax$, and the second one is the best fit between $0.89cax$ and $0.95cax$.

The largest uncertainties of the maximum resolved turbulent kinetic energy are present in the fully turbulent region after $0.96cax$. The maximum differs from $k/||uoutlet¯||2=0.127$, when averaging from $18.7tc$ to $100tc$, to $k/||uoutlet¯||2=0.137$ using an averaging window from $100tc$ to $200tc$. Nevertheless, the $68%$ 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 $100tc$ is already suitable to draw conclusions about the basic transition phenomena.

4.2.6 Integral Total Pressure Loss.

From an engineering perspective, it is highly relevant to be able to estimate simulation runtimes which are required to drive the statistical error for a blade efficiency criterion below a certain threshold. One such criterion is the integral total pressure loss, which can be evaluated from area averaged total pressure at inlet and outlet as
$ζ¯outlet=1−p0,outlet(t)p0,inlet(t)¯$
(18)
Note that this differs slightly from the definition given in Eq. (16) in terms of averaging order but enables the direct application of Eq. (12). The former definition can also be used in combination with error propagation and leads to similar results.

We choose a short section of $2tc$ of the signal starting at the estimated end of the initial transient of $7.76tc$ and use it to estimate the number of correlated samples as well as the variance. For a relative error on $ζ¯$ of $1%$, Eq. (12) yields a required averaging time of $297tc$. Indeed, the statistical error obtained from the simulation over $297tc$ (excluding the initial transient) is estimated very consistently as $eN(ζ¯)≈1.09%$. In light of the measurement uncertainty [32], a specified relative error of $10%$ yields a required averaging time of $3tc$, 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.

Fig. 10
Fig. 10
Close modal
To investigate the effect of periodic content superimposed on turbulent signals, we perform a test based on fully developed turbulent channel flow at $Reτ=395$ computed with an eighth order Discontinuous Galerkin scheme with Kennedy–Gruber split formulation on a mesh of 64 × 64 × 64 as described in Ref. [40]. This dataset should contain no significant periodic content, apart from that induced by the periodic domain of 2πδ × 2δ × πδ with the half channel height δ. The simulation ran for 40 eddy turnover times from which 20 were used to obtain converged averages of mean velocity and Reynolds stresses. We use a time signal of the spanwise velocity w from a probe in the logarithmic part of the boundary layer at y+ ≈ 151 whose mean has to be exactly 0 due to the symmetry of the configuration. On top of that we impose a synthetic periodic signal with a varying amplitude a and a single frequency f
$w*(t)=w(t)+asin(2πft)$
(19)
The PSD of the original in Fig. 10 (a = 0, f = 0) shows a clear inertial range with low frequency peaks not significantly greater than the noise. With the proposed synthetic signal, a discrete peak at a frequency slightly below the inertial range is generated, which stands out of the noise. The concrete frequency is arbitrary in this case, but the same trends can be observed nearly independent of the frequency.

As discussed before, the estimator for the MSE with a given confidence interval of $68%$ 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 $68%$ 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.

For the time signal described above, following Eq. (13) we define a series of non-overlapping windowed means as
$g¯m:=g¯tm→tm+Δtm$
(20)
for a window m with size Δtm with tm+1tm + Δtm as uncorrelated experiments with a known true mean. Note that the number of independent samples per experiment of roughly 30 and number of experiments of roughly 50 are still relatively low, so we cannot expect a perfect agreement with the statement in the previous paragraph. What we can indeed observe, is an indication of what happens with increasing amplitude of periodic signal. For M uncorrelated random experiments, we obtain a standard deviation of the sample means
$σsamplemeans=VarM(g¯m)$
(21)
This should be roughly equal to the standard deviation of the mean estimated from Nm samples averaged over all windows m
$σestimated=1M∑m=0M−1eNm(g¯m)$
(22)
by Eq. (8). The number of samples which lie within the estimated confidence interval of the true mean is given by
$Minside=|{g¯m|μ−σestimated
(23)

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.

Table 3

Quality of the error estimator for synthetic signals with increasing amplitude superimposed on the spanwise velocity of turbulent channel flow at y+ ≈ 151

$aU0$$fhU0$$MinsideM$σestimatedσsample means$σestimatedσsamplemeans$
0.0000.769.91 · 1037.45 · 1031.33
0.0250.8010.32 · 1037.45 · 1031.38
0.0450.8611.50 · 1037.52 · 1031.53
0.0850.9415.66 · 1037.83 · 1032.00
$aU0$$fhU0$$MinsideM$σestimatedσsample means$σestimatedσsamplemeans$
0.0000.769.91 · 1037.45 · 1031.33
0.0250.8010.32 · 1037.45 · 1031.38
0.0450.8611.50 · 1037.52 · 1031.53
0.0850.9415.66 · 1037.83 · 1032.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 $ζ¯>0.01$ 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 $pt¯100→355$ 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 $5%$ 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.

Nomenclature

Abbreviations

• CFD =

computational fluid dynamics

•
• CI =

confidence interval

•
• CMA =

cumulative moving average

•
• DNS =

direct numerical simulation

•
• IID =

independent and identically distributed

•
• LE =

•
• LES =

large eddy simulation

•
• LPT =

low-pressure turbine

•
• MSE =

mean-square error

•
• MSER =

marginal standard error rule

•
• PSD =

power spectral density

•
•     RANS =

Reynolds-averaged Navier–Stokes

•
• SMA =

simple moving average

•
• SRS =

scale-resolving simulations

•
• TE =

trailing edge

•
• WALE =

Latin Symbols

• c =

chord length

•
• f =

frequency

•
• g =

example signal from an LES

•
• k =

resolved turbulent kinetic energy

•
• p =

pressure

•
• t =

time

•
• E =

expected value of a random variable

•
• N =

number of samples

•
• T =

temperature

•
• $m˙$

mass flow

•
• cf =

skin friction coefficient

•
• eN =

estimated error based on N samples

•
• tc =

characteristic convective time ($c/||uoutlet¯||$)

•
• ttransient =

estimated end time of the initial transient

•
• $uζ$ =

tangential velocity component

•
• xtrans =

axial coordinate of the transition onset

•
• DN =

number of samples between two independent observations

•
• Mis =

isentropic Mach number

•
• $F=(Fx,Fy,Fz)T$ =

vector of the integral forces on the blade surface

•
• VarN =

sample variance

•
• $u=(u,v,w)T$ =

vector of the Cartesian velocity components

•
• $t=(tx,ty,0)T$ =

tangential vector on the blade surface pointing in the positive mean flow direction

•
• $x=(x,y,z)T$ =

vector of the Cartesian coordinates

Greek Symbols

• α =

flow angle

•
• Δ =

cell size

•
• γ =

heat capacity ratio

•
• ζ =

total pressure loss coefficient

•
• η =

•
• μ =

population mean

•
• ρ =

autocorrelation function

•
• σ =

standard deviation

•
• τ =

lag between samples

•
• $τ$ =

wall shear stress vector

Subscripts

• 0 =

stagnation value

•
• ax =

axial component

•
• inlet, □outlet =

reference values at the inlet and outlet

•
• SS, □PS =

values on the suction side and pressure side

Superscripts

• $□′$ =

fluctuations based on time average

•
• $□+$ =

non-dimensionalized value based on friction velocity

Averaging operators

• $□¯$ =

time average

•
• $□¯tstart→tend$ =

time average from tstart to tend

References

1.
Frey
,
C.
,
Ashcroft
,
G.
,
Kersken
,
H.-P.
,
Schönweitz
,
D.
, and
Mennicken
,
M.
,
2018
, “
Simulation of Indexing and Clocking With Harmonic Balance
,”
Int. J. Turbomach. Propul. Power
,
3
(
1
), p.
1
.
2.
Frey
,
C.
,
Ashcroft
,
G.
,
Kersken
,
H.-P.
, and
Schlüß
,
D.
,
2019
, “
Flutter Analysis of a Transonic Steam Turbine Blade With Frequency and Time-domain Solvers
,”
Int. J. Turbomach. Propul. Power
,
4
(
2
), p.
15
.
3.
Leggett
,
J.
,
Priebe
,
S.
,
Sandberg
,
R.
,
Michelassi
,
V.
, and
Shabbir
,
A.
,
2016
, “
Detailed Investigation of RANS and LES Predictions of Loss Generation in An Axial Compressor Cascade At Off Design Incidences
,”
Proceedings of the ASME Turbo Expo 2016: Turbomachinery Technical Conference and Exposition, Vol. 2A: Turbomachinery
,
Seoul, South Korea
,
June 13–17
, Paper No. GT2016-57972
.
4.
Tucker
,
P.
,
2011
, “
Computation of Unsteady Turbomachinery Flows: Part 1–Progress and Challenges
,”
Prog. Aerosp. Sci.
,
47
(
7
), pp.
522
545
.
5.
Tucker
,
P.
,
2011
, “
Computation of Unsteady Turbomachinery Flows: Part 2—LES and Hybrids
,”
Prog. Aerosp. Sci.
,
47
(
7
), pp.
546
569
.
6.
Tyacke
,
J.
,
Tucker
,
P.
,
Jefferson-Loveday
,
R.
,
,
N.
,
Watson
,
R.
,
Naqavi
,
I.
, and
Yang
,
X.
,
2013
, “
Large Eddy Simulation for Turbines: Methodologies, Cost and Future Outlooks
,”
ASME J. Turbomach.
,
136
(
6
), p.
061009
.
7.
Hillewaert
,
K.
,
Carton de Wiart
,
C.
,
Verheylewegen
,
G.
, and
Arts
,
T.
,
2014
, “
Assessment of a High-Order Discontinuous Galerkin Method for the Direct Numerical Simulation of Transition At Low-Reynolds Number in the T106c High-Lift Low Pressure Turbine Cascade
,”
Proceedings of the ASME Turbo Expo 2014: Turbine Technical Conference and Exposition, Vol. 2B: Turbomachinery
,
Düsseldorf, Germany
,
June 16–20
, Paper No. GT2014-26739
.
8.
Wheeler
,
A. P. S.
,
Sandberg
,
R. D.
,
Sandham
,
N. D.
,
Pichler
,
R.
,
Michelassi
,
V.
, and
,
G.
,
2016
, “
Direct Numerical Simulations of a High-Pressure Turbine Vane
,”
ASME J. Turbomach.
,
138
(
7
), p.
071003
.
9.
Pichler
,
R.
,
Zhao
,
Y.
,
Sandberg
,
R. D.
,
Michelassi
,
V.
,
Pacciani
,
R.
,
Marconcini
,
M.
, and
Arnone
,
A.
,
2018
, “
LES and RANS Analysis of the End-Wall Flow in a Linear LPT Cascade: Part I – Flow and Secondary Vorticity Fields Under Varying Inlet Condition
,”
Proceedings of the ASME Turbo Expo 2018: Turbomachinery Technical Conference and Exposition, Vol. 2B: Turbomachinery of Turbo Expo: Power for Land, Sea, and Air
,
Oslo, Norway
,
June 11–15
, Paper No. GT2018-76233
.
10.
Bergmann
,
M.
,
Morsbach
,
C.
, and
Ashcroft
,
G.
,
2020
, “Assessment of Split Form Nodal Discontinuous Galerkin Schemes for the LES of a Low Pressure Turbine Profile,”
Direct and Large Eddy Simulation XII
,
M.
García-Villalba
,
H.
Kuerten
, and
M. V.
Salvetti
, eds.,
Springer International Publishing
,
Basel
, pp.
365
371
.
11.
Tyacke
,
J. C.
, and
Tucker
,
P. G.
,
2015
, “
Future Use of Large Eddy Simulation in Aeroengines
,”
ASME J. Turbomach.
,
137
(
8
), p.
081005
.
12.
Tyacke
,
J.
,
,
N.
,
Trojak
,
W.
,
Watson
,
R.
,
Ma
,
Y.
, and
Tucker
,
P.
,
2019
, “
Turbomachinery Simulation Challenges and the Future
,”
Prog. Aerosp. Sci.
,
110
(
3
), p.
100554
.
13.
Talnikar
,
C.
,
Blonigan
,
P. J.
,
Bodart
,
J.
, and
Wang
,
Q.
,
2015
,
Optimization With LES – Algorithms for Dealing With Sampling Error of Turbulence Statistics
,
American Institute of Aeronautics and Astronautics
.
14.
Oliver
,
T. A.
,
Malaya
,
N.
,
Ulerich
,
R.
, and
Moser
,
R. D.
,
2014
, “
Estimating Uncertainties in Statistics Computed From Direct Numerical Simulation
,”
Phys. Fluids
,
26
(
3
), p.
035101
.
15.
Trenberth
,
K. E.
,
1984
, “
Some Effects of Finite Sample Size and Persistence on Meteorological Statistics. Part I: Autocorrelations
,”
Mon. Weather Rev.
,
112
(
12
), pp.
2359
2368
.
16.
Ries
,
F.
,
,
K.
,
Dressler
,
L.
,
Janicka
,
J.
, and
,
A.
,
2018
, “
Evaluating Large Eddy Simulation Results Based on Error Analysis
,”
Theor. Comput. Fluid Dyn.
,
32
(
6
), pp.
733
752
.
17.
Mockett
,
C.
,
Knacke
,
T.
, and
Thiele
,
F.
,
2010
, “
Detection of Initial Transient and Estimation of Statistical Error in Time Resolved Turbulent Flow Data
,”
Proceedings of the 8th International Symposium on Engineering Turbulence Modelling and Measurements
,
Marseille, France
,
June 9–11
.
18.
Morsbach
,
C.
, and
Bergmann
,
M.
,
2020
, “Critical Analysis of the Numerical Setup for the Large-Eddy Simulation of the Low-Pressure Turbine Profile T106C,”
Direct and Large Eddy Simulation XII
,
García-Villalba
,
M.
,
Kuerten
,
H.
,
Salvetti
,
M. V.
, eds.,
Springer International Publishing
,
Basel
, pp.
343
348
.
19.
White
,
K. P. J.
,
1997
, “
An Effective Truncation Heuristic for Bias Reduction in Simulation Output
,”
Simulation
,
69
(
6
), pp.
323
334
.
20.
Tennekes
,
H.
, and
Lumley
,
J. L.
,
1972
,
A First Course in Turbulence
,
MIT Press
,
Cambridge, MA
.
21.
Knapp
,
B.
,
Frantal
,
S.
,
Cibena
,
M.
,
Schreiner
,
W.
, and
Bauer
,
P.
,
2011
, “
Is An Intuitive Convergence Definition of Molecular Dynamics Simulations Solely Based on the Root Mean Square Deviation Possible?
,”
J. Comput. Biol.
,
18
(
8
), pp.
997
1005
.
22.
,
K.
,
Robinson
,
S.
, and
Davies
,
R.
,
2010
, “
Automating Warm-Up Length Estimation
,”
J. Oper. Res. Soc.
,
61
(
9
), pp.
1389
1403
.
23.
White
,
K. P. J.
,
Cobb
,
M. J.
, and
Spratt
,
S. C.
,
2000
, “
A Comparison of Five Steady-state Truncation Heuristics for Simulation
,”
Proceedings of the 32nd Conference on Winter Simulation, WSC ’00
,
Orlando, FL
,
Dec. 10–13
,
Society for Computer Simulation International
, pp.
755
760
.
24.
White
,
K. P.
, and
Robinson
,
S.
,
2010
, “
The Problem of the Initial Transient (again), Or why MSER Works
,”
J. Simul.
,
4
(
4
), pp.
268
272
.
25.
Spratt
,
S. C.
,
1998
, “
Heuristics for the Startup Problem
,”
Master’s thesis
,
Department of Systems Engineering
,
University of Virginia
.
26.
Yousefi
,
S.
,
2011
, “
MSER-5Y: An Improved Version of MSER-5 With Automatic Confidence Interval Estimation
,”
Master’s thesis
,
North Carolina State University
.
27.
Donzis
,
D. A.
,
Yeung
,
P. K.
, and
Sreenivasan
,
K. R.
,
2008
, “
Dissipation and Enstrophy in Isotropic Turbulence: Resolution Effects and Scaling in Direct Numerical Simulations
,”
Phys. Fluids
,
20
(
4
), p.
045108
.
28.
Leith
,
C. E.
,
1973
, “
The Standard Error of Time-average Estimates of Climatic Means
,”
J. Appl. Meteorol.
,
12
(
6
), pp.
1066
1069
.
29.
O’Neill
,
P.
,
Nicolaides
,
D.
,
Honnery
,
D.
, and
Soria
,
J.
,
2004
, “
Autocorrelation Functions and the Determination of Integral Length with Reference to Experimental and Numerical Data
,”
Proceedings of the Fifteenth Australasian Fluid Mechanics Conference
,
M.
Behnia
,
W.
Lin
,
G. D.
McBain
, eds., Vol.
1
,
University of Sydney
.
30.
Lepage
,
P.
,
Gohlke
,
C.
, and
Hackett
,
D.
,
2020
, gplepage/gvar: gvar version 11.7.
31.
Benedict
,
L. H.
, and
Gould
,
R. D.
,
1996
, “
Towards Better Uncertainty Estimates for Turbulence Statistics
,”
Exp. Fluids
,
22
(
2
), pp.
129
136
.
32.
Michálek
,
J.
,
Monaldi
,
M.
, and
Arts
,
T.
,
2012
, “
Aerodynamic Performance of a Very High Lift Low Pressure Turbine Airfoil (T106C) At Low Reynolds and High Mach Number with Effect of Free Stream Turbulence Intensity
,”
ASME J. Turbomach.
,
134
(
6
), p.
061009
.
33.
Carton de Wiart
,
C.
,
Hillewaert
,
K.
,
Lorriaux
,
E.
, and
Verheylewegen
,
G.
,
2015
, “
Development of a Discontinuous Galerkin Solver for High Quality Wall-resolved/modelled DNS and LES of Practical Turbomachinery Flows on Fully Unstructured Meshes
,”
Proceedings of the ASME Turbo Expo 2015: Turbine Technical Conference and Exposition, Vol. 2B: Turbomachinery
,
,
June 15–19
, Paper No. GT2015-43428
.
34.
Hillewaert
,
K.
,
Hartmann
,
R.
,
Leicht
,
T.
,
Couaillier
,
V.
,
Wang
,
Z. J.
, and
Cagnone
,
J. S.
,
2016
, “
Summary and Conclusions of the 4th International Workshop on High Order CFD Methods
,”
ECCOMAS Congress 2016 VII European Congress on Computational Methods in Applied Sciences and Engineering
,
Crete Island, Greece
,
June 5–10
.
35.
van Leer
,
B.
,
1979
, “
Towards the Ultimate Conservative Difference Scheme. V. A Second-Order Sequel to Godunov’s Method
,”
J. Comput. Phys.
,
32
(
1
), pp.
101
136
.
36.
Roe
,
P. L.
,
1981
, “
Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes
,”
J. Comput. Phys.
,
43
(
2
), pp.
357
372
.
37.
Nicoud
,
F.
, and
Ducros
,
F.
,
1999
, “
Subgrid-Scale Stress Modelling Based on the Square of the Velocity Gradient Tensor
,”
Flow Turbul. Combust.
,
62
(
3
), pp.
183
200
.
38.
Schlüß
,
D.
,
Frey
,
C.
, and
Ashcroft
,
G.
,
2016
, “
Consistent Non-Reflecting Boundary Conditions for Both Steady and Unsteady Flow Simulations in Turbomachinery Applications
,”
ECCOMAS Congress 2016 VII European Congress on Computational Methods in Applied Sciences and Engineering
,
Crete Island, Greece
,
June 5–10
, pp. 7403–7422
.
39.
Weber
,
A.
, and
Sauer
,
M.
,
2016
,
Pymesh – Template Documentation. Technical report DLR-IB-AT-KP-2016-34, German Aerospace Center (DLR), Institute of Propulsion Technology, Linder Hoehe, Cologne, Germany
.
40.
Bergmann
,
M.
,
Gölden
,
R.
, and
Morsbach
,
C.
,
2018
, “
Numerical Investigation of Split Form Nodal Discontinuous Galerkin Schemes for the Implicit LES of a Turbulence Channel Flow
,”
Proceedings of the 7th European Conference on Computational Fluid Dynamics
,
Glasgow, UK
,
June 11–15
.