4.2.1 Mode

To analyse the simulation results conveniently, we set three planes in the domain perpendicular to the chamber axis shown in Fig. 15, with each plane having four probes. Probes of the 1st to 4th were located counter-clockwise on plane 1, probes of the 5th to 8th were located counter-clockwise on plane 2 and probes of the 9th to 12th were located counter-clockwise on plane 3.

Figure 15. Schematic of the analysis planes and the probes.

Section 4.1 has demonstrated the mean flow field. Here, mode shapes are investigated. After the limit cycle has formed, the waveforms at probe 1, probe 2, probe 3 and probe 4 were similar, and they exhibited a phase difference for the adjacent probes. To show details, plots of unfiltered pressure time histories taken at probes 1, 2, 3 and 4 are shown in Fig. 16, which gives the waveforms of limit cycle. The travelling compression wave can be seen by examining the peaks in the pressure data. Combining the positions of the probes, it is expected that the spinning first tangential (1T) mode applies to cases (a), (b) and (c). Note that the spinning direction may be arbitrary. For case (d), as shown in Fig. 16(d), the signals taken at probe 2 and probe 4 are out of phase with signals taken at probe 1, the reference signals. Likewise, the signals taken at probe 3 are in phase with signals of probe 1. This indicates that the dominant instability has the features of a standing second tangential (2T) mode for case (d). These features will be highlighted later.

Figure 16. Pressure time histories through probe 1 to probe 4.

Time snapshots of the entire combustion flow field at specific instants of time are useful for deducing the global character of the entire flow field. When these time snapshots are spaced sufficiently closely together in time, they can be converted to flow movies, which are helpful for understanding dynamics. Figure 17 shows the temporal behaviour of the pressure field over an entire cycle for different cases. As discussed earlier, 1T oscillations are observed for cases (a), (b) and (c), and their modes are found to be spinning, while a standing 2T wave is clearly observed for case (d). It is worth noting that azimuthal modes can appear as standing or spinning modes, and both are observed in gas turbines. Bifurcations between standing and turning modes were explained previously. Schuermans et al.^{(Reference Schuermans, Paschereit and Monkiewitz36–Reference Wolf, Staffelbach, Gicquel, MÜller and Poinsot38)} attribute these bifurcations to non-linear effects. They propose a non-linear theoretical approach showing that standing wave modes can be found at low oscillation amplitudes, but only one spinning mode is found for large-amplitude limit cycles. Similarly, we simulated many different cases having unique heat release zones but different maximum $Q^{\prime}$ (Equation (8)). Unfortunately, these bifurcations cannot be observed in our calculations. Some researchers explained these bifurcations using linear approaches: spinning modes would appear only in perfectly axisymmetric configurations, while any symmetry modification would lead to standing modes^{(Reference Wolf, Staffelbach, Gicquel, MÜller and Poinsot38)}. Results from Ref. (Reference Evesque, Polifke and Pankiewitz39) suggested that the initial flow conditions might trigger either standing or spinning modes. Some researchers observed that turbulence can cause random mode switching between standing and spinning structures^{(Reference Schuermans, Paschereit and Monkiewitz36,Reference Noiray, Bothien and Schuermans40)} . According to present calculations, it was found that modulation of heat release length might cause mode switching between standing and spinning structures. Increment of heat release length tended to excite higher-order standing modes.

Figure 17. Temporal behaviour of the pressure field over an entire cycle on the slice ( $x = - 9.375$ ). Arrowhead denotes the spinning direction of the pressure wave. Conditions correspond to Table 2.

Transverse oscillation modes have been discussed. In what follows, attention is focused on the longitudinal modes. A clear picture of the underlying longitudinal modes emerges in the present simulation, and the resulting wall pressure is plotted at different instants in time over an entire cycle. Such longitudinal mode shapes are shown in Fig. 18. Superimposed on Fig. 18 are the labels of the heat release zone. Examining this figure, it is apparent that the mode shape in the heat release zone appears to be approaching the half-wave mode for cases (a), (b) and (c). As the heat release zone is extended to $x = - 0.625$ from the faceplate (case (d)), the first longitudinal mode becomes the dominant mode. For the rest of the chamber region, the case (a) mode is the third longitudinal mode and the case (b) mode is the second longitudinal mode, whereas the half-wave mode is clearly observed for case (c) and case (d). All cases considered in this context showed a complex wave shape at any given instant in time, and the pressure reached a peak at different nodal locations at different instances in time. As Smith^{(Reference Smith, Ellis, Xia, Sankaran, Anderson and Merkle29)} explained, it is the combined effects of the forward-running $u + c$ and the backward-running $u - c$ waves that cause this inherent phase difference in the wave shape. In all, when the heat release zone locates near the faceplate where it is the pressure antinode of the longitudinal mode, increment of heat release length tends to excite lower-order longitudinal mode.

