H ∞ loop-shaping control of azimuthal combustion instabilities in annular combustors

Abstract

Lean combustion technique and annular geometries are preferred in aero-engines and gas turbines, which however may lead to azimuthal combustion instabilities. Active control can be used to stabilize combustion instabilities. Due to its easy use, linear feedback controllers embedded with linear flame response models under weak perturbation amplitude are typically preferred. However, flame responses to oncoming disturbances are typically nonlinear; such controllers are not guaranteed to stabilize. Model-based control strategies generally focus on axisymmetric cases, even though symmetry breaking of azimuthal thermoacoustic modes often occurs in annular combustors. This work uses an improved thermoacoustic model to simulate combustion instabilities within annular combustors, providing a platform on which control strategies development can be performed. The improved model takes into account the nonlinear flame response and the symmetry breaking of azimuthal modes. Single-input single-output control strategies targeting on these nonlinear instabilities are developed in this work. Such controllers can achieve stability for linear and nonlinear fluctuations as well as in symmetric and non-axisymmetric cases. The controllers adopt the $H_{\infty }$ loop-shaping control strategy and a satisfactory robust performance is obtained.

Introduction

In order to reduce ${\mathrm{NO}}_{\mathrm{x}}$ emissions, modern gas turbines and aero-engines are desirable to operate under lean premixed combustion conditions. However, lean premixed systems are highly susceptible to undesirable combustion instabilities caused by a coupling between unsteady heat release and acoustic waves within the combustion chamber. Because of this coupling, the pressure and heat release oscillation amplitudes may successively grow, until nonlinear effects cause saturation into limit cycle oscillations. ¹ These instabilities always lead to significant noise and structural damages. ² To suppress combustion instabilities, the coupling between the acoustic waves and the unsteady heat release must be interrupted. Available methods for this aim may be achieved using either passive or active control.^3–6 Active control technology combines combustion instabilities with advanced control theories, so it has the characteristics of flexibility, wide application, and high efficiency. The open-loop transfer function (OLTF) represents the dynamic response of combustion system to external actuations, which is required before the controller design process.^7,8

In order to achieve uniform temperature distribution at the turbine inlet and combustion chamber geometry compactness, modern gas turbine combustion chambers usually have an annular geometry with multiple azimuthally arranged burners. ⁹ The largest size of the annular combustor is often the circumference, and the lowest acoustic resonance frequency is then associated with modes in the azimuthal direction. Annular combustors often suffer from azimuthal thermoacoustic instability, ¹⁰ and its model-based control is then more challenging than that for longitudinal modes because of the high complexity of its physical process and the modeling of annular combustors.

The design of traditional linear feedback controllers for annular combustion chambers is typically based on the OLTF corresponding to small and linear perturbations. ¹¹ However, such feedback controllers are always activated from within nonlinear limit cycle oscillations since the dynamics of real unstable combustors become dominated by nonlinear mechanisms once perturbation amplitudes grow sufficiently, during which the linear OLTF will be invalid. So the effectiveness of these controllers is not guaranteed. ⁷ There has been no systematic research work to solve this problem until now.

When the flame response inside multiple burners is different, the rotating symmetry of azimuthal mode will be broken, ¹² and such mode is defined as the non-axisymmetric case in this study. The case that the flame responses inside multiple burners are identical is considered axisymmetric. The design of the controller is a big challenge when the azimuthal mode of the annular combustor exhibits axial asymmetry. A multiple-input-multiple-output feedback controller has been designed for modal annular combustors with different burner types. ¹¹ In practice, it would be highly desirable to reduce the number of feedback loops and sensors/actuators as far as possible. When considering nonlinear flame responses, the form of symmetry breaking becomes more complicated. It is thus much more attractive to develop a single-input single-output (SISO) controller to effectively suppress the unstable mode exhibiting asymmetry in annular combustors.

This work improves the annular combustor model developed by Bauerheim et al. ¹³ The modified model takes into account the effects of flame nonlinearity, bulk fluid flow, and viscothermal damping inside the burners which are ignored in the original model for the sake of simplicity. The model combustor utilizes a pressure transducer as the controller sensor, and a valve on the fuel feed line which can modulate the fuel flow rate of all burners simultaneously and equally is used as the actuator.

In this work, $H_{\infty }$ loop-shaping controllers are designed in order to guarantee stabilizations for all cases, from small linear fluctuations to large limit cycle oscillations, and from the axisymmetric cases to the non-axisymmetric cases. The modern robust control concept of an $\nu \text{-gap}$ among open-loop plants is the most helpful measure of distance between them. ⁷ The $\nu \text{-gap}$ metric provides a bound on the minimum required “robustness stability margin” for $H_{\infty }$ loop-shaping controller. ¹⁴ $H_{\infty }$ loop-shaping, one of the modern robust control methodologies, has been successfully applied to a variety of fields, ¹⁵ such as control of combustion instabilities,^16,17 voltage converter design, ¹⁸ and flight control system. ¹⁹ The $H_{\infty }$ loop-shaping design is a combination of loop-shaping and modern $H_{\infty }$ optimization which can guarantee the robust stability and desired performance of closed-loop controller simultaneously. ²⁰ It would be highly desirable to pursue low-order controllers in practice, although $H_{\infty }$ loop-shaping can be readily applied to plants with high dimensions. ⁷ In this work, a low-order rational transfer function is used to approximate the OLTF for the sake of controller design since the controllers generally have orders comparable to those of the OLTF.

