Uncertainty quantification and sensitivity analysis of thermoacoustic stability with non-intrusive polynomial chaos expansion

Abstract

In this paper, non-intrusive polynomial chaos expansion (NIPCE) is used for forward uncertainty quantification and sensitivity analysis of thermoacoustic stability of two premixed flame configurations. The first configuration is a turbulent swirl combustor, modeled by the Helmholtz equation with an $n-\tau$ flame model. Uncertain input parameters are the gain and the time delay of the flame, as well as the magnitude and the phase of the outlet reflection coefficient. NIPCE is successfully validated against Monte Carlo simulation. It is observed that the first order expansion suffices to yield accurate results. The second configuration under investigation is a low order network model of a laminar slit burner, with the flame transfer function identified from weakly compressible CFD simulations of laminar reacting flow. Firstly the uncertainty and sensitivity of the growth rate due to three uncertain input parameters of the CFD model – i.e., flow velocity, burner plate temperature and equivalence ratio – are analyzed. A Monte Carlo simulation is no longer possible due to the computational cost of the CFD simulations. Secondly, two additional uncertain parameters are taken into account, i.e., the respective magnitudes of inlet and outlet reflection coefficients. This extension of the analysis does not entail a considerable increase in computational cost, since the additional parameters are included only in the low order network model. In both cases, the second order expansion is sufficient to model the uncertainties in growth rate.

Introduction

Thermoacoustics deals with a coupling between acoustics and heat release. This coupling may lead to a self-induced instability with excessive fluctuations in pressure, velocity, and temperature. Such instabilities occur, e.g., in gas turbines or rocket engines, where they can lead to a catastrophic system failure. To guarantee safety margins and normal operating conditions, the thermoacoustic behavior of a combustor should be studied. Indeed, a comprehensive thermoacoustic stability analysis is an important part of combustor design.

In computational analysis, boundary and operating conditions as well as model parameters of a system model are in general not known exactly, but instead are uncertain. Uncertainties propagate through the system model and affect the prediction of quantities of interest, making them uncertain. Forward uncertainty quantification (UQ) strives to characterize in a quantitative manner the impact of uncertain input or model parameters on the reliability of model predictions. An important aspect of the study of uncertainties is sensitivity analysis, which investigates the influence of each uncertain parameter on a quantity of interest. Sensitivity analysis helps to identify the most important uncertain parameters, which should be accounted for. Due to the fast development of data-driven methodologies in recent years, UQ and sensitivity analysis are becoming important topics in all engineering fields [1], [2], [3], [4], [5].

The main task of UQ in linear thermoacoustics stability analysis is to investigate how uncertainties in geometry, operating and boundary conditions as well as modeling parameters affect the growth rates of the thermoacoustics eigenmodes. Despite the fact that thermoacoustic systems are in general very sensitive to such uncertainties, until now only a small number of studies have been devoted to this topic, which shall be reviewed briefly in the following.

The most common and straightforward UQ method is Monte Carlo simulation [6]. The method numerically generates random samples of uncertain input parameters. For each sample, the quantity of interest is computed. The ensemble of obtained results is assumed to faithfully represent the variability of the quantity of interest. Monte Carlo simulation requires a large number of samples and is feasible only for system models with fast evaluation. If a single evaluation is computationally expensive, more sophisticated UQ methods are required. Fundamentally, there are two strategies to overcome a high computational cost: (1) simplify the original system model to reduce the evaluation time, or (2) lower the number of samples required to achieve converged statistics. An example of the latter approach are smart sampling techniques, such as Latin hypercube sampling [7].

Bauerheim et al. [8] investigated a simplistic model of an annular combustor with 19 burners using a Helmholtz equation coupled with $n-\tau$ flame models. The uncertain input parameters were gain n and time lag τ of each single flame, resulting in a total of 38 uncertain parameters. The authors applied an active subspace approach [9] that reduced the 38 input parameters to three active variables. Fifty samples, i.e., model evaluations, were required to identify the active variables. Linear, quadratic and cubic reduced order algebraic models were fitted, using a few dozen samples. Finally, a Monte Carlo simulation was performed on the reduced order models to quantify the risk factor, i.e., the probability of an unstable state. The quadratic and cubic models showed accurate results in comparison to the original model. The active subspace approach is well suited for problems with a large number of uncertain parameters, since it eliminates the curse of dimensionality. However, in order to firstly identify active variables, and secondly fit a surrogate model, the method requires several tens or even hundreds of system evaluations.

Ndiaye et al. [10] also coupled a Helmholtz equation solver with an $n-\tau$ flame model in order to assess the thermoacoustic stability of a single premixed swirled burner. Again gain n and time lag τ were considered as uncertain parameters. Relying on a multiple linear regression technique, a bilinear algebraic surrogate model was tuned to reproduce growth rates of the dominant thermoacoustic mode, as predicted by the full model under variation of the uncertain parameters. Then the surrogate model was used to compute the risk factor with a Monte Carlo approach. The authors carefully assessed the number of full model evaluations required to tune the surrogate model and conclude that ten full model evaluations suffice to accurately estimate the risk factor.