Figure 18. Computed mode shapes for the longitudinal mode. Conditions correspond to Table 2.

Overall, combined longitudinal-transverse acoustic modes occurred in the present air heater, illustrated in Fig. 19.

Figure 19. Pressure mode shapes for different heat release zones.

When the system reached the limit cycle, the combustion response function was turned off. Figure 20 shows the time histories of pressure oscillations for case (a) after its combustion response function was turned off. A full view of the decay process for pressure fluctuation is shown in the inset of Fig. 20. The fluctuations are observed to be damped, and the flow field gradually becomes stable. In the absence of combustion response function, it was found that the numerical and real viscosity was sufficient to damp out any oscillations within a short time. Accordingly, the unsteady results presented in this paper were predominately acoustic in nature.

Figure 20. Time histories of pressure oscillations for case (a). The combustion response function is turned off at 0.577s.

4.2.2 Frequency

We next turn our attention to the frequency response predictions of the present air heater. Figure 21 depicts a waterfall plot for the pressure signals at probe 1 after the limit cycle has developed, which shows the transition of the frequencies and amplitudes as time proceeds. The sampling duration is 80ms, and the frequency resolution is 12.5Hz. As anticipated, in Fig. 21 the first three harmonics can be clearly observed, and these frequencies do not change considerably during the period of the limit cycle. The higher peaks represent harmonics of the primary mode occurring at integer multiples of its frequency, but with a steadily decreasing amplitude.

Table 5 Frequency comparisons using the classical acoustic model and the present computational model

Figure 21. Waterfall plot for the pressure signals at probe 1 for $\sim$ 80ms.

The resonant mode frequencies of the air heater may be estimated using Equation (13) based on the following assumptions: the chamber of the air heater was viewed as a closed cylindrical duct with no nozzles and injectors; inviscid uniformly distributed ideal gas was considered; the condition in the combustor was assumed to satisfy linear acoustics, namely small-amplitude oscillations; and there was low-speed mean flow. Table 5 summarises the frequency using the classical acoustic model and the present computational model. Note that the chamber consists of two cylindrical sections with different diameters and lengths. In addition, there exists a notable temperature gradient along the flow direction in the first domain of the combustor. Thus, only the conditions in the last domain of the combustor satisfy the above assumptions. The results obtained from the present model were extracted from Fig. 21. It is evident that the computed frequency agrees fairly well with the theoretical value.

(13) \begin{equation}{f_{mnq}} = \frac{c}{2}\sqrt {{{\left( {\frac{{{\alpha _{mn}}}}{{{R_c}}}} \right)}^2} + {{\left( {\frac{q}{{{L_c}}}} \right)}^2}}\end{equation}

where ${\alpha _{mn}}$ is ${n_{th}}$ root of the ${J_n}$ Bessel function of the first kind, c is the speed of sound in the medium, ${L_c}$ is the length of cylinder, ${R_c}$ is the radius of the cylinder and m, n and q are the wavenumbers for the tangential, radial and longitudinal modes, respectively.

