2.1 Acoustic model

The acoustic network we consider is a duct of length $L=1$  m with an inlet of area $A_{1}=1.96\times 10^{-3}~\text{m}^{2}$ . We prescribe the inlet mean flow $\overline{u}_{1}=0.4~\text{m}~\text{s}^{-1}$ , mean pressure $\overline{p}_{1}=1.01\times 10^{5}~\text{Pa}$ and temperature $T_{1}=300~\text{K}$ . The gas is considered to be ideal, obeying the equation of state $p=\unicode[STIX]{x1D70C}RT$ , where $\unicode[STIX]{x1D70C}$ is the air density and $R=287~\text{J}~\text{kg}^{-1}~\text{K}^{-1}$ the air gas constant. Across the heater, located at a distance $x_{h}$ downstream, a temperature jump $\unicode[STIX]{x1D6E5}\equiv (T_{2}/T_{1})^{1/2}=c_{2}/c_{1}=1.4$ is determined by the heater mean heat release and an area change $\unicode[STIX]{x1D6E9}\equiv A_{2}/A_{1}=1.1$ . Here we have defined the speed of sound $c\equiv \sqrt{\unicode[STIX]{x1D6FE}RT}$ , where $\unicode[STIX]{x1D6FE}=1.4$ is the specific heat ratio of air. Subscripts $_{1}$ and $_{2}$ denote variables upstream and downstream the heater, respectively. We decompose the acoustic velocity, pressure and density fluctuations into downstream ( $f$ ) and upstream ( $g$ ) travelling acoustic waves and an entropy wave. When entropy waves are accelerated, for example at a choked exit, they in turn generate indirect acoustic waves. In our system, however, this phenomenon is not modelled, and entropy waves are simply convected by the uniform mean flow out of the domain.

Mass, momentum and energy fluxes are conserved through the heater via the Rankine–Hugoniot jump conditions (Dowling Reference Dowling1995; Stow & Dowling Reference Stow and Dowling2001). The reflection coefficients provide the relations $f_{1}=R_{1}g_{1}\text{e}^{-s\unicode[STIX]{x1D70F}_{1}}$ and $g_{2}=R_{2}f_{2}\text{e}^{-s\unicode[STIX]{x1D70F}_{2}}$ , where $s=\unicode[STIX]{x1D70E}+\text{i}\unicode[STIX]{x1D714}$ is the Laplace variable, $\unicode[STIX]{x1D70F}_{1}\equiv 2x_{h}c_{1}/(c_{1}^{2}-\overline{u}_{1}^{2})$ and $\unicode[STIX]{x1D70F}_{2}\equiv 2(L-x_{h})c_{2}/(c_{2}^{2}-\overline{u}_{2}^{2})$ . The inlet and outlet reflection coefficients are chosen to be $R_{1}=R_{2}=-0.9$ . This value is chosen to be larger than $-1$ to account for the fact that the presence of a mean flow changes the expression for the conservation of acoustic energy at the boundaries (Polifke Reference Polifke2011). Because we are not interested in calculating explicitly the effect of entropy waves, we substitute the mass equation into the momentum and energy equations (Dowling Reference Dowling1997). The remaining equations describe the evolution of acoustic waves only, without neglecting entropy waves.

By following a procedure analogous to that described in Dowling (Reference Dowling1995), Orchini et al. (Reference Orchini, Illingworth and Juniper2015), we calculate the linear acoustic response to heat release fluctuations $q^{\prime }$ , which is given by the equations:

which represent momentum and energy conservation across the heater element.

The coefficients in the matrix $\unicode[STIX]{x1D648}$ are reported in appendix A. By setting the determinant of $\unicode[STIX]{x1D648}$ equal to zero (nonlinear eigenvalue problem) we find the acoustic eigenfrequencies. Solving for the wave amplitudes as a function of $q^{\prime }$ provides the acoustic transfer functions of pressure, velocity and density fluctuations with respect to heat release oscillations at any point in the duct. We are interested in the velocity response $u^{\prime }$ just upstream the flame for the coupling with the heater, which we measure as the frequency response $u^{\prime }/q^{\prime }\equiv H(\text{i}\unicode[STIX]{x1D714})$ , where $\unicode[STIX]{x1D714}$ is the oscillation angular frequency. Further details on the acoustic modelling are provided in Orchini et al. (Reference Orchini, Illingworth and Juniper2015).