The rest of this paper is arranged as follows. The next section describes the annular combustor model used in this work, including the principle of the network reduction methodology and nonlinear flame transfer function (NFTF). Analysis of the azimuthal thermoacoustic modes is carried for both the axisymmetric and the non-axisymmetric cases in the “Analysis of the unstable azimuthal thermoacoustic mode” section. The “Feedback control of azimuthal combustion instabilities” section presents the design processes of $H_{\infty }$ loop-shaping controllers which are guaranteed to be stabilizing from within limit cycle oscillations for the axisymmetric and the non-axisymmetric cases. A standard high-order controller and a low-order controller are obtained in this section. Conclusions are given in the final section.

Combustor model

Description of the network model

A thermoacoustic network model for PBC (plenum+burners+chamber) configuration ¹³ is used to develop model-based feedback controller strategies. This work improves the previous model by taking into account the flame nonlinear response, bulk fluid flow, and thermal viscosity inside the burners.

The model of the annular combustor is shown schematically in Figure 1(a). The two ends of the burners are, respectively, connected to the annular plenum and chamber. Closed-end boundary conditions are used here at the inlet of the plenum and the outlet of the chamber. The half-perimeter and cross-sectional area of the cavities are represented by $L=\pi R$ and $S$ , respectively (we assume that the plenum and the chamber have the same dimensions). $L_k$ and $S_k$ denote the length and cross-section area of the kth burner. (Note that each burner has the same sizes.) The azimuthal position is given by the angle θ associated with an abscissa $x=R\theta$ for the plenum and chamber. The position of the flame in the kth burners is indicated by the normalized abscissa $\alpha =y_{F,k}/L_k$ . Mean pressure in the annular cavities is denoted by $\overline{p}$ ( ${\overline{p}}_u$ in the plenum and ${\overline{p}}_b$ in the chamber). The subscripts u and b stand for unburned and burned gases respectively. The flow can be taken to be composed of a steady uniform mean flow (denoted by overbars) and a small perturbation (denoted by primes). For example, the pressure, $p=\overline{p}+p^{\prime }$ .

Figure 1.

Schematic diagram of (a) the PBC configuration and (b) annular network reduction. Source: taken from Bauerheim et al. ¹³

The model developed in Bauerheim et al. ¹³ did not account for the bulk fluid flow and any damping mechanisms. According to mass conservation equation and equation (3), these assumptions are acceptable for annular cavities since their sections are relatively large. However, the cross-section area of the burners can be very small in some real combustors, thus the thermal viscosity and bulk fluid flow inside the burners should be considered.

Following the approach given in Bauerheim et al., ¹³ the full configuration is split into M sectors (Figure 1(b)). Then, for each individual sector, the acoustic problem is separated into three parts:

Propagation in the annular cavities (plenum and chamber) can be described by a 4 × 4 rotation matrix $R_A$ , similar as equation (23) in Bauerheim et al. ¹³

Propagation in the kth burner without passing through the flame:

The cross-section area of the burners is much smaller than that of the annular cavities (including the plenum and the chamber). It should be noted that the viscothermal damping effect should be accounted for if the burner channel radius is smaller than 5 cm for room temperature. ⁵ Therefore, the bulk fluid flow and viscothermal damping effects inside the burners are taken into consideration. The pressure perturbation $p_k^{\prime }$ and the longitudinal velocity perturbation $v_k^{\prime }$ inside the kth burner can be written as

$$p_k^{\prime }(y,t)=\left(A_{\ell ,k}{e}^{-s_{\ell }\frac{y}{{\overline{c}}_{\ell }+{\overline{v}}_{\ell }}}+B_{\ell ,k}{e}^{s_{\ell }\frac{y}{{\overline{c}}_{\ell }-{\overline{v}}_{\ell }}}\right)e^{s_{\ell }t}$$

(1)

and

$${\overline{\rho }}_{\ell }{\overline{c}}_{\ell }{v}_k^{\prime }(y,t)=\left(A_{\ell ,k}{e}^{-s_{\ell }\frac{y}{{\overline{c}}_{\ell }+{\overline{v}}_{\ell }}}-B_{\ell ,k}{e}^{s_{\ell }\frac{y}{{\overline{c}}_{\ell }-{\overline{v}}_{\ell }}}\right)e^{s_{\ell }t}$$

(2)

where $\rho$ and $c$ denote the density and the speed of sound, respectively. $A_{\ell ,k}$ and $B_{\ell ,k}$ are the amplitudes of the upward and downward propagating acoustic waves, respectively, where the subscript $\ell$ can be u or b, meaning the propagation before or after the flames. $s_{\ell }$ is complex frequency. $s_{\ell }=s+\Delta s_{\ell }$ , where $\Delta s_{\ell }$ is the complex frequency correction for the viscothermal damping effects, and $\Delta s_{\ell }={\Lambda }_{\ell }\sqrt{s}$ , with the coefficient ${\Lambda }_{\ell }$ expressed as ⁷

$${\Lambda }_{\ell }={\left(\frac{\pi {\mu }_{\ell }}{{\overline{\rho }}_{\ell }{S}_k}\right)}^{1/2}\left(1+\frac{\gamma -1}{P_r{\hspace{0.17em}}^{1/2}}\right)$$

(3)

where $\gamma$ , ${\mu }_{\ell }$ , and $P_r$ are the ratio of specific heats, dynamic viscosity, and Prandtl number, respectively. The values of these parameters can be found from the appendix of Li and Morgans. ⁷

The acoustic propagation relation in the upstream and downstream regions of the flame inside the kth burner can be described by a 2 × 2 matrix $R_B(s_{\ell }\Delta y/{\overline{c}}_{\ell })$ (see Appendix 1 for the form of $R_B(s_{\ell }\Delta y/{\overline{c}}_{\ell })$ ).

3. Coupling between the plenum and the chamber via the burners:

As shown in Figure 1(b), an H-shaped junction is used to link the acoustic states of the input (corresponding to the end of the sector k−1) and the output (located at the beginning of the sector k). This relation can be expressed as