It is interesting to compare the oscillation characters found in experiment with those obtained from computation and theoretical analysis. As already noted, strong peaks around 427.2Hz, 1,112Hz, 1,322Hz and 1,570Hz are observed on hot firing test (Fig. 3). The theoretical frequency of the first longitudinal mode is about 448Hz, estimated by Equation (13). Interestingly, the first experimental dominant frequency (427.2Hz) matches this value. The second dominant frequency (1,112Hz) shows good match with the frequency found in case (c). The third dominant frequency (1,322Hz) matches well with the frequency found in case (b). Also, the match between the last dominant frequency (1,570Hz) and that found in case (a) is good. Note that, according to our calculations, every case has its own oscillation modes and the corresponding frequencies, and these modes and frequencies, were insensitive to the changes of interaction index between heat release and acoustic modes, and to the maximum $Q^{\prime}$ (Equations (7) and (8)), but sensitive to the heat release length. It was logical to assume that the flame might exhibit a significant unsteady motion; thus, the heat release length varied during the experiment, which led to different oscillation modes and frequencies. A waterfall plot of chamber pressure may support this assumption (Fig. 3), the peak frequency exhibiting maximal amplitude shifts among the dominant frequencies as the experiment time progressed, indicating that the heat release length is changing. However, this assumption needs more tests. Overall, the present model successfully predicted the dominant frequencies observed on hot firing test, except the first pure longitudinal mode.

4.3 Theoretical analysis of the modal behaviour in air heater combustor, i.e. standing or spinning waves

Considering the geometry of the present air heater combustor, the wave equations can be written in polar coordinates:

(14) \begin{equation}\frac{{{\partial ^2}p^{\prime}}}{{\partial {\theta ^2}}} - \alpha \frac{{\partial p^{\prime}}}{{\partial t}} - \frac{{{\partial ^2}p^{\prime}}}{{\partial {t^2}}} = - \frac{{\partial Q^{\prime}}}{{\partial t}}\end{equation}

The Galerkin method is used to solve Equation (14), assuming that the pressure field can be expressed as the sum of a series of orthogonal basis functions $p^{\prime}\left( {{\rm{\theta }},{\rm{t}}} \right) = \mathop \sum \nolimits_1^\infty {\eta _n}\left( t \right){\psi _n}\left( \theta \right)$ . When considering the situation where the thermoacoustic coupling in the air heater results from the combination of the eigenmodes sharing the nth eigenvalue, the pressure field can be approximated by Equation (15).

(15) \begin{equation}p^{\prime}(\theta ,t) = {\eta _1}(t) \cdot \cos (n \cdot \theta ) + {\eta _2}(t) \cdot \sin (n \cdot \theta )\end{equation}

It is assumed that the modal amplitudes can be expressed as Equation (16), in which A and B vary slowly in time.