2.2 Heat release model

King’s law (King Reference King1914) expresses the heat transferred from a hot-wire to the flow under steady flow conditions. With a quasi-steady argument, (i.e. assuming that the instantaneous heat transfer is determined by the instantaneous velocity), an unsteady, nonlinear model for the heat release fluctuations is obtained:

where $k$ is the air thermal conductivity, $c_{p}$ the specific heat per unit volume, assumed to be constant, $T_{w}$ , $L_{w}$ and $d_{w}$ are the temperature, total length and diameter of the hot-wire respectively, and $\unicode[STIX]{x1D70F}$ a time delay. The heat release time-delayed response with respect to acoustic velocity fluctuations models the low-frequency response of the heater, as found by Lighthill (Reference Lighthill1954). The non-dimensional fluctuations are given by:

where we have defined the non-dimensional parameter $K\equiv \sqrt{2\unicode[STIX]{x03C0}c_{p}d_{w}\overline{u}/k}$ .

In Subramanian et al. (Reference Subramanian, Sujith and Wahi2013) a Taylor expansion of the heat release fluctuations (2.3) around $u^{\prime }=0$ was considered for a weakly nonlinear analysis. However, due to the presence of an absolute value in (2.3), such a Taylor series converges to the original heat release expression only for $u^{\prime }>-1$ , as shown in figure 2(a). This is problematic because the fact that $q^{\prime }$ increases when $u^{\prime }<-1$ decreases is the only saturation mechanism, which we want to model in this study. The weakness of a Taylor expansion of this model can be quantitatively shown by calculating the describing function (DF) of King’s model. Imposing harmonic velocity oscillations $u^{\prime }=A\cos (\unicode[STIX]{x1D714}t)$ , the DF is defined as (Gelb & Velde Reference Gelb and Velde1968):

Its gain is shown in figure 2(b) together with the gain of the DF of the Taylor expansion, which does not saturate when the amplitude increases. Note that the unsteady King’s law DF gain first increases up to an amplitude of $A=1$ and then decreases. The initial increase in gain is probably a non-physical feature of the model. It arises upon the introduction of a time-delayed unsteady velocity dependence in King’s law, which was originally derived under steady flow assumptions. Although time-delayed unsteady flow effects provide the correct connection between King’s law at small amplitudes with the linear dynamical theory of Lighthill (Reference Lighthill1954) at low frequencies, it is questionable if these effects correctly capture the nonlinear dynamics of the heater at high amplitudes. Indeed, computational fluid dynamics (CFD) simulations performed by Selimefendigil, Föeller & Polifke (Reference Selimefendigil, Föeller and Polifke2012) on the fully nonlinear unsteady equations of a pulsating flow around a hot-wire showed that the gain of the heat release DF monotonically decreases when the forcing amplitude increases (see figure 2). Furthermore, Witte & Polifke (Reference Witte and Polifke2015) reported that the Reynolds number, $\mathit{Re}$ , based on the wire diameter has a great influence on the linear response of the heater, which, for low values of $\mathit{Re}$ , is very different from the one predicted by Lighthill (Reference Lighthill1954).

Figure 2.

Nonlinear heat release models. (a) Comparison between unsteady King’s law, fifth-order Taylor expansion and polynomial fit approximations. $K$ is fixed to 1. (b) The Describing Function gain of King’s model and its approximations. The saturation mechanism starts at $A=1$ and is due to the abrupt change of sign of the heat release derivative. The Taylor expansion around $u^{\prime }=0$ cannot capture this mechanism. CFD results from Selimefendigil et al. (Reference Selimefendigil, Föeller and Polifke2012) (scaled to match the linear gain) are in qualitative agreement with the least-square fit model that we use in this paper.