$${\mathrm{X}}_k(x=0,t)=H_k\cdot {\mathrm{X}}_{k-1}(x=2L/M,t)$$

(4)

where ${\mathrm{X}}_k(x,t)$ stands for the acoustic disturbances at location x in the kth sector: ${\mathrm{X}}_k(x,t)={[\hspace{0.17em}{\mathit{p}}_{\mathit{P}}^{\prime },{\overline{\rho }}_u{\overline{c}}_u{u}_P^{\prime },p_C^{\prime },{\overline{\rho }}_b{\overline{c}}_b{u}_C^{\prime }]}_k^T$ . The expression and its derivation of matrix $H_k$ can be found in Bauerheim et al. ¹³ Inside each burner, the jump condition across the flame is updated as the bulk fluid flow in the burners is taken into account, with more details presented as follows.

The 1D Euler equations (continuity, momentum, and energy equation) in the kth burner have the form ¹

$${[\rho v]}_{y_{F^{-}},k}^{y_{F^{+}},k}=0$$

(5)

$${[p+\rho v^2]}_{y_{F^{-}},k}^{y_{F^{+}},k}=0$$

(6)

$${\left[\left(c_{p}T+\frac{1}{2}{v}^2\right)\rho v\right]}_{y_{F^{-}},k}^{y_{F^{+}},k}=q_k$$

(7)

The subscripts $F^{-}$ and $F^{+}$ indicate the locations immediately upstream and downstream the flame, and $q_k$ is the heat release rate per unit area of the kth flame. Linearizing equations (5) to (7), and substituting the NFTF ${\Gamma }_k$ into $q_k^{\prime }$ (see equation (18)), the jump condition is then obtained

$$\left[\begin{array}{c}{p}_k^{\prime }(y_{F^{-}})\\ v_k^{\prime }(y_{F^{+}})\end{array}\right]={\vartheta }_k\left[\begin{array}{c}{p}_k^{\prime }(y_{F^{-}})\\ v_k^{\prime }(y_{F^{+}})\end{array}\right]$$

(8)

The detailed derivation process of the jump condition and the exact form of ${\vartheta }_k$ are presented in Appendix 2.

The transfer matrix of the kth sector is $R_A(2L{s}_{\ell }/M{\overline{c}}_{\ell })\cdot H_k$ . Connecting M sectors and making use of the periodicity condition, one gets the dispersion relation for the entire system

$$\det \left({\displaystyle \prod _{k=M}^1{R}_k\left(\frac{2L}{M}\right)H_k-I_4}\right)=0$$

(9)

where $I_4$ is a 4 × 4 identity matrix. Equation (9) is an implicit equation for the complex frequency s with the real part being the growth rate and the imaginary part being the frequency. Note that ${\Gamma }_k$ (k = 1−M) is included in equation (9) but has not been given yet. The flame transfer function used in this work is described in the next section.

The NFTF

In this work, the flame nonlinear response is considered through a gain/time delay flame model with saturation bounds. The mean heat release per unit area for the burner k is denoted by ${\overline{q}}_k$ . The fractional heat release fluctuation at the kth burner, $q_k^{\prime }(t)/{\overline{q}}_k$ , is related to the fractional equivalence ratio perturbation within that burner, ${\varphi }_k^{\prime }(t)/\overline{\varphi }$ , using the nonlinear flame model^21,22