(16) \begin{equation}\left\{ {\begin{array}{*{20}{c}}{{\eta _1} = A(t)\sin (nt)}\\\\[-8pt] {{\eta _2} = B(t)\cos (nt)}\end{array}} \right.\end{equation}

Thus, the expression for pressure field $p^{\prime}$ can be decomposed as Equation (17).

(17) \begin{equation}\begin{array}{l}p(\theta ,t) = {\eta _1}(t) \cdot \cos (n \cdot \theta ) + {\eta _2}(t) \cdot \sin (n \cdot \theta )\\{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = A\left( t \right)\sin \left( {nt} \right) \cdot \cos (n \cdot \theta ) + B(t)\cos \left( {nt} \right) \cdot \sin (n \cdot \theta )\\{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \underbrace {B\left( t \right)\sin \left( {nt + n\theta } \right)}_{Spinning{\kern 1pt} {\kern 1pt} {\kern 1pt} mode} + \underbrace {\left[ {A\left( t \right) - B(t)} \right]\sin \left( {nt} \right)\cos \left( {n\theta } \right)}_{Standing{\kern 1pt} {\kern 1pt} mode}\end{array}\end{equation}

From Equation (17), it can be concluded that, when the modal amplitudes A and B are different, the result is a combination of standing and spinning waves. When the modal amplitudes A and B are identical, the situation corresponds to a spinning wave.

The heat release response function is written as Equation (18).

(18) \begin{equation}Q^{\prime} = F\left[ {p^{\prime}\left( {t - \tau } \right)} \right] = \beta \cdot p^{\prime}\left( {t - \tau } \right) - k \cdot {p{\prime}^3}\left( {t - \tau } \right),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \beta = {\left. {\frac{{dF}}{{dp^{\prime}}}} \right|_{p^{\prime}=0}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} k{\kern 1pt} {\kern 1pt} = - \frac{1}{6}{\left. {\frac{{{d^3}F}}{{d{{p{\prime}}^3}}}} \right|_{p^{\prime}=0}}\end{equation}

Considering a general case in which the heat release response is azimuthal non-uniform, the heat release distribution parameter β can be expressed as Equation (19).

(19) \begin{equation}\beta = {\beta _0} \cdot \left\{ {1 + \sum\limits_{m=1}^M {\left[ {{C_m}\cos \left( {m\theta } \right) + {D_m}\sin \left( {m\theta } \right)} \right]} } \right\}\end{equation}

Differential equation for ${\eta _1}$ can be obtained by substituting Equations (15) and (18) into Equation (14), multiplying both sides with ${\rm{cos}}\left( {n\theta } \right)$ and integrating over $\theta$ ranging from 0 to 2Π. A similar expression for ${\eta _2}$ is obtained by multiplying both sides with sin(n $\theta$ ) before integration. These two differential equations are in the following form. For details of the mathematical development, the reader should consult Ref. (Reference Noiray, Bothien and Schuermans40).

(20) \begin{equation}\begin{array}{l} {\ddot{\eta }_{1} \left(t\right)+\alpha \dot{\eta }_{1} \left(t\right)+n^{2} \eta _{1} \left(t\right)=\beta _{0} \left(1+\frac{C_{2n} }{2} \right)\cdot \dot{\eta }_{1} \left(t-\tau \right)} \\ \\[-7pt]{-3k{\left\{\dot{\eta }_{1} \left(t\right)\left[3\eta _{1}^{2} \left(t\right)+\eta _{2}^{2} \left(t\right)\right]+2\dot{\eta }_{2} \left(t\right)\eta _{1} \left(t\right)\eta _{2} \left(t\right)\right\}\mathord{\left/ {\vphantom {\left\{\dot{\eta }_{1} \left(t\right)\left[3\eta _{1}^{2} \left(t\right)+\eta _{2}^{2} \left(t\right)\right]+2\dot{\eta }_{2} \left(t\right)\eta _{1} \left(t\right)\eta _{2} \left(t\right)\right\} 4}} \right. \kern-\nulldelimiterspace} 4} } \end{array}\end{equation}

(21) \begin{equation}\begin{array}{l}{\ddot{\eta }_{2} \left(t\right)+\alpha \dot{\eta }_{2} \left(t\right)+n^{2} \eta _{2} \left(t\right)=\beta _{0} \cdot \left(1-\frac{C_{2n} }{2} \right)\cdot \dot{\eta }_{2} \left(t-\tau \right)} \\ \\[-7pt] {-3k{\left\{\dot{\eta }_{2} \left(t\right)\left[3\eta _{2}^{2} \left(t\right)+\eta _{1}^{2} \left(t\right)\right]+2\dot{\eta }_{1} \left(t\right)\eta _{1} \left(t\right)\eta _{2} \left(t\right)\right\}\mathord{\left/ {\vphantom {\left\{\dot{\eta }_{2} \left(t\right)\left[3\eta _{2}^{2} \left(t\right)+\eta _{1}^{2} \left(t\right)\right]+2\dot{\eta }_{1} \left(t\right)\eta _{1} \left(t\right)\eta _{2} \left(t\right)\right\} 4}} \right. \kern-\nulldelimiterspace} 4} } \end{array}\end{equation}

The MATLAB solver was used to solve the delay differential equations (Equations (20, 21)). Figure 22 shows the calculation results of the limit cycle amplitudes ${\eta _1}$ and ${\eta _2}$ for different ${C_{2n}}$ and τ. It can be seen that when ${C_{2n}}=0$ , i,e. the heat release distribution is azimuthally uniform, there will be a spinning mode, since the modal amplitude A is equal to B. As the index ${C_{2n}}$ of the heat release azimuthal distribution increases, the modal amplitudes A and B become different; the situation is a combination of standing and spinning waves. When ${C_{2n}}$ becomes larger than its critical value, the standing wave is established. According to this analysis, when the length of the heat release zone is increased, the azimuthally inhomogeneous degree of heat release is increased. Thus, the spinning mode is established for the shorter heat release zone case, while the standing mode is established for the longer heat release case.

Figure 22. The limit cycle amplitudes ${\eta _1}$ and ${\eta _2}$ for different ${C_{2n}}$ and τ.