Because a smooth, nonlinear, dynamical model for a heater response is not available in the literature, we perform a least-square fit of (2.3) onto a fifth-order polynomial of the form:

In order to ensure that the linear stability of the thermoacoustic system is not affected by the fitting, we constrain the linear coefficient to be $\unicode[STIX]{x1D6FC}_{1}=1/2$ , which is the first-order Taylor coefficient of King’s law (2.3). The $\unicode[STIX]{x1D6FC}$ coefficients depend on the range chosen for the fit. We choose the range $u^{\prime }\in [-2,2]$ , which is wide enough to capture the saturation mechanism up to amplitudes that are much larger than those we will consider in this study. The coefficients obtained with the fit are $\unicode[STIX]{x1D6FC}_{2}=-0.108$ , $\unicode[STIX]{x1D6FC}_{3}=-0.044$ , $\unicode[STIX]{x1D6FC}_{4}=+0.059$ , $\unicode[STIX]{x1D6FC}_{5}=-0.012$ . From figure 2(a), one can see that the least-square fit model has the same linear behaviour as the unsteady King’s law and a qualitatively similar nonlinear behaviour, with a smooth saturation mechanism. The saturation starts at values smaller than $u^{\prime }=1$ , which is consistent with the experimental observations of Heckl (Reference Heckl1990). Also, figure 2(b) shows that the gain of the fitted model decreases monotonically with the amplitude forcing, which is consistent with the nonlinear unsteady calculations of Selimefendigil et al. (Reference Selimefendigil, Föeller and Polifke2012). In the following, we will find that oscillations saturate with amplitudes $|u^{\prime }|<1$ , for which our fit is in good agreement with CFD results. If an appropriate model for the nonlinear dynamic response of the heater is provided, the $\unicode[STIX]{x1D6FC}$ coefficients can be obtained via a Taylor expansion up to the desired order.

Figure 3.

(a) Gain of the hot-wire linear dynamic response as derived by Lighthill (Reference Lighthill1954) (dashed lines) and fitted response (solid line). (b) Block diagram representation of the modelled thermoacoustic dynamics. The linear response $L(s)$ , which combines the acoustic and linear heat responses, is cast into state-space form. When combined with the nonlinear heat release saturation mechanism, it forms a closed-loop nonlinear thermoacoustic system.

Note that the gain of (2.3) is constant for all frequencies, since it is a static heat release model (i.e. $q^{\prime }$ does not depend on $\text{d}u^{\prime }/\text{d}t$ ). A physical mechanism that damps high-frequency modes has to be included, otherwise non-physical high-frequency thermoacoustic oscillations are likely to arise. Lighthill (Reference Lighthill1954) derived analytically the linear dynamic response, $G$ , of a hot-wire subject to velocity perturbations in the low- and high-frequency limits, showing that its behaviour is analogous to that of a low-pass filter. Lighthill’s results for the gain amplitude against frequency are reported in figure 3(a). The low-pass filter cutoff frequency is a function of the wire diameter and mean flow intensity. A gain decrease with frequency provides the necessary stabilisation of high-frequency modes. The total response of the heater is therefore given by the product between a linear dynamic response and a nonlinear amplitude saturation, which is a Wiener–Hammerstein model. Note that here we are retaining the heater phase response due to the time delay in the nonlinear heating element. This is a necessary condition here for the model to exhibit subcritical bifurcations. The latter would not be possible if the phase response had no amplitude dependence, because the gain monotonically decreases with the amplitude.

