Sensitivity analysis of thermoacoustic instability with adjoint Helmholtz solvers

Gas turbines and rocket engines sometimes suffer from violent oscillations caused by feedback between acoustic waves and flames in the combustion chamber. These are known as thermoacoustic oscillations and they often occur late in the design process. Their elimination usually requires expensive tests and redesign. Full scale tests and laboratory scale experiments show that these oscillations can usually be stabilized by making small changes to the system. The complication is that, while there is often just one unstable natural oscillation (eigenmode), there are many possible changes to the system. The challenge is to identify the optimal change systematically, cheaply, and accurately. This paper shows how to evaluate the sensitivities of a thermoacoustic eigenmode to all possible system changes with a single calculation by applying adjoint methods to a thermoacoustic Helmholtz solver. These sensitivities are calculated here with finite-difference and finite-element methods, in the weak form and the strong form, with the discrete adjoint and the continuous adjoint, and with a Newton method applied to a nonlinear eigenvalue problem and an iterative method applied to a linear eigenvalue problem. This is a detailed comparison of adjoint methods applied to thermoacoustic Helmholtz solvers. matlab codes are provided for all methods and all figures so that the techniques can be easily applied and tested. This paper explains why the finite difference of the strong form equations with replacement boundary conditions should be avoided and why all of the other methods work well. Of the other methods, the discrete adjoint of the weak form equations is the easiest method to implement; it can use any discretization and the boundary conditions are straightforward. The continuous adjoint is relatively easy to implement but requires careful attention to boundary conditions. The summation by parts finite difference of the strong form equations with a simultaneous approximation term for the boundary conditions is more challenging to implement, particularly at high order or on nonuniform grids. Physical interpretation of these results shows that the well-known Rayleigh criterion should be revised for a linear analysis. This criterion states that thermoacoustic oscillations will grow if heat release rate oscillations are sufficiently in phase with pressure oscillations. In fact, the criterion should contain the adjoint pressure rather than the pressure. In self-adjoint systems the two are equivalent. In non-self-adjoint systems, such as all but a special case of thermoacoustic systems, the two are different. Finally, the sensitivities of the growth rate of oscillations to placement of a hot or cold mesh are calculated, simply by multiplying the feedback sensitivities by a number. These sensitivities are compared successfully with experimental results. With the same technique, the influence of the viscous and thermal acoustic boundary layers is found to be negligible, while the influence of a Helmholtz resonator is found, as expected, to be considerable.

This article appears in the following collection:

Figure 1

Distributions of density $\overline{\rho }\left(x\right)$, heat release region $h\left(x\right)$, and measurement region $w_{\overline{\rho }}$ for (a) the Rijke tube in [33] and (b) a simplified model of a ${\mathrm{H}}_2/\mathrm{LOx}$ rocket engine. The dimensional, reference, and nondimensional quantities are in Table 1.

Figure 2

Direct pressure $p$ and velocity $u$ eigenfunctions calculated with four different spatial discretizations (FDS, FEW, FDW, and SBP ) and two iteration procedures (nonlin and linear) for the (a) Rijke tube and (b) model rocket engine, for $N=100$. All eight lines lie on top of each other and all eigenvalues are the same to the requested tolerances. This figure was created with Fig_002.m. The agreement between the methods can be improved by increasing $N$.

Figure 3

Adjoint pressure $p^{†}$ eigenfunctions calculated with four different spatial discretizations (FDS, FEW, FDW, and SBP) with either the discrete adjoint method DA or the continuous adjoint method CA for the (a) Rijke tube and (b) model rocket engine, for $N=100$. The bottom two frames show a close-up around the boundaries, revealing oscillations in the discrete adjoint of the finite difference of the strong form equations with replacement boundary conditions (FDS_DA). For the rocket engine, the CA and DA eigenfunctions and eigenvalues differ slightly from each other because $w_{\overline{\rho }}$ is not quite zero at the boundary (order ${10}^{-3}$), so the neglected term in (24) has a slight influence. This figure was created with Fig_003.m. By running this with larger $N$, one can check that the eigenvalues are the same to the requested tolerance.

Figure 4

Absolute values of the base state sensitivities ($\partial s/\partial \tau , \partial s/\partial h, \partial s/\partial w_{\overline{\rho }}, \partial s/\partial \overline{v}, \partial s/\partial n, \partial s/\partial k_u$, and $\partial s/\partial k_d$) calculated with the discrete adjoint method (DA) and continuous adjoint method (CA) for four spatial discretizations (FDS, FEW, FDW, and SBP) and for the (a) Rijke tube and (b) model rocket engine, for $N=100$. All lines except FDS_DA lie on top of each other and all sensitivities to $n, k_u$, and $k_d$ are the same to the requested precision. This figure was created with Fig_004.m. Higher precision can be achieved by increasing $N$.