Silva et al. [11], [12] investigated the same swirled burner as Ndiaye et al. [10], considering four uncertain parameters, i.e., gain n and time lag τ of the flame model as well as magnitude and phase of the outlet reflection coefficient. Direct and adjoint eigenvectors were used to construct a first and second order expansions of the nonlinear eigenvalue problem around a reference eigenvalue. Using these expansions, deviations from the reference eigenvalue were computed at reduced computational cost. The results of the second order expansion were in good agreement with the ones obtained by solving the nonlinear eigenvalue problem. Recently, the same burner was investigated by Mensah et al. [13]. The authors used a flame transfer function fitted from the experiment instead of a n–τ flame model. The uncertainty in the flame transfer function was modeled with two parameters, relative error in gain and absolute error in phase, which were assumed to be constant in the entire frequency range. Magri et al. [14] used the adjoint approach in combination with the active subspace approach to compute the risk factor of the annular combustor investigated previously by Bauerheim et al. [8].

With the adjoint method one can easily build a reduced order model and compute sensitivities. One should keep in mind, however, that the accuracy of results will decrease with increasing variances in the input uncertainties, since the adjoint method constructs an expansion around a single reference point. Furthermore, the commonly used matrix-based method requires access to the state matrices, which may require substantial effort, or may – e.g., for a commercial CFD solvers – not be possible at all. In that case one could use matrix-free methods, as proposed by Waugh et al. [15].

In the present study, non-intrusive polynomial chaos expansion (NIPCE) is used for UQ in linear thermoacoustic stability analysis. The polynomial chaos expansion approximates uncertain output as a polynomial function of uncertain input parameters and provides statistical moments of the output. The advantage of the non-intrusive variant of the method is the treatment of an investigated system as a “black box”, so the method can be applied to any system model without any special code modifications. Another advantage is that NIPCE provides a polynomial approximation of the output quantity, which can be used as a reduced order model, or for sensitivity analysis.

NIPCE has hitherto not been used in thermoacoustics – except for the analysis of Nair et al. [16], which is discussed below – but it has been applied frequently in computational fluid dynamics [3] and also in simulations of reacting flows [17]. There are several possibilities for constructing the NIPCE: Hosder et al. [18], [19], [20] used linear regression, Reagan et al. [1], [21] used the Latin hypercube sampling technique, Tritschler et al. [22] used the Gauss quadrature. We choose the Gauss quadrature, since it synchronizes the expansion order with the number of quadrature points. Hence, a large number of system evaluations leads to a higher expansion order.

Nair et al. [16] explored the use of NIPCE to quantify the uncertainty of subcritical Hopf bifurcations predicted by a simplistic model of a Rijke tube. However, a subcritical bifurcation represents a special challenge, since it exhibits a discontinuity, which with standard polynomials cannot represent well. Indeed, NIPCE was originally developed for continuous uncertainties, thus we use this method for continuous problems. Instead of NIPCE, Nair and co-workers eventually employed a sampling technique based on equal probabilities to compute the failure probability for several input parameter values. Then they reconstructed the response in the entire parameter space by an interpolation technique that captures discontinuities. The UQ was then performed by a Monte Carlo simulation of the interpolated model.

We have discussed briefly the use of active subspace and adjoints for UQ in thermoacoustics. NIPCE should not be seen as a competitor to the above mentioned methods, but rather as a complementary tool to solve problems, where neither active subspace nor adjoints are feasible. For instance, if access to solver state matrices is not feasible, then the adjoint approach cannot be employed, unless the matrix-free methods are used. On the contrary, NIPCE can easily be applied in that case, since it treats the solver as a “black box”. If a single system evaluation is so expensive that the random sampling required for identification of active variables is impractical, then NIPCE with Gauss quadrature may still solve the problem, since the method requires only a few quadrature points. On the other hand, the application of the NIPCE with standard polynomials is constrained to continuous uncertainties, and the number of uncertain parameters should not exceed roughly ten. For discontinuous or steep functions one may try NIPCE with Haar polynomials [23] or the above-mentioned method proposed by Nair et al. [16].

In this study, we employ NIPCE to perform UQ and sensitivity analysis for thermoacoustic systems. Moreover, it is shown that probability density functions of the sensitivities are suitable for analysis of systems with multiple uncertain input parameters. The paper is structured as follows: in the next section, the fundamentals of the NIPCE method are described briefly. Then the NIPCE is validated against a Monte Carlo simulation for a simplified combustor, which is modeled by the Helmholtz equation with the n–τ flame model. The uncertain parameters are the gain n and the time delay τ, as well as the magnitude and phase of an outlet reflection coefficient. Afterward, a thermoacoustic network model is investigated using NIPCE. The network model includes flame transfer functions identified from CFD simulations. In this case, validation against the Monte Carlo simulation is no longer possible, since a single CFD simulation for the flame identification takes several hours. Therefore this is a typical application case for NIPCE. The uncertain input parameters are the flow velocity, the burner plate temperature, and the equivalence ratio. In the next part of this study, we show that the uncertainty analysis can be extended by additional uncertain parameters without significant computational effort, provided that these parameters are included only in the network model and independent of the other uncertain parameters, in particular the flame transfer function. In the present study, two additional uncertain parameters, i.e., the magnitude of the inlet and outlet reflection coefficients, are introduced. A summary of results concludes the paper.