Lighthill’s dynamic response can be combined with the linear acoustic response $H(s)$ , by defining the total linear response $L(s)\equiv H(s)|G(s)|$ , as shown in figure 3(b). Thus, the thermoacoustic system is modelled as a Lur’e system (a linear component in a feedback loop with a nonlinear component). We extend the linear response $L$ to the entire Laplace space via the rational approximation $L(s)\approx \boldsymbol{C}(s\unicode[STIX]{x1D644}-\unicode[STIX]{x1D63C})^{-1}\boldsymbol{B}$ , where $\unicode[STIX]{x1D63C}$ is an $N\times N$ matrix, $\unicode[STIX]{x1D644}$ is the identity matrix and $\boldsymbol{B}$ , $\boldsymbol{C}$ are the input (column, $N\times 1$ ) and output (row, $1\times N$ ) vectors, respectively. The size of the state-space model, $N$ , determines the order of approximation. The coefficients of the state-space model are determined in a least-square sense. The location of the poles is optimised recursively starting from the acoustic eigenvalues as initial guesses, by following the procedure described in Gustavsen & Semlyen (Reference Gustavsen and Semlyen1999). We then cast this response in state-space form:

where we have defined the state-space vector $\boldsymbol{x}$ , which can be thought of as a linear combination of the acoustic waves $f$ and $g$ evaluated at the heater location. Further details on the least-square procedure described above are provided in Orchini et al. (Reference Orchini, Illingworth and Juniper2015). If an expression for the heat release is provided, the equations above describe the dynamics of thermoacoustic oscillations.

In the following, variables are presented in a non-dimensional form. We scale lengths with the duct length $L$ , velocities with the mean flow velocity $\overline{u}_{1}$ and time with the characteristic time scale $L/c_{1}$ , so that:

4.1 $O(\unicode[STIX]{x1D716})$ : eigenvalue problem

At order $\unicode[STIX]{x1D716}$ we retrieve the homogeneous linear equations (4.6) for the evolution of $\boldsymbol{x}_{1}$ :

Because the left-hand structure of the equations will be the same at all $\unicode[STIX]{x1D716}$ orders, it is convenient to define the spectral operator:

so that the nonlinear eigenvalue problem in the frequency domain can be rewritten as $\unicode[STIX]{x1D648}_{s}\boldsymbol{x}_{1}=\mathbf{0}$ . This is the same eigenvalue problem as in (3.3). Also, because here we have fixed $K=K^{c}$ , we know that the system is marginally stable (its spectrum is shown in figure 4 b with crosses).

We can simplify the evolution of the dynamical system to the evolution of the marginally stable thermoacoustic eigenmode only, i.e. by ignoring the contribution of the eigenmodes with a negative growth rate, because they will quickly be damped (Sipp & Lebedev Reference Sipp and Lebedev2007). Close to the Hopf bifurcation, we expect the dynamical system to saturate to limit cycle oscillations at the slow time scales, therefore we can write:

where $\unicode[STIX]{x1D714}^{c}$ is the angular frequency of the marginally stable eigenmode, $\boldsymbol{x}_{1}^{W}$ the corresponding right-hand eigenvector and $W(t_{2},t_{4})$ a complex-valued variable which depends on the slow time scales only. $W$ contains information on the amplitude saturation and frequency shift effects caused by nonlinear effects. At the next odd orders, we will explicitly find the dependence of $W$ with respect to the slow time scales.

4.2 $O(\unicode[STIX]{x1D716}^{2})$ : mean shift and second harmonic

At this order we obtain the equations for the evolution of $\boldsymbol{x}_{2}$ . From this order on, forcing terms will appear in the right-hand side of the equations. In general, the forcing terms of order $\unicode[STIX]{x1D716}^{N}$ are due to the nonlinear interactions between the solutions $\boldsymbol{x}_{k}$ at orders $k<N$ , which are known. The only forcing term at this order is $\unicode[STIX]{x1D6FC}_{2}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))^{2}$ , which is due to the interaction of $\boldsymbol{x}_{1}$ with itself. By using the expression (4.8) for $\boldsymbol{x}_{1}$ , we obtain the inhomogeneous linear equation:

where