Figure 5

Absolute values of the feedback sensitivities ($\partial s/\partial {\dot{m}}_{\overline{\rho },u}, \partial s/\partial {\dot{m}}_{\overline{\rho },p}, \partial s/\partial {\mathbf{f}}_{\overline{\rho },u}, \partial s/\partial {\mathbf{f}}_{\overline{\rho },p}, \partial s/\partial {\dot{q}}_{\overline{p},u}$, and $\partial s/\partial {\dot{q}}_{\overline{p},p}$) calculated with the discrete adjoint method (DA) and continuous adjoint method (CA) for four spatial discretizations (FDS, FEW, FDW, and SBP) for the (a) Rijke tube and (b) model rocket engine, for $N=100$. All lines lie on top of each other except FDS_DA, which oscillates. These results were created with Fig_005.m.

Figure 6

Absolute values of the receptivity of the eigenvalue to periodic injection of ${\dot{m}}_{\overline{\rho }}\left(x\right)e^{s_{f}t}$ into the mass equation (2c), ${\mathbf{f}}_{\overline{\rho }}\left(x\right)e^{s_{f}t}$ into the momentum equation (2d), and ${\dot{q}}_{\overline{p}}\left(x\right)e^{s_{f}t}$ into the energy equation (2e), where $s_f$ is an imaginary number equal to the angular frequency of forcing. From (6), the receptivity of the eigenvalue to mass injection into the mass equation is ${\dot{m}}_{\overline{\rho }}^{†}/\overline{\rho }$, momentum injection into the momentum equation is ${\mathbf{f}}_{\overline{\rho }}^{†}/\overline{\rho }$, and heat injection into the energy equation is ${\dot{q}}_{\overline{p}}^{†}/\overline{p}$. These are calculated with the discrete adjoint method (DA) and continuous adjoint method (CA) for four spatial discretizations (FDS, FEW, FDW, and SBP) for the (a) Rijke tube and (b) model rocket engine, for $N=100$. All lines lie on top of each other except FDS_DA, which oscillates. This figure was created with Fig_006.m.

Figure 7

Real and imaginary components of the feedback sensitivities in Fig. 5, calculated with discrete adjoint (DA) of the finite-element method (FEW) for the (a) Rijke tube and (b) model rocket engine, for $N=100$. The wider lines show the same sensitivities calculated when the adjoint pressure $p^{†}$ is replaced by the direction pressure $p$, i.e., as if the system were self-adjoint. There is little difference for the Rijke tube, which is weakly non-self-adjoint, but a large difference for the rocket engine, which is strongly non-self-adjoint. This figure was created with Fig_007.m.

Figure 8

Real and imaginary components of the base state sensitivities in Fig. 4, calculated with discrete adjoint (DA) of the finite-element method (FEW) for the (a) Rijke tube and (b) model rocket engine, for $N=100$. As in Fig. 7, the wider lines show the base state sensitivities formed when the adjoint pressure $p^{†}$ is replaced by the direction pressure $p$, i.e., if the system is assumed to be self-adjoint. This figure was created with Fig_008.m.

Figure 9

Growth rates (color contours) and nondimensional frequencies (black line contours) for the Rijke tube with a variable-diameter iris placed at the downstream boundary and a variable heater position $(X_h,X_w)$. A quasi-Newton algorithm uses adjoint-based gradient information (white arrows) to converge to the point with maximum growth rate at a user-defined frequency, in this case 3.2. This figure was created with Fig_009.m.

Figure 10

Diagram of a Rijke tube with perimeter ${\mathrm{\Gamma }}_c$ and cross-sectional area $S_c$ containing a passive feedback device at position $X$ with width $\delta X$.

Figure 11

Growth rate shift caused by drag from an adiabatic mesh placed at position $X_m$.

Figure 12

Growth rate shift caused by drag from a secondary hot mesh placed at position $X_m$. The offset is caused by a change to the base state induced by the secondary hot mesh [43].

Figure 13

Growth rate shift caused by the viscous acoustic boundary layer, estimated with the feedback sensitivity.

Figure 14

Growth rate shift caused by the thermal acoustic boundary layer, estimated with the feedback sensitivity.

Figure 15

Growth rate shift of an oscillation at $\omega =1160\mathrm{rad}{\mathrm{s}}^{-1}$ due to a Helmholtz resonator with natural frequency ${\omega }_R=1020\mathrm{rad}{\mathrm{s}}^{-1}$, as a function of the resonator position $X_n$.