$$\frac{q_k^{\prime }(t)}{{\overline{q}}_k}=\{\begin{array}{l}{n}_{\varphi }\frac{{\varphi }_k^{\prime }(t-{\tau }_{\varphi })}{{\overline{\varphi }}_k},\text{if}\left|q_k^{\prime }(t)/{\overline{q}}_k\right|<{\alpha }_{\varphi }\hfill \\ {\alpha }_{\varphi }\text{sign}(q_k^{\prime }(t)),\text{otherwise}\hfill \end{array}$$

(10)

where $n_{\varphi }$ represents the gain at low perturbation amplitude stage, ${\tau }_{\varphi }$ represents the time delay from fuel input to the combustion, and ${\alpha }_{\varphi }$ represents the saturation bound of the fractional change in heat release rate. The saturation-bound nonlinear flame model is a very classical describing model in nonlinearity control problem and is chosen in the current research.

The definition of $\varphi$ is the mass fuel-to-air ratio $(\hspace{0.17em}f/a)$ relative to the fuel-to-air ratio at stoichiometric conditions. Expressing it as mean and fluctuating components, and linearizing it, one then finds that

$$\frac{{\varphi }_k^{\prime }}{{\overline{\varphi }}_k}\text{=}\frac{f_k^{\prime }}{{\overline{f}}_k}-\frac{a_k^{\prime }}{{\overline{a}}_k}$$

(11)

The mass flow of air immediately before the flame location inside the kth burner can be linearized to obtain

$$\frac{a_k^{\prime }}{{\overline{a}}_u}=\frac{v_{F^{-},k}^{\prime }}{{\overline{v}}_u}+\frac{{\rho }_{F^{-},k}^{\prime }}{{\overline{\rho }}_u}$$

(12)

$v_{F^{-},k}^{\prime }/{\overline{v}}_u$ is always much larger than ${\rho }_{F^{-},k}^{\prime }/{\overline{\rho }}_u$ in practice. ²³ Using this assumption and combining equations (11) and (12) give

$$\frac{{\varphi }_k^{\prime }}{\overline{\varphi }}\text{=}\frac{f_k^{\prime }}{{\overline{f}}_u}-\frac{v_{F^{-},k}^{\prime }}{{\overline{v}}_u}$$

(13)

The mean heat release rate per unit area of the kth flame is given by ²⁴

$${\overline{q}}_k=\eta \overline{\varphi }\Delta H_{\mathit{st}}\frac{\mathit{FA}{R}_{\mathit{st}}}{1+\overline{\varphi }FA{R}_{\mathit{st}}}{\overline{\rho }}_u{\overline{v}}_u$$

(14)

In equation (14), $\eta$ is the combustion efficiency (it is assumed herein that all flames have the same combustion efficiency), $\mathit{FA}{R}_{\mathit{st}}$ is the stoichiometric fuel-to-air equivalence ratio, and $\Delta H_{\mathit{st}}$ is the calorific value of the fuel. Substituting equations (13) and (14) into equation (10) and taking Fourier transform give a new form of flame transfer function ²⁴

$${\widehat{q}}_k(s)=N_{\varphi ,k}\left[J{\widehat{f}}_k(s)-G{\widehat{v}}_{F^{-},k}(s)\right]e^{-s{\tau }_{\varphi }}$$

(15)

where

$$N_{\varphi ,k}=\{\begin{array}{l}1\text{for }{\beta }_k\le 1\hfill \\ 1-\frac{2{\psi }_k}{\pi }+\frac{2{(1-1/{\beta }_k{\hspace{0.17em}}^2)}^{1/2}}{\pi {\beta }_k}\text{for }{\beta }_k>1\hfill \end{array}$$

(16)

In equation (16), ${\beta }_k=\left|\left({\widehat{f}}_k/{\overline{f}}_u-{\widehat{v}}_{F^{-},k}/{\overline{v}}_u\right)n_{\varphi }/{\alpha }_{\varphi }\right|$ , ${\psi }_k=\arccos (1/{\beta }_k)$ , $J=n_{\varphi }\eta \Delta H_{\mathit{st}}$ , and $G=n_{\varphi }\eta \Delta H_{\mathit{st}}\overline{\varphi }FA{R}_{\mathit{st}}{\overline{\rho }}_u/(1+\overline{\varphi }FA{R}_{\mathit{st}})$ . The superscript “^” denotes the perturbation in the frequency domain.

To avoid unrealistic unstable modes at high frequencies, the NFTF, equation (15), may be filtered by a first-order low-pass filter^23,25

$${\widehat{q}}_k(s)=\frac{N_{\varphi ,k}}{{\tau }_{c}s+1}\left[J{\widehat{f}}_k(s)-G{\widehat{v}}_{F^{-},k}(s)\right]e^{-s{\tau }_{\varphi }}$$

(17)

where $1/{\tau }_c$ indicates the corner frequency (one chooses 50 Hz in the present work). It is generally assumed that the perturbation in fuel flow rate is close to zero ( $f_k^{\prime }\approx 0$ ) when the controller is not turned on and ${\Gamma }_k$ becomes

$${\Gamma }_k=\frac{{\widehat{q}}_k(s)}{{\widehat{v}}_{F^{-},k}(s)}=-\frac{G{N}_{\varphi ,k}}{{\tau }_{c}s+1}{e}^{-s{\tau }_{\varphi }}$$

(18)

And β _k becomes

$${\beta }_k=\left|\frac{n_{\varphi }}{{\alpha }_{\varphi }}\frac{{\widehat{v}}_{F^{-},k}}{{\overline{v}}_{u,k}}\right|$$

(19)

Analysis of the unstable azimuthal thermoacoustic mode

By inserting equation (18) into equation (9), a dispersion relation of azimuthal thermoacoustic modes is obtained. Note that ${\Gamma }_k$ depends on the fractional velocity perturbation $|{\widehat{v}}_{F^{-},k}/{\overline{v}}_u|$ when the controller is turned off (see equations (15) to (19)). Therefore, the dispersion equation (9) can be solved only if $|{\widehat{v}}_{F^{-},k}/{\overline{v}}_k|(k=1-M)$ are all given. Axisymmetric and non-axisymmetric cases are studied in the “Axisymmetric case of the unstable azimuthal thermoacoustic mode” and “Symmetry breaking of the unstable azimuthal thermoacoustic mode” sections, respectively.

Axisymmetric case of the unstable azimuthal thermoacoustic mode

One assumes that flames’ parameters and responses in each burner are the same. If the oscillation develops gradually from a weak disturbance, the fractional velocity perturbation at the location immediately upstream the flame inside each burner $|{\widehat{v}}_{F^{-},k}/{\overline{v}}_u|(k=1-M)$ grows simultaneously and equally, and the azimuthal thermoacoustic mode of the annular combustor is perfectly axisymmetry. One focuses on this case in this section.

One then solves equation (9) and gets the value of $s={\omega }_r+i{\omega }_i$ , where ${\omega }_i$ and ${\omega }_r$ indicate the eigenfrequency and growth rate of azimuthal mode, respectively. When ${\omega }_r>0$ , the corresponding mode is unstable. It is mentioned before that the bulk fluid flow and the viscothermal damping effects inside the burners are taken into consideration. One solves equation (9) under three conditions: (1) considering the viscothermal damping effects alone and ignoring the bulk fluid flow, (2) considering the bulk fluid flow alone and ignoring the viscothermal damping effects, and (3) considering both two factors. Figure 2 depicts the variation of growth rate of the first mode with ${\widehat{v}}_{F^{-},1}/{\overline{v}}_u$ for these three conditions. It can be seen from Figure 2 that the bulk fluid flow does have a great influence on the first mode, while the influence of viscothermal damping effects is relatively small but still exists. For naturally growing modes and their limit cycle states, only certain mode shapes over the circumference could be the solution of the system. However, the current work also deals with the situation that a random manner of acoustic disturbances which may also occur, e.g. a sudden fuel flow rate rise in a certain feeding line of the burner. Current research thus covers the mode shapes for a naturally growing disturbance (and their limit cycle states) and the random manners mentioned above. It should also be noted that random acoustic disturbances could lead to transient effect and mode dynamics. Current research, which only deals with the frequency domain model, is not sufficient to capture and predict these behaviors. A time domain framework is needed in order to study that. The present way to incorporate nonlinear flame models into this network model does not aim to provide predictive tools for thermoacoustic modes in annular combustors.

Figure 2.

Evolutions of the growth rate of the first mode with ${\widehat{v}}_k/{\overline{v}}_u$ for three conditions: result 1: one considers the viscothermal damping effect but ignores the bulk fluid flow; result 2: one considers the bulk fluid flow but ignores the viscothermal damping effect; result 3: both two factors are taken into account.

It is found that only the first mode is unstable, as shown in Figure 3. $V=\left({\widehat{v}}_{F^{-},1},{\widehat{v}}_{F^{-},2},{\widehat{v}}_{F^{-},3},{\widehat{v}}_{F^{-},4}\right)/{\overline{v}}_u$ denotes the fractional velocity perturbations immediately before four flames. Figure 4 shows the trajectory of the frequency and corresponding growth rate of the unstable mode. Part of the trajectory of this mode has a positive growth rate. The growth rate equals 33.5  ${\text{s}}^{\text{-1}}$ when ${\widehat{v}}_{{\text{F}}^{-},k}/{\overline{v}}_u\hspace{0.17em}=\hspace{0.17em}0$ (k could be any number from 1 to M since the fractional velocity perturbation at the flame location inside each burner is identical), meaning that the weak disturbances increase exponentially with an envelope of $e^{33.5t}$ . The growth rates remain unchanged when ${\widehat{v}}_{F^{-},k}/{\overline{v}}_u\le 0.1$ and becomes 0 when ${\widehat{v}}_{F^{-},k}/{\overline{v}}_{u,k}=0.19$ , hence this is the fractional velocity perturbation amplitude at which limit cycle oscillation occurs. Meanwhile, the eigenfrequency also changes slightly with increasing ${\widehat{v}}_{{\text{F}}^{-},k}/{\overline{v}}_u$ after entering nonlinear response stage. The parameters used for these frequency domain analyses are listed in Table 1.

Figure 3.

Frequencies and growth rates of the first three modes for three axisymmetric cases.

Figure 4.

Evolutions of the eigenfrequency and corresponding growth rate of the first mode with ${\widehat{v}}_k/{\overline{v}}_u$ for the perfect axisymmetric case.

Table 1.

Parameters used in the frequency domain analysis for both the axisymmetric and the non-axisymmetric cases.

$n_{\varphi }(-)$	${\alpha }_{\varphi }(-)$	${\tau }_{\varphi }(\text{s})$	${\tau }_c(\text{s})$	$\overline{\varphi }(-)$
−2.5	0.25	0.0025	0.003183	0.7
$\mathit{FA}{R}_{\mathit{st}}(-)$	$\Delta H_{\mathit{st}}(\text{J/kg})$	$\eta (-)$	$L(\text{m})$	$S({\text{m}}^{\text{2}})$
0.0671	44,316,000	0.8	0.54	0.00785
$L_k(\text{m})$	$S_k({\text{m}}^{\text{2}})$	$M(-)$	${\overline{\rho }}_u({\text{kg/m}}^{\text{3}})$	${\overline{c}}_u(\text{m/s})$
0.04	0.00028	4	10.93	519
${\overline{M}}_u(-)$	${\overline{\rho }}_b({\text{kg/m}}^{\text{3}})$	${\overline{c}}_b(\text{m/s})$	$\gamma (-)$	$\mathit{Rg}(\text{J/(kg}\text{K)})$
0.03	4.7	832	1.4	287.14

Symmetry breaking of the unstable azimuthal thermoacoustic mode

Thermoacoustic oscillations have a low amplitude initially during which the linearity of the flame response plays a leading role. Unstable mode will grow in perturbation amplitude until nonlinear effects become sufficiently important to achieve a limit cycle. At the nonlinear stage, if acoustic systems inside different burners are not identical due to some random factors (e.g. a large intensity disturbance occurs somewhere inside the annular combustor), the azimuthal mode will be strongly affected by this symmetry modification (for the saturation model, the difference of velocity disturbance immediately before each flame will not lead to the symmetry breaking when the heat release disturbances are in the linear response stage). The sum of the fractional velocity perturbations immediately before M flames is denoted as W which can be considered as an indicator of the oscillation development stage, expressed as

$$W={\displaystyle \sum _{k=1}^M\frac{{\widehat{v}}_{F^{-},k}}{{\overline{v}}_u}}$$

(20)

Taking $W=0.6$ as an example, the stability analysis is carried out for several perfect axisymmetry and non-axisymmetry cases. Figure 5 shows the calculation result. W is fixed in these cases, but the difference lies in the intensity of symmetry breaking. The growth rate of unstable mode increases with the intensity of symmetry breaking, implying that symmetry breaking would deteriorate the stability of annular combustors. The difference between modal growth rates of the axisymmetry case corresponding to $V=(0.15,0.15,0.15,0.15)$ and of the non-axisymmetry case corresponding to $V=(0.10,0.20,0.10,0.20)$ is about 15%, which is a relatively large value. It is worth noting that configuration ways of $V$ are theoretically infinite due to various random factors. A few of representative examples are selected herein to illustrate the effect of asymmetry on the unstable mode. Obviously, when symmetry breaking occurs, the state of the limit cycle oscillation is also diverse.

Figure 5.

Modal frequencies and corresponding growth rates of the unstable mode of axisymmetry and non-axisymmetric cases for a fixed value of $W$ .

Feedback control of azimuthal combustion instabilities

OLTF for control purposes

The OLTF for control purposes is from the voltage signal $\widehat{U}$ produced by the valve-driving machine to the signal ${\widehat{p}}_r$ measured by the pressure transducer, as shown in Figure 6. One assumes that the valve on the fuel feed line can modulate the fuel flow rate of all M burners simultaneously and equally. A microphone which has the same azimuthal angle with one of the burners is mounted on the plenum as the pressure sensor, as shown in Figure 1(a). A SISO controller is achieved in this work. The microphone electrical output signal is assumed to be the measured pressure—no sensor transfer function is used. The applied voltage $\widehat{U}$ is associated with the resulting fractional change in fuel flow rate within the fuel feed line ( ${\widehat{f}}_{\mathit{total}}(s)/{\overline{f}}_{\mathit{total}}$ ) by the valve transfer function. The form of this transfer function is the same as that used in Morgans and Dowling ¹¹

$$L_t(s)=\frac{{\widehat{f}}_{\mathit{total}}(s)/{\overline{f}}_{\mathit{total}}}{\widehat{U}}={\kappa }_t{e}^{-s{\gamma }_t}$$

(21)

where ${\kappa }_t=0.1{\text{V}}^{-1}$ , ${\gamma }_t=2\times {10}^{-4}\text{s}$ , and $f_{\mathit{total}}(t)$ is the total fuel flow rate. This form captures the feature of the fuel valve that its phase lag increases with increasing the frequency. It is recommended to refer to Morgans and Dowling ¹¹ for more details on this valve transfer function.

Figure 6.

Block diagram of the feedback control system.

The NFTFs have the form expressed as equation (17) when the controller is turned on. Substituting equation (17) into equation (13) in Bauerheim et al. ¹³ and after a series of derivations, the following expression can be obtained

$${\widehat{p}}_{u,k}(x=0,s)=E(s){\widehat{f}}_k(s)$$

(22)

where $E(s)$ is the transfer function from ${\widehat{f}}_k$ to ${\widehat{p}}_r$ (note that ${\widehat{p}}_r(s)={\widehat{p}}_{u,k}(x=0,s)$ ) and ${\widehat{f}}_k(s)$ is the fluctuating component of the fuel flow rate per unit area for the burner k. Obviously, the following relationship exists since each burner has the same geometry and fuel flow rate

$$\frac{{\widehat{f}}_{\mathit{total}}(s)}{{\widehat{f}}_k(s)}=M\cdot S_k$$

(23)

Combing equations (21) to (23) then gives the modal OLTF $O(s)$ needed for the controller design

$$O(s)=\frac{{\widehat{p}}_r(s)}{\widehat{U}(s)}=L_t(s)E(s)\frac{{\overline{f}}_{\mathit{total}}}{M\cdot S_k}$$

(24)

Due to the flame nonlinear response and the symmetry breaking of the azimuthal thermoacoustic modes, $O(s)$ varies depending on fractional velocity perturbations immediately before the M flames $V=\left({\widehat{v}}_{F^{-},1},{\widehat{v}}_{F^{-},2},{\widehat{v}}_{F^{-},3},{\widehat{v}}_{F^{-},4}\right)/{\overline{v}}_u$ . Figure 7 shows the bode plot of $O(s)$ for a series of plants with different $V$ .

Figure 7.

Bode plot of the $O(s)$ for five different V. (a) transfer function gain. (b) transfer function phase.

By assuming each mode to be described as a second-order transfer function, the stability of the system can be assessed by checking the phase change across the peaks. ³ A 180° phase increasing across a peak indicates an unstable conjugate pair of poles while a 180° phase decrease indicates a stable conjugate pair of poles. ¹¹ From Figure 7, the plant of $V=(0.19,0.19,0.19,0.19)$ corresponds to a stable OLTF, while the rest of five plants correspond to unstable OLTFs. It can be seen that the controller designed to stabilize the linear OLTF cannot guarantee to stabilize all possible OLTFs since due to the nonlinear feature. In order to address this problem, we need a controller which could both reduce the disturbance amplitude when the system is in limit cycle (including the axisymmetric and the non-axisymmetric cases) and obtain closed-loop stability for all possible plants (it is noted that when the system is in limit cycle stage, the growth rate equals zero but the system is not stable). To fulfill these requirements, a controller with enough robustness is required.

$H_{\infty }$ loop-shaping and the $\nu \text{-}\mathit{gap}$ metric

$H_{\infty }$ loop-shaping

The $H_{\infty }$ loop-shaping method is used in this work to design robust controllers. This method which combines the loop-shaping design and the robust stabilization was proposed by McFarlane and Glover ²⁶ and is described as follows ²⁷ :

The OLTF $O(s)$ is shaped with frequency-dependent pre-or post-compensator $W_1$ and $W_2$ so that they have a high gain where disturbance rejection is important, and a low gain where robustness is required. The shaped plant becomes $O_W=W_{1}O{W}_2$ . For SISO systems, only one of the compensators is needed and the other one thus equals unity. Then the shaped plant becomes $O_W=W_{1}O$ (or $O{W}_2$ ).

A stabilizing controller $K_{\infty }$ is then synthesized via the Matlab package ncfsyn, which minimizes the $H_{\infty }$ norm of the closed-loop transfer function (see equation (25)). The norm $b_{\infty }$ provides a measure of the system’s robustness ( $b_{\infty }$ is also called stability margin) ⁷

$$b_{\infty }\mathrm{=}{\Vert \left[\begin{array}{l}I\\ O_W\end{array}\right]{(I-K_{\infty }{O}_W)}^{-1}[\begin{array}{cc}I& K_{\infty }\end{array}]\Vert }_{\infty }^{-1}$$

(25)

If $b_{\infty }$ is suitably small, then a robust controller that meets the requirements is obtained. The phase of the shaped transfer function then will be altered so as to achieve stability but the gain will not be significantly affected. ⁷ If the norm is not suitably small, then another iteration on the choice of the compensators $W_1$ or $W_2$ must be performed.

$\nu \text{-}\mathit{gap}$ metric

The difference between two open-loop plants can be measured by the $\nu \text{-}\mathit{gap}$ metric. ¹⁴ For SISO systems, the $\nu \text{-}\mathit{gap}$ between plants $O_1$ and $O_2$ , denoted by ${\delta }_{\nu }$ , is given by equation (26) supposing that the winding number condition is satisfied ⁷

$${\delta }_{\nu }(O_1,O_2)=\underset{{\omega }_i}{\sup }\psi (O_1,O_2)=\underset{{\omega }_i}{\sup }\frac{|O_1({\omega }_{i}i)-O_2({\omega }_{i}i)|}{{(1+|O_1({\omega }_{i}i){|}^2)}^{1/2}{(1+|O_2({\omega }_{i}i){|}^2)}^{1/2}}$$

(26)

The $\nu \text{-}\mathit{gap}$ metric fits highly naturally into the $H_{\infty }$ loop-shaping synthesis framework. ⁷ If the $\nu \text{-}\mathit{gap}$ between two plants (e.g. $O_1\&O_2$ ) meet the condition $b_{\infty }>{\delta }_{\nu }$ , then the controller which guarantees the robust stability for one plant ( $O_1$ ) must robustly stabilize the other one ( $O_2$ ).

Controller design

Due to the flame nonlinearities introduced by the NFTF, the OLTF varies with the fractional velocity perturbations immediately before the M flames $V$ . A nominal linear plant need to be sought as the input argument of the Matlab command ncfsyn which performs the $H_{\infty }$ loop-shaping synthesis. The difference between the nominal plant and all possible OLTF $O_i(V,s)$ can be measured using the $\nu \text{-}\mathit{gap}$ metric. The design procedure is interpreted by the form of a block diagram, as shown in Figure 8.

Figure 8.

Block diagram of the feedback controller design process. NFTF: nonlinear flame transfer function; OLTF: open-loop transfer function; RTF: rational transfer function.

From Figure 4, flame nonlinearity acts when $V>(0.1,0.1,0.1,0.1)$ , and the limit cycle oscillation occurs when $V=(0.19,0.19,0.19,0.19)$ for the perfectly axisymmetric case. Considering the monotonicity of the saturation flame model ²⁸ (see Figure 3 in Stow and Dowling ²⁸ for more details), one chooses $O(V=(0.15,0.15,0.15,0.15),s)$ as the nominal OLTF $O_N(s)$ . The system is infinite dimensional because of the presence of the time delays in the OLTF (e.g. the time delays ${\tau }_c\&{\tau }_{\varphi }$ in the NFTF and the measurement time delay, see equations (18) and (21)). The Matlab fitting procedure fitfrd is adopted in this work to acquire a rational transfer function (RTF). We fit the OLTF over [0, 1000] Hz since the NFTF is valid in low-frequency range. ²⁵ The resulting RTF $O_{\mathit{NR}}(s)$ has an order of 16. Figure 9 shows the comparison between the original plant and the fitted plant. The error between $O_N(s)$ and $O_{\mathit{NR}}(s)$ can be measured by the $\nu \text{-}\mathit{gap}$ metric, as shown in Figure 10. The value of the $\nu \text{-}\mathit{gap}$ over [0, 1000] Hz is always smaller than $1\times {10}^{-5}$ which is a sufficiently small value, implying that the fitted plant matches perfectly the original plant.

Figure 9.

Comparison between the original OLTF and fitted OLTF (high-order fitting). OLTF: open-loop transfer function.

Figure 10.

The $\nu \text{-}\mathit{gap}$ metric between the original OLTF and the high-order RTF.

The fractional velocity perturbations $|{\widehat{v}}_{F^{-},k}/{\overline{v}}_u|(k=1-M)$ are randomly selected within the range [0, 0.23] and 10,000 OLTFs $O_i(V,s)(i=1-10,000)$ are obtained (one assumes that all possible plants have been covered since the data size is large enough). The $\nu \text{-gaps}$ among $O_i(V,s)(i=1-10000)$ and $O_{\mathit{NR}}(s)$ are calculated by equation (26) and among which the maximum value is $\max ({\delta }_{\nu })=0.31$ . Figure 11(a) shows the evolution of $\psi (O_{\mathit{NR}}(s),O_i(s))$ with frequency for several typical cases of $V$ .

Figure 11.

Plot of $\psi (O_{\mathit{NR}}(s),O_i(s))$ with frequency for different $V$ . (a) High-order fitting and (b) low-order fitting. The black dot-dash line represents the robust stability margin ϵ of the controller.

The selection of the compensator governs the $H_{\infty }$ loop-shaping design procedure as it affects the entire process of the loop-shaping design. A pre-compensator with the form

$$W_1(s)=2.5\times {10}^{-2}\frac{{3581}^2}{s^2+7.162s+{3581}^2}$$

(27)

is chosen. It consists of a proportional element and a second-order low-pass filter. It has a gain fall-off at high frequencies for robustness.

After shaping the singular value of the fitted nominal plant $O_{\mathit{NR}}(s)$ by the selected weight $W_1$ , the $H_{\infty }$ loop-shaping synthesis is then performed via the Matlab command ncfsyn. An 18 order controller $K_{\infty }$ is generated with a stability margin value of 0.37. The stability margin $b_{\infty }$ is larger than $\max ({\delta }_{\nu })$ calculated above, as shown in Figure 11(a). This implies that the controller could stabilize all the possible plants.

The closed-loop transfer function has the form as equation (28). If any isolated closed-loop pole locates in the right half plane, the system is unstable. On the contrary, the system is stable when there is no isolated closed-loop pole located in the right half plane. ¹¹ Then same open-loop plants $O_i(s)$ as shown Figure 11(a) are used here to verify the effectiveness of the controller. As shown in Figure 12, all these plants are stabilized

$$C_i(V,s)=\frac{O_i(V,s)}{1-K_{\infty }(s)O_i(V,s)}$$

(28)

Figure 12.

Pole-zero plot of the closed-loop transfer function for five typical cases of V. (a) High-order controller and (b) low-order controller. The circle represents the zeros and the cross represents the poles.

When the heat release rate oscillations grow sufficiently large, limit cycle oscillations are established and the corresponding OLTF seems to be stable even without feedback control. Under this condition, sensitivity transfer functions could be used to measure the performance of the controller ⁷

$$S(s)=\frac{1}{1+K_{\infty }(s)O_i(s)}$$

(29)

When the sensitivity function is smaller than 1, the feedback controller can attenuate the disturbances. The limit cycle oscillation states are not unique due to the consideration of symmetry breaking. Figure 13(a) shows plots of $\left|S\right|$ versus frequency for several $V$ corresponding to different limit cycle oscillation states. The frequency range calculated is 130–150 Hz. (The frequencies of the unstable modes are always very close to 138 Hz, as shown in Figures 4 and 5.) The fact that $\left|S\right|<1$ at these frequencies implies that the feedback controller reduces the oscillation amplitude from within the limit cycles, the plant appears unstable, and the controller then stabilizes it.

Figure 13.

Plots of $\left|S\right|$ with frequency for several limit cycle oscillation states. (a) High-order controller and (b) low-order controller.

Low-order $H_{\infty }$ controller design method

A $H_{\infty }$ loop-shaping controller has been designed in the “Controller design” section to effectively suppress the azimuthal combustion instabilities from within the limit cycles. The order of the controller is 18. This part introduces a method to design a low-order $H_{\infty }$ loop-shaping controller proposed by Li and Morgans ⁷ and applies this method to the control of the azimuthal combustion instabilities in the annular combustors. Different from the traditional model reduction algorithm, ²⁹ this method reduces the order of the open-loop plant to get the low-order controller directly since the order of the designed controller is equal to the sum of the order of the shaped nominal plant and the weights.

A four order RTF is used to fit the original plant. As shown in Figure 14, the low-order RTF is well fitted at the low-frequency range and has larger deviations at high frequencies than the high-order fitted RTF. Same as high-order controller design procedure, 10,000 OLTFs $O_i(V,s)(i=1-10000)$ are obtained. The $\nu \text{-gaps}$ ${\delta }_{\nu }$ between $O_i(V,s)(i=1-10,000)$ and $O_{\mathit{NR}}(s)$ are calculated by equation (26) and among which the maximum value is $\max ({\delta }_{\nu })=0.303.$ Figure 11(b) shows the evolution of $\psi (O_{\mathit{NR}}(s),O_i(s))$ with frequency for several typical cases of $V$ .

Figure 14.

Comparison between the original OLTF and fitted OLTF (low-order fitting). OLTF: open-loop transfer function.

The same weights as equation (27) are used to reshape the low-order fitted nominal plant. The final controller has an order of 6 in the denominator and 4 in the numerator

$$K_{\infty }=\frac{9.643\times {10}^5(s^2+0.0068s+12.82)(s+2.0596)(s+0.0512)}{(s^2-0.0648s+12.9)(s^2+0.076s+12.74)(s+1.9989)(s+0.4639)}$$

(30)

The robust stability margin $b_{\infty }$ is 0.352, larger than $\max ({\delta }_{\nu })$ , guaranteeing the system’s stability. Figure 12(b) is the pole-zero plot of the closed-loop transfer function for several typical cases of $V$ . It shows that all these plants obtain closed-loop stability. Figure 13(b) shows the plot of $\left|S\right|$ versus frequency for several $V$ corresponding to different limit cycle oscillation states. The gain of sensitivity function is always smaller than 1. The result is very close to that of the higher order controller. From the above fact, it is inferred that the low-order controller could successfully suppress the oscillation from within the limit cycles.

Conclusions

In this work, a theoretical annular combustion chamber model is improved by taking into account the effects of the flame nonlinearity, the bulk fluid flow, and the viscothermal wall damping mechanisms inside the burners which are ignored in the original model for the sake of simplicity. This improved model is used to predict the azimuthal combustion instability and develop the corresponding model-based active control strategies for the annular combustors. A pressure transducer mounted on the plenum is used as the sensor and a valve on the fuel feed line as the actuator.

The nonlinear response of the flame to the flow disturbances is modeled using a simple gain/time delay model with saturation bounds. The NFTF may vary depending on the flow disturbance levels, and the NFTF within each burner may therefore be different since the flow state in each burner can be different. It is obvious that the OLTF of the controlled system is not unique, and one of them is selected and then fitted to obtain a nominal OLTF (RTF) for controller design.

The distances between all the possible OLTFs and the nominal OLTF are measured using the $\nu \text{-}\mathit{gap}$ metric, among which the maximum value sets the minimum robustness stability margin required of the controller. The $H_{\infty }$ loop-shaping synthesis is then performed and a SISO controller meeting the robustness requirement is obtained. A standard high-order controller and a low-order controller have been designed. Both controllers are guaranteed to stabilize the annular combustor for all the possible OLTFs.