The superscripts are used to classify the forcing terms by their dependence on the complex amplitudes $W$ . This forcing is composed of two components: a steady forcing with zero frequency (due to the interaction between the eigenmode $\boldsymbol{x}_{1}^{W}$ and its complex conjugate) and second harmonic contributions at frequency $2\unicode[STIX]{x1D714}^{c}$ (due to the interaction between the eigenmode $\boldsymbol{x}_{1}^{W}$ and itself). These are not resonant terms, because the spectrum of the linear operator does not contain $2\unicode[STIX]{x1D714}^{c}$ or $0$ as eigenvalues (see figure 4 b) and (4.9) can be readily solved.

We look for a steady-state solution $\boldsymbol{x}_{2}$ which has the same shape as the forcing, by using the ansatz

Substituting the latter into (4.9), taking the Laplace transform (with respect to the fast time scale $t_{0}$ ), and matching the terms according to their amplitude dependence, we obtain the sets of linear equations:

The matrices $\unicode[STIX]{x1D648}_{2\text{i}\unicode[STIX]{x1D714}^{c}}$ and $\unicode[STIX]{x1D648}_{0}$ are non-singular and can be inverted, yielding the solutions at the various amplitude levels of $\boldsymbol{x}_{2}$ .

In particular, $\boldsymbol{x}_{2}^{W^{2}}$ (and its complex conjugate (c.c.)) describes second harmonic oscillations, whereas $\boldsymbol{x}_{2}^{|W|^{2}}$ , having zero frequency, will cause a shift in the mean acoustic level, caused by the presence of even terms in the expansion of the nonlinear heater element response. This is a well-known effect in hydrodynamics, where zero frequency corrections are due to quadratic terms arising from the nonlinear convective term of the Navier–Stokes equations. A weakly nonlinear expansion allows distinction between the base flow (solution of the steady Navier–Stokes equations) and the mean flow (time-averaged solution of the unsteady Navier–Stokes equations) (Sipp & Lebedev Reference Sipp and Lebedev2007; Meliga, Chomaz & Sipp Reference Meliga, Chomaz and Sipp2009). These effects are not found if the nonlinearity expansion contains only odd terms. This is often the case for low-order thermoacoustic modelling (Noiray, Bothien & Schuermans Reference Noiray, Bothien and Schuermans2011; Noiray & Schuermans Reference Noiray and Schuermans2013; Ghirardo et al. Reference Ghirardo, Juniper and Moeck2015), although some experimental evidence of acoustic level mean shifts can be found in the literature (Flandro, Fischbach & Majdalani Reference Flandro, Fischbach and Majdalani2007).

4.3 $O(\unicode[STIX]{x1D716}^{3})$ : third harmonic and saturation

At this order we obtain the equations for the evolution of $\boldsymbol{x}_{3}$ :

The explicit expressions of the forcing terms are reported in § C.1. The derivatives of the slow time scales may be rewritten as:

Resonant forcings arise at this order, with angular frequency $\unicode[STIX]{x1D714}^{c}$ , which act on two different amplitude levels, $W$ and $|W|^{2}W$ . A solvability condition, known as the Fredholm alternative (Oden & Demkowicz Reference Oden and Demkowicz2010), needs to be satisfied for a solution to exist. It requires the sum of the resonant forcing terms to be orthogonal to the kernel (nullspace) of the (singular) adjoint operator $\unicode[STIX]{x1D648}_{\,\text{i}\unicode[STIX]{x1D714}^{c}}^{\dagger }$ (Sipp & Lebedev Reference Sipp and Lebedev2007; Meliga et al. Reference Meliga, Chomaz and Sipp2009). This generalises the idea of cancelling the secular terms used for weakly nonlinear analysis of scalar problems.

The adjoint matrix $\unicode[STIX]{x1D648}_{\,\text{i}\unicode[STIX]{x1D714}^{c}}^{}$ is defined through the scalar product:

and corresponds to the Hermitian of the direct matrix $\unicode[STIX]{x1D648}_{\,\text{i}\unicode[STIX]{x1D714}^{c}}^{}$ . The latter has a zero eigenvalue, which corresponds to the direct eigenvector $\boldsymbol{x}_{1}^{W}$ . Because an adjoint matrix has eigenvalues which are complex conjugate of those of its direct matrix, $\unicode[STIX]{x1D648}_{\,\text{i}\unicode[STIX]{x1D714}^{c}}^{\dagger }$ has a zero eigenvalue, and its kernel is spanned by the adjoint eigenvector $\boldsymbol{x}_{1}^{\dagger }$ only. This can be calculated as the Hermitian of the left eigenvector of the operator $\unicode[STIX]{x1D648}_{\,\text{i}\unicode[STIX]{x1D714}^{c}}$ corresponding to the eigenvalue 0. The solvability condition therefore requires:

where the right-hand terms in the bracket are all the resonant forcings, and we have defined the matrix $\unicode[STIX]{x1D64B}\equiv (\unicode[STIX]{x1D644}+\unicode[STIX]{x1D70F}\boldsymbol{B}\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{C}\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{c}\unicode[STIX]{x1D70F}})$ . By rearranging (4.16), we obtain:

where the complex values $\unicode[STIX]{x1D706}_{3}$ , $\unicode[STIX]{x1D708}_{3}$ , known as the Landau coefficients, are defined by:

The values we found for the Landau coefficients when $K<K^{c}$ are $\unicode[STIX]{x1D706}_{3}=-0.0659-0.1523\text{i}$ and $\unicode[STIX]{x1D708}_{3}=0.0007-0.0048\text{i}$ . The dependence of the Landau coefficients on the number of modes retained in the acoustic state-space model is shown in figure 6 a. Note that, in this case, we find that the sign of $\text{Re}(\unicode[STIX]{x1D708}_{3})$ with one mode is different from the sign of the saturated value, containing multiple modes, whereas the sign of $\text{Re}(\unicode[STIX]{x1D706}_{3})$ does not change. The sign of the ratio of the latter values is important, as it distinguishes between sub- and supercritical bifurcations. This shows that the nature of the bifurcation predicted with a single-mode approximation may be different from the actual response of the system.

We then seek a solution $\boldsymbol{x}_{3}$ via the ansatz:

By using the relation (4.17), we match the solution and forcing terms of (4.13) by their dependence on the amplitude $W$ , and obtain the sets of linear equations:

Although the matrix $\unicode[STIX]{x1D648}_{\,\text{i}\unicode[STIX]{x1D714}^{c}}$ is singular, the values of the Landau coefficients guarantee that solutions for $\boldsymbol{x}_{3}^{W}$ and $\boldsymbol{x}_{3}^{|W|^{2}W}$ exist. They can be calculated, e.g. by using the pseudo-inverse matrix of $\unicode[STIX]{x1D648}_{\,\text{i}\unicode[STIX]{x1D714}^{c}}$ . These solutions provide a nonlinear correction to the shape of the linearly unstable mode. Equation (4.20c ), instead, can be readily solved by inverting the matrix $\unicode[STIX]{x1D648}_{3\text{i}\unicode[STIX]{x1D714}^{c}}$ , which is non-singular. $\boldsymbol{x}_{3}^{W^{3}}$ accounts for third harmonic contributions to the oscillatory solution.

4.4 $O(\unicode[STIX]{x1D716}^{4})$ : mean shift and fourth harmonic

At this order we obtain the equations for the evolution of $\boldsymbol{x}_{4}$ :

None of the forcing terms resonate. Their expressions are provided in § C.2. Using the ansatz:

we can readily calculate the solutions

which provide contributions to the mean acoustic level shift and second and fourth harmonic oscillations.

4.5 $O(\unicode[STIX]{x1D716}^{5})$ : fifth harmonic and saturation

At this order we obtain the equations for the evolution of $\boldsymbol{x}_{5}$ .

In this section, we will not explicitly calculate the solution $\boldsymbol{x}_{5}$ , but we will only derive the dependence of the amplitude $W$ with respect to the slow time scale $t_{4}$ . This is achieved by applying the Fredholm alternative solvability condition on the resonant forcings at frequency $\unicode[STIX]{x1D714}^{c}$ which appear on the right-hand side of (4.29).

By imposing that the resonant terms are orthogonal to the kernel of the adjoint operator $\unicode[STIX]{x1D648}_{\,\text{i}\unicode[STIX]{x1D714}^{c}}^{\dagger }$ , we obtain:

The resonant forcing expressions can be found in § C.3. Equation (4.30) can be simplified into the Stuart–Landau equation:

where the Landau coefficients are defined by:

When using 12 acoustic modes in the state space, the values we obtain for these coefficients when $K<K^{c}$ are $\unicode[STIX]{x1D706}_{5}=-0.0202-0.0184\text{i}$ , $\unicode[STIX]{x1D708}_{5}=0.0072+0.0024\text{i}$ , $\unicode[STIX]{x1D707}_{5}=-0.0014-0.0007\text{i}$ . The coefficient $\unicode[STIX]{x1D708}_{5}$ changes sign across the Hopf location.

5.1 Time domain simulations

Time domain simulations are performed by initialising the integration to a state which is slightly perturbed from the fixed point solution. Because the equations are time delayed, the initial state covers the history of the system for a time $-\unicode[STIX]{x1D70F}\leqslant t\leqslant 0$ . We start from a value of $K<K^{c}$ , and then increase it in steps of $\unicode[STIX]{x0394}K$ . Until $K\leqslant K^{c}$ , the initial perturbations are damped and the system converges to fixed points solutions. At $K=1.43$ , the oscillations start growing and converge towards a limit cycle attractor with a large amplitude (see figure 8 a). This solution is used to initialise the subsequent integrations, for which the control parameter $K$ is varied in both directions in steps of $\pm \unicode[STIX]{x0394}K$ . The amplitude of the oscillations gets larger as $K$ increases. On the other hand, the amplitude of the velocity fluctuations decreases smoothly for $K<K^{c}$ , until we reach the fold location at $K^{F}=1.421$ . At this location, the initialised highly perturbed initial state decays to the fixed point solution, as shown in figure 8(b). The largest value of the control parameter we consider is $K_{max}=1.424$ , for which the parameter expansion of the weakly nonlinear analysis is $\unicode[STIX]{x1D716}_{max}=\sqrt{K_{max}-K_{c}}=0.039\ll 1$ .

To validate the linear analysis, the growth/decay rates close to the Hopf/fold locations are extracted from the time series by using linear regression on the logarithm of the fluctuations amplitude. The latter is obtained using the Hilbert transform intervals of the time series, corresponding to the linear amplitude regime. The growth and decay rates, reported in figure 8, can be compared with those we obtain when solving the eigenvalue problem (3.3) at the same locations. These values are $\unicode[STIX]{x1D70E}=4.92\times 10^{-4}$ and $\unicode[STIX]{x1D70E}=-9.97\times 10^{-5}$ respectively, and are in good agreement with the time-marching results.

5.2 Numerical continuation

To validate the weakly nonlinear analysis, we compare the oscillations’ amplitudes at the various frequency components with those predicted by numerical continuation of limit cycles. This is because the latter also yields information about the unstable limit cycles, which are more difficult to investigate with time-marching simulations. The bifurcation diagram shown in figure 9(a) shows the amplitude of the oscillations at the fundamental frequency as predicted from the weakly nonlinear analysis at fifth order and limit cycle continuation. As one can see, no significant difference between the third and fifth order is observed close the Hopf point $K^{C}$ . However, the fifth-order analysis significantly improves the predictions of the third-order analysis as (i) it predicts the existence of a bistable region; (ii) it accurately locates the fold point; (iii) it predicts the amplitude of the fundamental harmonic oscillations both for stable and unstable limit cycles.

Figure 9.

Comparison between the weakly nonlinear analysis at various orders (lines) and numerical continuation results (circles). The bistable region is highlighted in grey. Solid and dashed lines indicate stable and unstable solutions, respectively. (a) Bifurcation diagram showing the amplitude of the oscillations at the resonant frequency. The fifth-order analysis is in good agreement with the exact solution. (b) Frequency shift $\unicode[STIX]{x1D6E5}\unicode[STIX]{x1D714}$ of the oscillations with respect to the marginally stable frequency at the Hopf bifurcation.

We highlight that the weakly nonlinear analysis is accurate in the vicinity of the critical point, meaning $\unicode[STIX]{x1D716}_{max}\ll 1$ . Note that it could be possible for the width of the subcritical region to be wider, so that larger values of $\unicode[STIX]{x1D716}$ would be needed to locate the fold point. In this case, a higher-order expansion would be needed for the weakly nonlinear analysis to reasonably approximate the solution away from the critical point. This is not shown in this manuscript, but is discussed in Orchini (Reference Orchini2016), where the weakly nonlinear expansion is performed up to seventh order.

Figure 10.

Magnitude of the harmonic contributions in the spectrum of the stable oscillatory solutions. Results from weakly nonlinear analysis (lines) and amplitude of fast Fourier transform of time domain simulations (circles) are compared.

Further information can be extracted from the weakly nonlinear analysis, in particular the frequency shift of limit cycle oscillations with respect to that of the marginally stable eigenvalue (figure 9 b), and the amplitude of the harmonics in the spectrum of the oscillatory solution (figure 10). The zeroth harmonic corresponds to the mean shift in the acoustic level previously discussed, which is indeed observed in time domain simulations, due to the fact that the nonlinearity is not odd. All these results show a good agreement with those obtained from the exact solution of the full nonlinear system.

5.3 Computational cost

One of the attractive characteristics of the weakly nonlinear analysis is the low computational cost compared to time-marching solutions or continuation algorithms. The expansion of the equations can be tedious at high orders, but the process can be automatised with symbolic solvers. The information needed for the analysis is (i) a marginally stable solution, (ii) the direct and adjoint eigenvectors corresponding to the marginally stable solution and (iii) the Taylor expansion of the nonlinear heat release fluctuations. Once these components have been found, the calculation of a bifurcation diagram reduces to a series of matrix inversions and multiplications. This is affordable even for systems with a large number of degrees of freedom as shown by Sipp & Lebedev (Reference Sipp and Lebedev2007), Meliga et al. (Reference Meliga, Chomaz and Sipp2009) for the analysis of hydrodynamic bifurcations. Thus, the entire bifurcation diagram can be constructed by means of weakly nonlinear analysis, which would require many more measurements if determined experimentally or numerically. Moreover, the dynamic evolution of the oscillations associated with the global mode can be obtained for various values of the bifurcation parameter, by direct integration of the Stuart–Landau equation.

Techniques such as numerical continuation (Waugh, Kashinath & Juniper Reference Waugh, Kashinath and Juniper2014) require time marching for several cycles in order to converge to a limit cycle. Direct numerical simulations require even more cycles and can converge only to stable oscillations. For our low-order model and the results presented above, we find that the weakly nonlinear analysis is 20 times faster than numerical continuation, and three orders of magnitudes faster than brute force time marching. Note that, in the present study, the numerical continuation method takes only a few minutes on a desktop computer because we are using a low-order network model. However, numerical continuation quickly becomes very expensive when, e.g. a combustion process is introduced in the model. For example, Waugh et al. (Reference Waugh, Kashinath and Juniper2014) considered a thermoacoustic system with a simple model for a flame, the $G$ -equation. It required 14 000 CPU hours to build a bifurcation diagram with numerical continuation. Using weakly nonlinear analysis, one only has to calculate the direct and adjoint eigenstates around the Hopf point, and perform a small number of matrix inversions (ten matrix inversions were needed for the analysis presented here). These operations can be expensive as well for systems with many degrees of freedom, but this is affordable nowadays (Sipp & Lebedev Reference Sipp and Lebedev2007) and it is much cheaper than any method involving time marching of the equations.