1 Introduction
Thermoacoustic oscillations, which are caused by the coupling between unsteady acoustics and unsteady heat release in combustion chambers, can threaten the operability of combustion systems such as rocket engines and gas turbines. If the system is susceptible to thermoacoustic instabilities, the amplitude of oscillations initially increases exponentially and finally saturates to self-sustained oscillations. In many situations, thermoacoustic oscillations are undesirable: they produce noise, structural vibrations and impose limits on the operating conditions, which can reduce the system efficiency (Lieuwen & Yang Reference Lieuwen and Yang2005; Culick Reference Culick2006).
Linear stability analysis can predict the onset of thermoacoustic oscillations when a control parameter of the system is varied. The change in behaviour from linear stability (small perturbations decay) to linear instability (small perturbations grow) is termed bifurcation. The linear analysis cannot predict the final amplitude of the bifurcated state of the system and the transition process, because it describes only the initial phase of the perturbations’ evolution. Nonlinear and non-normal effects can lead to a variety of complex behaviours, such as triggering and bistability (Juniper Reference Juniper2011). Moreover, nonlinear effects determine the nature of the final state, which can be a fixed point, a limit cycle or a more complex solution. All these phenomena have been observed experimentally in thermoacoustic systems (Noiray et al. Reference Noiray, Durox, Schuller and Candel2008; Gotoda et al. Reference Gotoda, Nikimoto, Miyano and Tachibana2011; Kabiraj & Sujith Reference Kabiraj and Sujith2012; Jegadeesan & Sujith Reference Jegadeesan and Sujith2013; Rigas et al. Reference Rigas, Jamieson, Li and Juniper2015a ).
Weakly nonlinear analysis is a nonlinear method capable of tracking the evolution of the amplitude of the oscillations. It is based on an asymptotic expansion of the governing equations in the vicinity of a bifurcation point. This approach provides an analytical description of the perturbation dynamics, which is exact up to the order of the truncated expansion. In general, the solution is calculated as a superposition of one or more spatial modes with a time-dependent amplitude. The temporal evolution of the amplitude is reduced to an ordinary differential equation (ODE) – appearing as a Stuart–Landau equation – for every linearly unstable mode. Solving the latter equation is much faster than time marching the full nonlinear system, and also provides physical insight into the nonlinear interactions between the modes.
Weakly nonlinear analysis has been widely used in hydrodynamics to study the transition of globally unstable flows and derive low-order models of the Navier–Stokes equations around bifurcation points (Chomaz Reference Chomaz2005). In the simple case of cylinder flow, which undergoes a supercritical Hopf bifurcation at a diameter-based Reynolds number
$\mathit{Re}\approx 46$
, the Stuart–Landau equation accurately captures the amplitude of the most unstable global shedding mode close to the threshold of bifurcation (Landau Reference Landau1944; Provansal, Mathis & Boyer Reference Provansal, Mathis and Boyer1987). Sipp & Lebedev (Reference Sipp and Lebedev2007) showed how the Stuart–Landau equation can be derived from the Navier–Stokes equations using global stability analysis and a multiple time scale expansion. As a consequence of the global character of the analysis, adjoint methods were required to identify the Landau coefficients of the model.
Weakly nonlinear tools have been applied also to thermoacoustic systems. Culick (Reference Culick2006) used the method of averaging to derive the amplitude evolution for thermoacoustic models with one or two oscillating modes. In the same framework, Juniper (Reference Juniper2012) described how the averaged quantities can be connected to the Flame Describing Function methodology. Ghirardo, Juniper & Moeck (Reference Ghirardo, Juniper and Moeck2015) applied the method of averaging to azimuthal thermoacoustic instabilities, in which two counter-rotating azimuthal modes with the same frequency are simultaneously unstable. Subramanian, Sujith & Wahi (Reference Subramanian, Sujith and Wahi2013) used the method of multiple scales to derive a Stuart–Landau equation at third order describing the evolution of the oscillation amplitude in a Rijke tube. The Landau coefficients showed that the Hopf bifurcation was subcritical, and the low-amplitude limit cycles arising close to the Hopf point were unstable, in agreement with experimental studies.
The unstable limit cycles arising from subcritical bifurcations in a Rijke tube may undergo a fold bifurcation and create a region of bistability (Ananthkrishnan, Deo & Culick Reference Ananthkrishnan, Deo and Culick2005; Juniper Reference Juniper2011; Jegadeesan & Sujith Reference Jegadeesan and Sujith2013). In the weakly nonlinear analysis performed by Subramanian et al. (Reference Subramanian, Sujith and Wahi2013), however, the expansion was truncated at third order. Therefore, the fold point, the amplitude of stable limit cycles, and the bistable region could not be predicted by their weakly nonlinear methods, because this would require expansion to at least fifth order. In their case, however, even expansion at higher order would not have captured the fold point for reasons explained in § 2.2. It is also worth noticing that in the weakly nonlinear studies on the Rijke tube (Juniper Reference Juniper2012; Subramanian et al. Reference Subramanian, Sujith and Wahi2013) only one Galerkin mode was used to describe the dynamics. This can be a rough approximation for some thermoacoustic networks, because considering only one acoustic mode may alter the nature and amplitude of the oscillations, as was discussed by Jahnke & Culick (Reference Jahnke and Culick1994), Ananthkrishnan et al. (Reference Ananthkrishnan, Deo and Culick2005).
In this paper, we perform a high-order weakly nonlinear analysis of thermoacoustic oscillations in a Rijke tube close to a subcritical Hopf bifurcation. The high-order expansion allows us to obtain analytical expressions for the location of the bistable region and the amplitude of both unstable and stable limit cycles. In our analysis, this is achieved by using a wave-based approach when solving the linear acoustic equations. It provides a more general description of a thermoacoustic network, and enables us to straightforwardly include temperature and area variations in the analysis. We also retain the contribution of multiple acoustic modes to the dynamics of the thermoacoustic system. This corresponds to approaching the problem in a global framework, and adjoint methods are required to calculate the Landau coefficients through solvability conditions (Sipp & Lebedev Reference Sipp and Lebedev2007).
The paper is organised as follows. In § 2.1 we present the Rijke tube thermoacoustic set-up and the wave-based governing equations. In § 2.2 the nonlinear heat release model adopted for this study is discussed. In § 3 we perform a linear stability analysis of the system and we identify the location of Hopf bifurcations, when the heat release power is used as a control parameter. In § 4 we present the theoretical framework for the weakly nonlinear analysis, deriving in detail the equations for the amplitude evolution of the dominant mode up to fifth order. In § 5 we validate our weakly nonlinear results against the exact solutions of the fully nonlinear equations, obtained with a continuation algorithm method and time domain simulations. A good agreement between the weakly nonlinear and fully nonlinear analysis is observed in the oscillation amplitude, harmonics contributions and frequency shift. Finally, in § 6 we summarise our findings and discuss possible future applications.
2 Thermoacoustic modelling
Figure 1. Horizontal Rijke tube model. Subscripts
$_{1}$
and
$_{2}$
indicate flow and wave properties in the upstream and downstream ducts, respectively. The intensity of the heat release fluctuations will be used as a control parameter.
The configuration considered in this study is that of a horizontal Rijke tube, as shown in figure 1. It has been extensively considered by many authors (Matveev Reference Matveev2003; Juniper Reference Juniper2011; Magri & Juniper Reference Magri and Juniper2013; Subramanian et al. Reference Subramanian, Sujith and Wahi2013; Mariappan, Sujith & Schmid Reference Mariappan, Sujith and Schmid2015) for the analysis of thermoacoustic phenomena. However, the weakly nonlinear analyses of thermoacoustic oscillations presented in these studies have been approximated by considering a single pair of Galerkin modes for the acoustic response. It was shown by Jahnke & Culick (Reference Jahnke and Culick1994), Ananthkrishnan et al. (Reference Ananthkrishnan, Deo and Culick2005), Kashinath, Waugh & Juniper (Reference Kashinath, Waugh and Juniper2014) that considering only one acoustic mode may alter the amplitude and type of thermoacoustic oscillations. Also, mean flow and temperature jump effects are often neglected, although their presence affects the thermoacoustic eigenmodes and the stability of the system. Although our model remains low order, we consider a wave-based approach which naturally yields a more realistic description of the acoustic response of the system (Dowling Reference Dowling1995; Orchini, Illingworth & Juniper Reference Orchini, Illingworth and Juniper2015), and is easily scalable to more complex acoustic networks, which is important when considering, for example, gas turbines. As is customary for the analysis of thermoacoustic oscillations in gas turbines, we linearise the acoustic equations and retain the heater response as the only nonlinear element (Dowling Reference Dowling1997; Culick Reference Culick2006).
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:
$$\begin{eqnarray}\unicode[STIX]{x1D648}\left[\begin{array}{@{}c@{}}g_{1}\\ f_{2}\end{array}\right]=\left[\begin{array}{@{}cc@{}}M_{11} & M_{12}\\ M_{21} & M_{22}\end{array}\right]\left[\begin{array}{@{}c@{}}g_{1}\\ f_{2}\end{array}\right]=\left[\begin{array}{@{}c@{}}0\\ q^{\prime }/(A_{1}c_{1})\end{array}\right]\!,\end{eqnarray}$$
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:
$$\begin{eqnarray}Q=k(T_{w}-T_{1})L_{w}\left(1+\sqrt{\frac{2\unicode[STIX]{x03C0}c_{p}d_{w}}{k}|\overline{u}+u^{\prime }(t-\unicode[STIX]{x1D70F})|}\right)\!,\end{eqnarray}$$
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:
$$\begin{eqnarray}q^{\prime }=K\left(\sqrt{|1+u^{\prime }(t-\unicode[STIX]{x1D70F})|}-1\right)\!,\end{eqnarray}$$
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):
$$\begin{eqnarray}\text{DF}(q^{\prime })\equiv \frac{K}{A\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\left(\sqrt{|1+A\cos (\unicode[STIX]{x1D703}-\unicode[STIX]{x1D714}\unicode[STIX]{x1D70F})|}-1\right)\text{e}^{\text{i}\unicode[STIX]{x1D703}}\,\text{d}\unicode[STIX]{x1D703}.\end{eqnarray}$$
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:
$$\begin{eqnarray}q^{\prime }(t)=K\mathop{\sum }_{n=1}^{5}\unicode[STIX]{x1D6FC}_{n}{u^{\prime }}^{n}(t-\unicode[STIX]{x1D70F}).\end{eqnarray}$$
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:
$$\begin{eqnarray}\displaystyle \dot{\boldsymbol{x}} & = & \displaystyle \unicode[STIX]{x1D63C}\boldsymbol{x}+\boldsymbol{B}q^{\prime }\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle u^{\prime } & = & \displaystyle \boldsymbol{C}\boldsymbol{x},\end{eqnarray}$$
$\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:
$$\begin{eqnarray}x_{h}/L\rightarrow x_{h},\quad u^{\prime }/\overline{u}\rightarrow u^{\prime },\quad \unicode[STIX]{x1D70F}c_{1}/L\rightarrow \unicode[STIX]{x1D70F},\quad \unicode[STIX]{x1D714}L/c_{1}\rightarrow \unicode[STIX]{x1D714}.\end{eqnarray}$$
3 Linear stability analysis
By coupling the acoustic response (2.6b
) with the heat release model (2.5), we obtain a nonlinear, delayed dynamical system that describes the evolution of thermoacoustic oscillations. The dynamics of the state variables
$\boldsymbol{x}$
is coupled upon the introduction of the heat element through the acoustic velocity
$u^{\prime }=\boldsymbol{C}\boldsymbol{x}$
:
$$\begin{eqnarray}\dot{\boldsymbol{x}}(t)=\unicode[STIX]{x1D63C}\boldsymbol{x}(t)+K\boldsymbol{B}\mathop{\sum }_{n=1}^{5}\unicode[STIX]{x1D6FC}_{n}(\boldsymbol{C}\boldsymbol{x}(t-\unicode[STIX]{x1D70F}))^{n},\end{eqnarray}$$
where
$K$
is the control parameter, which determines the intensity of the heat release fluctuations. The matrix
$\unicode[STIX]{x1D63C}$
and the column and row vectors
$\boldsymbol{B}$
and
$\boldsymbol{C}$
are functions of the state-space model order of approximation and of the flame location
$x_{h}$
. The latter is fixed to
$x_{h}=0.13$
, which is the location at which our thermoacoustic network is most prone to thermoacoustic oscillations at its lowest eigenfrequency when the time delay is chosen to be
$\unicode[STIX]{x1D70F}=0.04$
. This was obtained by searching for the maximum value of the Rayleigh index
$\mathit{Ry}\equiv \int _{0}^{T}p_{x_{h}}^{\prime }(t)q^{\prime }(t)\,\text{d}t$
, assuming harmonic oscillations and using (3.1) in order to relate heat release fluctuations to velocity fluctuations.
Fixing all the other parameters, there are specific values of
$K$
, called critical points, at which the system is marginally stable, meaning that all the eigenvalues of the spectrum of the linear operator have a negative growth rate except for a complex conjugate pair, which has a zero growth rate. At these critical points Hopf bifurcations occur, and the linear stability of the system changes across them. In order to perform a weakly nonlinear analysis, we need to first locate critical points and then expand the governing equations around them. This will yield the amplitude and frequency of limit cycle oscillations that occur after the bifurcation.
The fixed point
$\boldsymbol{x}=\mathbf{0}$
is a solution of the dynamical system (3.1) for any value of
$K$
. Its stability is determined by the eigenvalues of the linearised system:
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\boldsymbol{x}}(t)=\unicode[STIX]{x1D63C}\boldsymbol{x}(t)+\unicode[STIX]{x1D6FC}_{1}K\boldsymbol{B}\boldsymbol{C}\boldsymbol{x}(t-\unicode[STIX]{x1D70F}). & \displaystyle\end{eqnarray}$$
Taking the Laplace transform we have:
$$\begin{eqnarray}\left(s\unicode[STIX]{x1D644}-\unicode[STIX]{x1D63C}-\unicode[STIX]{x1D6FC}_{1}K\boldsymbol{B}\boldsymbol{C}\text{e}^{-s\unicode[STIX]{x1D70F}}\right)\hat{\boldsymbol{x}}=\mathbf{0},\end{eqnarray}$$
where
$s\equiv \unicode[STIX]{x1D70E}+\text{i}\unicode[STIX]{x1D714}$
is the non-dimensional Laplace variable and
$\unicode[STIX]{x1D644}$
the identity matrix. Note that, because
$\boldsymbol{B}$
and
$\boldsymbol{C}$
are column and row vectors, respectively, their (outer) product results in a matrix. The values of
$s$
for which the determinant of the linear operator vanishes are the thermoacoustic eigenvalues. The nonlinear eigenvalue problem (3.3) is solved iteratively while varying
$K$
until a marginally stable solution is found.
Figure 4. (a) Neutral line in the
$K{-}\unicode[STIX]{x1D70F}$
plane. Along it, the eigenmode with the smallest eigenfrequency has zero growth rate. Colours refer to the frequency of the marginally stable mode. Black lines indicate that mode(s) with a higher frequency have a positive growth rate. (b) Spectrum of the acoustic (circles) and thermoacoustic (crosses) system for the set of parameters chosen for the weakly nonlinear analysis, marked with a circle in the left panel. The paths of the eigenvalues from their acoustic values (
$K=0$
, dots) to their thermoacoustic values (
$K=K^{c}$
, crosses) are shown as dotted lines.
We first perform a parametric study in the
$K{-}\unicode[STIX]{x1D70F}$
plane to identify a set of Hopf bifurcations of the thermoacoustic mode with the smallest eigenfrequency (see figure 4
a). This is achieved using the open-source package DDE-BIFTOOL (Engelborghs, Luzyanina & Roose Reference Engelborghs, Luzyanina and Roose2002; Sieber et al.
Reference Sieber, Engelborghs, Luzyanina, Samaey and Roose2015), tracking the critical values of the heater power
$K^{c}$
at which the growth rate associated with the smallest eigenfrequency is zero. However, we find that, for some of these solutions, mode(s) with a higher frequency have a positive growth rate along the neutral lines of the first mode. This is physically possible when the values of
$K^{c}$
or
$\unicode[STIX]{x1D70F}$
are large. We show these solutions with black lines in figure 4(a). At these locations the system is not marginally stable, and the theory presented in the following cannot be applied.
Figure 5. Convergence of the marginally stable angular frequency
$\unicode[STIX]{x1D714}^{c}$
(a) and critical power
$K^{c}$
(b) with respect to the number of modes considered in the state-space approximation. A single-mode approximation does not accurately capture the correct response of the system.
A typical value for the time delay can be estimated from Lighthill’s (1954) theory,
$\unicode[STIX]{x1D70F}=0.2\,d_{w}/\overline{u}$
. For our configuration, this yields a non-dimensional value of the time delay of order
$10^{-2}$
. In the following, we will fix
$\unicode[STIX]{x1D70F}=0.04$
, for which the critical power is
$K^{c}=1.42$
. Figure 4(b) shows the spectrum of the acoustic and thermoacoustic systems at this Hopf bifurcation.
Note that the frequencies of the thermoacoustic eigenmodes are close to the acoustic ones, but the growth rates are significantly shifted as
$K$
is increased. The order of approximation of the state space – which is twice the number of modes considered – has to be large enough in order to capture the thermoacoustic response correctly, as was shown by Ananthkrishnan et al. (Reference Ananthkrishnan, Deo and Culick2005), Kashinath et al. (Reference Kashinath, Waugh and Juniper2014). A single-mode approximation of the acoustic response (as considered by Juniper (Reference Juniper2011), Subramanian et al. (Reference Subramanian, Sujith and Wahi2013) for weakly nonlinear expansions) is not able to capture the system response accurately in this case. Figure 5 shows that the marginally stable eigenfrequency and heat power obtained with a single-mode approximation are different from the ones obtained with multiple modes. The number of modes used can also affect the type of bifurcation predicted by the weakly nonlinear expansion, as we will discuss in § 4.3. In this study, we will retain 12 modes when performing the weakly nonlinear analysis. Thus, our thermoacoustic system is composed of 24 coupled differential equations. As a consequence, the use of adjoint methods is required to obtain solvability conditions of the weakly nonlinear equations, as discussed in § 4.3.
4 Weakly nonlinear analysis
We now perturb the bifurcation parameter from the critical point,
$K=K^{c}+\unicode[STIX]{x0394}K$
, with
$|\unicode[STIX]{x0394}K|\ll 1$
. After the Hopf bifurcation, the system is linearly unstable and oscillations will grow and can saturate to limit cycle oscillations due to nonlinear effects. In the case in which the Hopf bifurcation is subcritical, a bistable region exists before the Hopf point in which, by triggering the system, stable solutions with a finite amplitude can be found (Ananthkrishnan et al.
Reference Ananthkrishnan, Deo and Culick2005; Juniper Reference Juniper2011). To calculate the amplitude and frequencies of these oscillations, nonlinear methods are required. We accomplish this with a weakly nonlinear analysis, by expanding the evolution of the dynamical system (3.1) around the Hopf location. We denote with
$0<\unicode[STIX]{x1D716}\ll 1$
a small quantity that quantifies the amplitude of the oscillations close to the Hopf point. We then seek solutions for
$\boldsymbol{x}$
expressed as a power series of
$\unicode[STIX]{x1D716}$
:
$$\begin{eqnarray}\boldsymbol{x}=\unicode[STIX]{x1D716}\boldsymbol{x}_{1}+\unicode[STIX]{x1D716}^{2}\boldsymbol{x}_{2}+\unicode[STIX]{x1D716}^{3}\boldsymbol{x}_{3}+\unicode[STIX]{x1D716}^{4}\boldsymbol{x}_{4}+\unicode[STIX]{x1D716}^{5}\boldsymbol{x}_{5}+O(\unicode[STIX]{x1D716}^{6}).\end{eqnarray}$$
For a subcritical bifurcation, expansion at third order yields unstable limit cycle solutions. However, such oscillations are typically not observable in self-excited experiments, although they can be studied experimentally by forcing the system close to the unstable solutions (Jegadeesan & Sujith Reference Jegadeesan and Sujith2013). For subcritical phenomena however, one is also interested in calculating the amplitude of stable solutions, and in identifying the width of the bistable region. This requires terms of at least order 5 in a weakly nonlinear expansion.
We choose to work with the method of multiple scales. With this method, one assumes that several independent time scales act on the system. One is the fast time scale
$t_{0}$
, at which the oscillations of the marginally stable frequency respond. The slow time scales
$\tilde{t}_{2}$
and
$\tilde{t}_{4}$
are associated with long time saturation or growth processes. The total time derivative therefore reads:
$$\begin{eqnarray}\frac{\text{d}}{\text{d}t}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t_{0}}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\tilde{t}_{2}}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\tilde{t}_{4}}+\cdots .\end{eqnarray}$$
By Taylor expanding the dynamical system (3.1) around the fixed point solution
$\boldsymbol{x}=\mathbf{0}$
at the critical point
$K=K^{c}$
we obtain:
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{x}}{\unicode[STIX]{x2202}t_{0}}+\frac{\unicode[STIX]{x2202}\boldsymbol{x}}{\unicode[STIX]{x2202}\tilde{t}_{2}}+\frac{\unicode[STIX]{x2202}\boldsymbol{x}}{\unicode[STIX]{x2202}\tilde{t}_{4}} & = & \displaystyle \unicode[STIX]{x1D63C}\boldsymbol{x}+\mathop{\sum }_{n=1}^{5}\unicode[STIX]{x1D6FC}_{n}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}(t-\unicode[STIX]{x1D70F}))^{n}\cdots \nonumber\\ \displaystyle & & \displaystyle +\,\mathop{\sum }_{n=1}^{5}\unicode[STIX]{x1D6FC}_{n}\unicode[STIX]{x0394}K\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}(t-\unicode[STIX]{x1D70F}))^{n}.\end{eqnarray}$$
It is important to note that the orders of magnitude of
$O(\boldsymbol{x})=\unicode[STIX]{x1D716}$
and those of
$O(\unicode[STIX]{x0394}K)$
,
$O(\tilde{t}_{2})$
, and
$O(\tilde{t}_{4})$
are not independent. Upon the expansion of the equations, one can show that, at odd orders larger than 1, secular terms (i.e. forcings at resonant frequencies) arise due to nonlinear interactions. Solvability conditions need to be imposed on these forcings, which have to be balanced by contributions arising from slow time scales and control parameter terms (Rosales Reference Rosales2004; Strogatz Reference Strogatz2015). This reasoning leads to balances between the orders of magnitudes of the various terms, which read:
$$\begin{eqnarray}O(\boldsymbol{x}\unicode[STIX]{x0394}K)=\unicode[STIX]{x1D716}^{3},\quad O\left(\frac{\unicode[STIX]{x2202}\boldsymbol{x}}{\unicode[STIX]{x2202}\tilde{t}_{2}}\right)=\unicode[STIX]{x1D716}^{3},\quad O\left(\frac{\unicode[STIX]{x2202}\boldsymbol{x}}{\unicode[STIX]{x2202}\tilde{t}_{4}}\right)=\unicode[STIX]{x1D716}^{5}.\end{eqnarray}$$
We shall then rewrite all the quantities in terms of
$\unicode[STIX]{x1D716}$
by defining
$\unicode[STIX]{x0394}K\equiv \unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D6FF}_{2}$
,
$t_{2}\equiv \unicode[STIX]{x1D716}^{2}\tilde{t}_{2}$
and
$t_{4}\equiv \unicode[STIX]{x1D716}^{4}\tilde{t}_{4}$
. The parameter
$\unicode[STIX]{x1D6FF}_{2}$
can take the values
$\pm 1$
depending on the side of the Hopf point we are investigating. From the definition of
$\unicode[STIX]{x0394}K$
, we also obtain a measure of the expansion parameter in terms of the distance from the Hopf location:
$$\begin{eqnarray}\unicode[STIX]{x1D716}=\sqrt{|K-K^{c}|}.\end{eqnarray}$$
Note that
$\unicode[STIX]{x1D716}$
needs to be
$\ll 1$
. In the following, we will show that the largest value of the bifurcation parameter we need to consider to locate the fold point of a subcritical bifurcation is
$\unicode[STIX]{x1D716}\approx 0.05$
, which satisfies this condition.
Lastly, the time delay contained in our system acts at all the time scales we are considering (Das & Chatterjee Reference Das and Chatterjee2002). Considering the
$\unicode[STIX]{x1D716}$
scalings just discussed for the slow time scales, delayed variables are therefore functions of
$t_{0}-\unicode[STIX]{x1D70F}$
,
$t_{2}-\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D70F}$
, and
$t_{4}-\unicode[STIX]{x1D716}^{4}\unicode[STIX]{x1D70F}$
. These terms are then expanded in series of
$\unicode[STIX]{x1D716}$
. For ease of notation, in the following we will adopt the short notation
$\boldsymbol{x}(t)$
for
$\boldsymbol{x}(t_{0},t_{2},t_{4})$
, and
$\boldsymbol{x}(t-\unicode[STIX]{x1D70F})$
for
$\boldsymbol{x}(t_{0}-\unicode[STIX]{x1D70F},t_{2},t_{4})$
.
We now substitute the relations (4.1), (4.4) into equation (4.3). The complete list of terms we obtain is given in appendix B. By matching these terms by their
$\unicode[STIX]{x1D716}$
order, we obtain a set of linear, inhomogeneous differential equations which have to be solved in ascending order. We perform the weakly nonlinear expansion up to
$O(\unicode[STIX]{x1D716}^{5})$
; in the following subsections we will solve and discuss the equations order by order.
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}$
:
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{1}}{\unicode[STIX]{x2202}t_{0}}-\unicode[STIX]{x1D63C}\boldsymbol{x}_{1}-\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))=\mathbf{0}.\end{eqnarray}$$
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:
$$\begin{eqnarray}\unicode[STIX]{x1D648}_{s}\equiv \left(s\unicode[STIX]{x1D644}-\unicode[STIX]{x1D63C}-\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{B}\boldsymbol{C}\text{e}^{-s\unicode[STIX]{x1D70F}}\right)\!,\end{eqnarray}$$
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:
$$\begin{eqnarray}\boldsymbol{x}_{1}\approx W(t_{2},t_{4})\boldsymbol{x}_{1}^{W}\text{e}^{\text{i}\unicode[STIX]{x1D714}^{c}t_{0}}+\text{c.c.},\end{eqnarray}$$
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:
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{2}}{\unicode[STIX]{x2202}t_{0}}-\unicode[STIX]{x1D63C}\boldsymbol{x}_{2}-\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{2}(t-\unicode[STIX]{x1D70F}))=|W|^{2}\boldsymbol{F}_{2}^{|W|^{2}}+\left(W^{2}\boldsymbol{F}_{2}^{W^{2}}\text{e}^{2\text{i}\unicode[STIX]{x1D714}^{c}t_{0}}+\text{c.c.}\right)\!,\end{eqnarray}$$
where
$$\begin{eqnarray}\left.\begin{array}{@{}l@{}}\boldsymbol{F}_{2}^{W^{2}}\equiv \unicode[STIX]{x1D6FC}_{2}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{2}\text{e}^{-2\text{i}\unicode[STIX]{x1D714}^{c}\unicode[STIX]{x1D70F}},\\ \boldsymbol{F}_{2}^{|W|^{2}}\equiv 2\unicode[STIX]{x1D6FC}_{2}K^{c}\boldsymbol{B}|\boldsymbol{C}\boldsymbol{x}_{1}^{W}|^{2}.\end{array}\right\}\end{eqnarray}$$
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
$$\begin{eqnarray}\boldsymbol{x}_{2}=|W|^{2}\boldsymbol{x}_{2}^{|W|^{2}}+\left(W^{2}\boldsymbol{x}_{2}^{W^{2}}\text{e}^{2\text{i}\unicode[STIX]{x1D714}^{c}t_{0}}+\text{c.c.}\right)\!.\end{eqnarray}$$
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:
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D648}_{2\text{i}\unicode[STIX]{x1D714}^{c}}\boldsymbol{x}_{2}^{W^{2}}=\boldsymbol{F}_{2}^{W^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D648}_{0}\boldsymbol{x}_{2}^{|W|^{2}}=\boldsymbol{F}_{2}^{|W|^{2}}. & \displaystyle\end{eqnarray}$$
$\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}$
:
$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{x}_{3}}{\unicode[STIX]{x2202}t_{0}}-\unicode[STIX]{x1D63C}\boldsymbol{x}_{3}-\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{3}(t-\unicode[STIX]{x1D70F}))\nonumber\\ \displaystyle & & \displaystyle \quad =-\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{1}}{\unicode[STIX]{x2202}t_{2}}-\unicode[STIX]{x1D70F}\boldsymbol{B}\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{C}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{1}}{\unicode[STIX]{x2202}t_{2}}(t-\unicode[STIX]{x1D70F})\cdots \nonumber\\ \displaystyle & & \displaystyle \qquad +\,\left(W\boldsymbol{F}_{3}^{W}\text{e}^{\text{i}\unicode[STIX]{x1D714}^{c}t_{0}}+|W|^{2}W\boldsymbol{F}_{3}^{|W|^{2}W}\text{e}^{\text{i}\unicode[STIX]{x1D714}^{c}t_{0}}+W^{3}\boldsymbol{F}_{3}^{W^{3}}\text{e}^{3\text{i}\unicode[STIX]{x1D714}^{c}t_{0}}+\text{c.c.}\right)\!.\end{eqnarray}$$
The explicit expressions of the forcing terms are reported in § C.1. The derivatives of the slow time scales may be rewritten as:
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{1}}{\unicode[STIX]{x2202}t_{2}}(t)=\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}t_{2}}\boldsymbol{x}_{1}^{W}\text{e}^{\text{i}\unicode[STIX]{x1D714}^{c}t}+\text{c.c.}\end{eqnarray}$$
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:
$$\begin{eqnarray}\langle \boldsymbol{y},\unicode[STIX]{x1D648}_{\,\text{i}\unicode[STIX]{x1D714}^{c}}^{}\boldsymbol{x}\rangle =\langle \unicode[STIX]{x1D648}_{\,\text{i}\unicode[STIX]{x1D714}^{c}}^{\dagger }\boldsymbol{y},\boldsymbol{x}\rangle\end{eqnarray}$$
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:
$$\begin{eqnarray}\left\langle \boldsymbol{x}_{1}^{\dagger },-\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}t_{2}}\unicode[STIX]{x1D64B}\boldsymbol{x}_{1}^{W}+W\boldsymbol{F}_{3}^{W}+|W|^{2}W\boldsymbol{F}_{3}^{|W|^{2}W}\right\rangle =0,\end{eqnarray}$$
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:
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}t_{2}}=\unicode[STIX]{x1D706}_{3}W+\unicode[STIX]{x1D708}_{3}|W|^{2}W,\end{eqnarray}$$
where the complex values
$\unicode[STIX]{x1D706}_{3}$
,
$\unicode[STIX]{x1D708}_{3}$
, known as the Landau coefficients, are defined by:
$$\begin{eqnarray}\unicode[STIX]{x1D706}_{3}\equiv \frac{\langle \boldsymbol{x}_{1}^{\dagger },\boldsymbol{F}_{3}^{W}\rangle }{\left\langle \boldsymbol{x}_{1}^{\dagger },\unicode[STIX]{x1D64B}\boldsymbol{x}_{1}^{W}\right\rangle },\quad \unicode[STIX]{x1D708}_{3}\equiv \frac{\left\langle \boldsymbol{x}_{1}^{\dagger },\boldsymbol{F}_{3}^{|W|^{2}W}\right\rangle }{\left\langle \boldsymbol{x}_{1}^{\dagger },\unicode[STIX]{x1D64B}\boldsymbol{x}_{1}^{W}\right\rangle }.\end{eqnarray}$$
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:
$$\begin{eqnarray}\boldsymbol{x}_{3}=W\boldsymbol{x}_{3}^{W}\text{e}^{\text{i}\unicode[STIX]{x1D714}^{c}t_{0}}+|W|^{2}W\boldsymbol{x}_{3}^{|W|^{2}W}\text{e}^{\text{i}\unicode[STIX]{x1D714}^{c}t_{0}}+W^{3}\boldsymbol{x}_{3}^{W^{3}}\text{e}^{3\text{i}\unicode[STIX]{x1D714}^{c}t_{0}}+\text{c.c.}\end{eqnarray}$$
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:
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D648}_{\,\text{i}\unicode[STIX]{x1D714}^{c}}\boldsymbol{x}_{3}^{W}=\boldsymbol{F}_{3}^{W}-\unicode[STIX]{x1D706}_{3}\unicode[STIX]{x1D64B}\boldsymbol{x}_{1}^{W}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D648}_{\,\text{i}\unicode[STIX]{x1D714}^{c}}\boldsymbol{x}_{3}^{|W|^{2}W}=\boldsymbol{F}_{3}^{|W|^{2}W}-\unicode[STIX]{x1D708}_{3}\unicode[STIX]{x1D64B}\boldsymbol{x}_{1}^{W}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D648}_{3\text{i}\unicode[STIX]{x1D714}^{c}}\boldsymbol{x}_{3}^{W^{3}}=\boldsymbol{F}_{3}^{W^{3}}. & \displaystyle\end{eqnarray}$$
$\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.3.1 Stuart–Landau equation:
$O(\unicode[STIX]{x1D716}^{3})$
Equation (4.17) is known as the Stuart–Landau equation. Its roots yield the amplitude of limit cycle solutions and the frequency shift of the nonlinear oscillation with respect to the marginally stable eigenfrequency. By using the polar representation
$W=r\text{e}^{\text{i}\unicode[STIX]{x1D703}}$
, and by splitting the real and the imaginary parts, we have:
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}t_{2}}=\text{Re}(\unicode[STIX]{x1D706}_{3})r+\text{Re}(\unicode[STIX]{x1D708}_{3})r^{3}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}t_{2}}=\text{Im}(\unicode[STIX]{x1D706}_{3})+\text{Im}(\unicode[STIX]{x1D708}_{3})r^{2}. & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D703}$
is valid only for solutions with non-zero amplitude. Steady-state solutions are reached when the amplitude
$r$
does not vary in time. The amplitude levels at which this happens are: 
$$\begin{eqnarray}r_{1}=0,\quad r_{2}=\sqrt{-\frac{\text{Re}(\unicode[STIX]{x1D706}_{3})}{\text{Re}(\unicode[STIX]{x1D708}_{3})}}.\end{eqnarray}$$
Figure 6. (a) Saturation of the Landau coefficients
$\unicode[STIX]{x1D706}_{3}$
,
$\unicode[STIX]{x1D708}_{3}$
with respect to the number of modes. (b) Subcritical Hopf bifurcation diagram. Solid and dashed lines represent stable and unstable solutions, respectively. The oscillation’s amplitude at the fundamental frequency with corrections up to
$\unicode[STIX]{x1D716}^{3}$
is shown.
From the definitions (4.18), (C 2a
), one can see that
$\unicode[STIX]{x1D706}_{3}$
, being proportional to
$\unicode[STIX]{x1D6FF}_{2}=\pm 1$
, changes sign across the critical point
$K^{c}$
, whereas
$\unicode[STIX]{x1D708}_{3}$
does not vary when we change the bifurcation parameter. Therefore, the solution
$r_{2}$
is real only on one side of the Hopf location. The stability of the solutions is connected to the sign of the eigenvalue of the Jacobian
$J\equiv \text{Re}(\unicode[STIX]{x1D706}_{3})+3\,\text{Re}(\unicode[STIX]{x1D708}_{3})r^{2}$
evaluated at the two solutions. These values are:
$$\begin{eqnarray}J(r_{1})=\text{Re}(\unicode[STIX]{x1D706}_{3}),\quad J(r_{2})=-2\,\text{Re}(\unicode[STIX]{x1D706}_{3}).\end{eqnarray}$$
Note that, in the region where two solutions coexist, their stability is different. This distinguishes between super- and subcritical Hopf bifurcations. For the set of parameters we have chosen, we find that the bifurcation is subcritical: before the Hopf bifurcation, stable fixed points and unstable limit cycle solutions exist, and after it only unstable fixed point solutions are found, as shown in figure 6(b). The present analysis could be also applied to configurations in which the bifurcation is supercritical.
The solution of the phase equation (4.21b ) on the unstable limit cycle solution reads:
$$\begin{eqnarray}\unicode[STIX]{x1D703}=\unicode[STIX]{x1D716}^{2}\left(\text{Im}(\unicode[STIX]{x1D706}_{3})-\text{Im}(\unicode[STIX]{x1D708}_{3})\frac{\text{Re}(\unicode[STIX]{x1D706}_{3})}{\text{Re}(\unicode[STIX]{x1D708}_{3})}\right)t_{0}\equiv \unicode[STIX]{x1D6E5}\unicode[STIX]{x1D714}t_{0},\end{eqnarray}$$
where we have used the scaling
$t_{2}=\unicode[STIX]{x1D716}^{2}t_{0}$
between the fast and slow time scales.
$\unicode[STIX]{x1D6E5}\unicode[STIX]{x1D714}$
represents the frequency shift between the fundamental oscillation frequency of limit cycles and the marginally stable frequency
$\unicode[STIX]{x1D714}^{c}$
.
Combining the power expansion (4.1), the weakly nonlinear solutions (4.8), (4.11), (4.19) and the solution of the Stuart–Landau equation, we obtain an analytical expression for the time evolution of the thermoacoustic states up to third order, which reads:
$$\begin{eqnarray}\displaystyle \boldsymbol{x} & = & \displaystyle \unicode[STIX]{x1D716}^{2}r^{2}\boldsymbol{x}_{2}^{|W|^{2}}+\left[\unicode[STIX]{x1D716}r\boldsymbol{x}_{1}^{W}+\unicode[STIX]{x1D716}^{3}r\boldsymbol{x}_{3}^{W}+\unicode[STIX]{x1D716}^{3}r^{3}\boldsymbol{x}_{3}^{|W|^{2}W}\right]\text{e}^{\text{i}(\unicode[STIX]{x1D714}^{c}+\unicode[STIX]{x1D6E5}\unicode[STIX]{x1D714})t_{0}}\cdots \nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D716}^{2}r^{2}\boldsymbol{x}_{2}^{W^{2}}\text{e}^{2\text{i}(\unicode[STIX]{x1D714}^{c}+\unicode[STIX]{x1D6E5}\unicode[STIX]{x1D714})t_{0}}+\unicode[STIX]{x1D716}^{3}r^{3}\boldsymbol{x}_{3}^{W^{3}}\text{e}^{3\text{i}(\unicode[STIX]{x1D714}^{c}+\unicode[STIX]{x1D6E5}\unicode[STIX]{x1D714})t_{0}}+\text{c.c.}+O(\unicode[STIX]{x1D716}^{4}).\end{eqnarray}$$
Figure 6(b) shows the subcritical bifurcation diagram we obtain at this order. It contains the amplitude level of the velocity fluctuations
$u^{\prime }=\boldsymbol{C}\boldsymbol{x}$
calculated from (4.25) at the fundamental frequency. We shall postpone the discussion of the response at other frequencies to § 5. Because these limit cycle oscillations are unstable, the amplitudes that the Stuart–Landau equation predicts at this order correspond approximately to the level of triggering which is required to excite finite-amplitude oscillations. In order to predict the amplitude of stable limit cycles, we need to extend the weakly nonlinear expansion to higher orders.
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}$
:
$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{x}_{4}}{\unicode[STIX]{x2202}t_{0}}-\unicode[STIX]{x1D63C}\boldsymbol{x}_{4}-\unicode[STIX]{x1D6FC}_{1}K\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{4}(t-\unicode[STIX]{x1D70F}))\nonumber\\ \displaystyle & & \displaystyle \quad =|W|^{4}\boldsymbol{F}_{4}^{|W|^{4}}+|W|^{2}\boldsymbol{F}_{4}^{|W|^{2}}\cdots \nonumber\\ \displaystyle & & \displaystyle \qquad +\,\left(W^{2}\boldsymbol{F}_{4}^{W^{2}}\text{e}^{2\text{i}\unicode[STIX]{x1D714}^{c}t_{0}}+|W|^{2}W^{2}\boldsymbol{F}_{4}^{|W|^{2}W^{2}}\text{e}^{2\text{i}\unicode[STIX]{x1D714}^{c}t_{0}}+W^{4}\boldsymbol{F}_{4}^{W^{4}}\text{e}^{4i\unicode[STIX]{x1D714}^{c}t_{0}}+\text{c.c.}\right).\end{eqnarray}$$
None of the forcing terms resonate. Their expressions are provided in § C.2. Using the ansatz:
$$\begin{eqnarray}\displaystyle \boldsymbol{x}_{4} & = & \displaystyle |W|^{2}\boldsymbol{x}_{4}^{|W|^{2}}+|W|^{4}\boldsymbol{x}_{4}^{|W|^{4}}\cdots \nonumber\\ \displaystyle & & \displaystyle +\,\left(W^{2}\boldsymbol{x}_{4}^{W^{2}}\text{e}^{2\text{i}\unicode[STIX]{x1D714}^{c}t_{0}}+|W|^{2}W^{2}\boldsymbol{x}_{4}^{|W|^{2}W^{2}}\text{e}^{2\text{i}\unicode[STIX]{x1D714}^{c}t_{0}}+W^{4}\boldsymbol{x}_{4}^{W^{4}}\text{e}^{4\text{i}\unicode[STIX]{x1D714}^{c}t_{0}}+\text{c.c.}\right)\!,\end{eqnarray}$$
we can readily calculate the solutions
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{x}_{4}^{|W|^{2}}=\unicode[STIX]{x1D648}_{0}^{-1}\boldsymbol{F}_{4}^{|W|^{2}} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{x}_{4}^{|W|^{4}}=\unicode[STIX]{x1D648}_{0}^{-1}\boldsymbol{F}_{4}^{|W|^{4}} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{x}_{4}^{W^{2}}=\unicode[STIX]{x1D648}_{2\text{i}\unicode[STIX]{x1D714}^{c}}^{-1}\boldsymbol{F}_{4}^{W^{2}} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{x}_{4}^{|W|^{2}W^{2}}=\unicode[STIX]{x1D648}_{2\text{i}\unicode[STIX]{x1D714}^{c}}^{-1}\boldsymbol{F}_{4}^{|W|^{2}W^{2}} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{x}_{4}^{W^{4}}=\unicode[STIX]{x1D648}_{4\text{i}\unicode[STIX]{x1D714}^{c}}^{-1}\boldsymbol{F}_{4}^{W^{4}}, & \displaystyle\end{eqnarray}$$
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}$
.
$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{x}_{5}}{\unicode[STIX]{x2202}t_{0}}-\unicode[STIX]{x1D63C}\boldsymbol{x}_{5}-\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{5})\nonumber\\ \displaystyle & & \displaystyle \quad =-\unicode[STIX]{x1D64B}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{1}}{\unicode[STIX]{x2202}t_{4}}\cdots +\left(W\boldsymbol{F}_{5}^{W}+|W|^{2}W\boldsymbol{F}_{5}^{|W|^{2}W}+|W|^{4}W\boldsymbol{F}_{5}^{|W|^{4}W}\right)\text{e}^{\text{i}\unicode[STIX]{x1D714}^{c}(t_{0}-\unicode[STIX]{x1D70F})}\cdots \nonumber\\ \displaystyle & & \displaystyle \qquad +\,\left(|W|^{2}W^{3}\boldsymbol{F}_{5}^{|W|^{2}W^{3}}+W^{3}\boldsymbol{F}_{5}^{W^{3}}\right)\text{e}^{3\text{i}\unicode[STIX]{x1D714}^{c}(t_{0}-\unicode[STIX]{x1D70F})}+W^{5}\boldsymbol{F}_{5}^{W^{5}}\text{e}^{5\text{i}\unicode[STIX]{x1D714}^{c}(t_{0}-\unicode[STIX]{x1D70F})}+\text{c.c.}\end{eqnarray}$$
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:
$$\begin{eqnarray}\left\langle \boldsymbol{x}_{1}^{\dagger },-\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}t_{4}}\unicode[STIX]{x1D64B}\boldsymbol{x}_{1}^{W}+W\boldsymbol{F}_{5}^{W}+|W|^{2}W\boldsymbol{F}_{5}^{|W|^{2}W}+|W|^{4}W\boldsymbol{F}_{5}^{|W|^{4}W}\right\rangle =0.\end{eqnarray}$$
The resonant forcing expressions can be found in § C.3. Equation (4.30) can be simplified into the Stuart–Landau equation:
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}t_{4}}=\unicode[STIX]{x1D706}_{5}W+\unicode[STIX]{x1D708}_{5}|W|^{2}W+\unicode[STIX]{x1D707}_{5}|W|^{4}W,\end{eqnarray}$$
where the Landau coefficients are defined by:
$$\begin{eqnarray}\unicode[STIX]{x1D706}_{5}\equiv \frac{\left\langle \boldsymbol{x}_{1}^{\dagger },\boldsymbol{F}_{5}^{W}\right\rangle }{\left\langle \boldsymbol{x}_{1}^{\dagger },\unicode[STIX]{x1D64B}\boldsymbol{x}_{1}^{W}\right\rangle },\quad \unicode[STIX]{x1D708}_{5}\equiv \frac{\left\langle \boldsymbol{x}_{1}^{\dagger },\boldsymbol{F}_{5}^{|W|^{2}W}\right\rangle }{\left\langle \boldsymbol{x}_{1}^{\dagger },\unicode[STIX]{x1D64B}\boldsymbol{x}_{1}^{W}\right\rangle },\quad \unicode[STIX]{x1D707}_{5}\equiv \frac{\left\langle \boldsymbol{x}_{1}^{\dagger },\boldsymbol{F}_{5}^{|W|^{4}W}\right\rangle }{\left\langle \boldsymbol{x}_{1}^{\dagger },\unicode[STIX]{x1D64B}\boldsymbol{x}_{1}^{W}\right\rangle }.\end{eqnarray}$$
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.
4.5.1 Stuart–Landau equation:
$O(\unicode[STIX]{x1D716}^{5})$
The overall slow time scale evolution of the amplitude
$W$
is obtained by combining the results at the two time scales
$t_{2}$
and
$t_{4}$
(Fujimura Reference Fujimura1991; Gambino, Lombardo & Sammartino Reference Gambino, Lombardo and Sammartino2012). By using the scaling
$t_{4}=\unicode[STIX]{x1D716}^{2}t_{2}$
, we obtain:
$$\begin{eqnarray}\frac{\text{d}W}{\text{d}t_{2}}=\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}t_{2}}+\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}t_{4}}\frac{\unicode[STIX]{x2202}t_{4}}{\unicode[STIX]{x2202}t_{2}}=\left(\unicode[STIX]{x1D706}_{3}+\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D706}_{5}\right)W+\left(\unicode[STIX]{x1D708}_{3}+\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D708}_{5}\right)|W|^{2}W+\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D707}_{5}|W|^{4}W.\end{eqnarray}$$
By using the polar representation
$W=r\text{e}^{\text{i}\unicode[STIX]{x1D703}}$
, this decouples into:
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}r}{\text{d}t_{2}}=\text{Re}(\unicode[STIX]{x1D706})r+\text{Re}(\unicode[STIX]{x1D708})r^{3}+\text{Re}(\unicode[STIX]{x1D707})r^{5}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D703}}{\text{d}t_{2}}=\text{Im}(\unicode[STIX]{x1D706})+\text{Im}(\unicode[STIX]{x1D708})r^{2}+\text{Im}(\unicode[STIX]{x1D707})r^{4}, & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D706}=\unicode[STIX]{x1D706}_{3}+\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D706}_{5}$
, and similarly for
$\unicode[STIX]{x1D708}$
and
$\unicode[STIX]{x1D707}$
. The fixed points of the amplitude equation are: 
$$\begin{eqnarray}r_{1}=0,\quad r_{2,3}=\sqrt{\frac{-\text{Re}(\unicode[STIX]{x1D708})\pm \sqrt{(\text{Re}(\unicode[STIX]{x1D708}))^{2}-4\,\text{Re}(\unicode[STIX]{x1D706})\,\text{Re}(\unicode[STIX]{x1D707})}}{2\,\text{Re}(\unicode[STIX]{x1D707})}}.\end{eqnarray}$$
Figure 7. Comparison between third-order (black) and fifth-order (red) bifurcation diagrams. The analytically calculated amplitude of the acoustic velocity at the fundamental frequency is shown, according to (4.25). At fifth order the limit cycle saturates, so the limit cycle amplitude and the location of the fold point
$K^{F}$
can be calculated.
The existence of real solutions for
$r_{2,3}$
depends only on the sign of the terms under the square roots. Because the bifurcation is subcritical, we have two solutions after the Hopf location (an unstable fixed point and a stable limit cycle), a bistable region between the Hopf and the fold points with three solutions (a stable fixed point, an unstable and a stable limit cycle), and a region with only one stable fixed point before the fold. This is shown in figure 7, where the path that the oscillations will follow when the bifurcation parameter is varied across the bistable region is shown with arrows. As expected, close to the Hopf point the fifth-order expansion correctly matches the third-order analysis.
The stability of the solutions is determined by the sign of the Jacobian
$J=\text{Re}(\unicode[STIX]{x1D706})+3\,\text{Re}(\unicode[STIX]{x1D708})r^{2}+5\,\text{Re}(\unicode[STIX]{x1D707})r^{4}$
evaluated at the solutions. These values are:
$$\begin{eqnarray}\displaystyle & \displaystyle J(r_{1})=\text{Re}(\unicode[STIX]{x1D706}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle J(r_{2,3})=-4\,\text{Re}(\unicode[STIX]{x1D706})+\frac{(\text{Re}(\unicode[STIX]{x1D708}))^{2}}{\text{Re}(\unicode[STIX]{x1D707})}\mp \text{Re}(\unicode[STIX]{x1D708})\sqrt{(\text{Re}(\unicode[STIX]{x1D708}))^{2}-4\,\text{Re}(\unicode[STIX]{x1D706})\text{Re}(\unicode[STIX]{x1D707})}. & \displaystyle\end{eqnarray}$$
The frequency shift
$\unicode[STIX]{x1D6E5}\unicode[STIX]{x1D714}$
on the limit cycles can be readily calculated from (4.34b
) by using the scaling
$t_{2}=\unicode[STIX]{x1D716}^{2}t_{0}$
:
$$\begin{eqnarray}\unicode[STIX]{x1D6E5}\unicode[STIX]{x1D714}_{2,3}\equiv \unicode[STIX]{x1D716}^{2}\left(\text{Im}(\unicode[STIX]{x1D706})+\text{Im}(\unicode[STIX]{x1D708})r_{2,3}^{2}+\text{Im}(\unicode[STIX]{x1D707})r_{2,3}^{4}\right)\!.\end{eqnarray}$$
5 Results validation
We now compare the weakly nonlinear analysis discussed in § 4 with the fully nonlinear results obtained by solving the nonlinear dynamical system (3.1) with no approximations. This is achieved with two methods: time marching the governing equations and numerical continuation of limit cycles. Time domain simulations are performed using MATLAB delay differential equations solver dde23 and the numerical continuation algorithm is based on the DDE-BIFTOOL package (Engelborghs et al. Reference Engelborghs, Luzyanina and Roose2002). Although the latter can only predict periodic oscillations and their stability, we find that, for the set of parameters we have investigated, the system always converges towards fixed points or periodic oscillations. Therefore, the two methods yield the same results for the steady-state response, as was verified.
Figure 8. (a) At
$K=K^{c}+\unicode[STIX]{x0394}K=1.43$
the fixed point solution becomes unstable and oscillations grow to a limit cycle. (b) At
$K=1.421$
no limit cycle solutions exist, and the system converges to fixed point solutions even when initialised to a highly perturbed state. The insets show the exponential growth and decay rates of the oscillations on a logarithmic scale. The growth rate is extracted by performing a fit (red line in the insets) onto the magnitude of the Hilbert transformed signal.
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.
6 Conclusions
In this study we have investigated nonlinear thermoacoustic oscillations in a Rijke tube using a fifth-order weakly nonlinear expansion of the governing equations. The framework we analysed, composed of a wave-based approach for the linear acoustic equations and a nonlinear model for the heat release response, is general and can be extended to more complex networks. We have shown how a weakly nonlinear expansion close to the Hopf bifurcations can be used to predict the amplitude and frequency of stable and unstable limit cycles and the region of bistability for subcritical bifurcations. We have compared our weakly nonlinear results with those obtained by time marching the fully nonlinear thermoacoustic equations, which does not introduce any further approximation in the governing equations, and we have shown that the method yields accurate results when the expansion is truncated at fifth order. Using this type of analysis, the numerical/experimental effort needed to construct a bifurcation diagram can be significantly reduced.
The method requires an analytical model for the description of the nonlinear heat release dynamics. It can be challenging, however, to obtain such a model when, for example, the dynamics of a turbulent flame is considered. To overcome this issue, nonlinear system identification procedures could be used to fit the
$\unicode[STIX]{x1D6FC}$
coefficients of the weakly nonlinear response onto the flame nonlinear response (measured with a series of experiments). The latter procedure can be performed only around frequencies at which thermoacoustic instabilities are expected. The latter can be easily determined via a network model if flame transfer function (FTF) measurements are available. All the information required for the weakly nonlinear analysis is then known and the theory can be applied.
The weakly nonlinear analysis we performed is valid for a deterministic system, close to a Hopf bifurcation, and assumes that only one thermoacoustic mode is marginally stable and all the others are damped. The influence of turbulence and combustion noise can be modelled phenomenologically as random forcing (Noiray & Schuermans Reference Noiray and Schuermans2013; Rigas et al. Reference Rigas, Morgans, Brackston and Morrison2015b ), the effects of which can be captured by including stochastic forcing in the Stuart–Landau equation. In certain problems, a second mode may become marginally stable in the vicinity of the first Hopf location. In this case, the theory can be extended by considering a codimension two bifurcation, as described in Meliga et al. (Reference Meliga, Chomaz and Sipp2009). It can also be extended to investigate the response of the system to external unsteady or steady forcing. In this case, the forcing terms appear explicitly in the amplitude equations, as shown by Sipp (Reference Sipp2012), and the weakly nonlinear interaction between the global mode and the forcing is taken explicitly into account. Moreover, the unknown Landau coefficients of the model can be identified from transient and steady-state measurements when the system is subjected to external forcing (Rigas, Morgans & Morrison Reference Rigas, Morgans and Morrison2016) even for turbulent cases. This method could provide a systematic way for devising open-loop control strategies for the regulation of thermoacoustic oscillations (Paschereit & Gutmark Reference Paschereit and Gutmark2008; Cósić et al. Reference Cósić, Bobusch, Moeck and Paschereit2012).
Acknowledgements
This study was founded by the European Research Council through project ALORS 2590620.
Appendix A. Acoustic matrix coefficients
The coefficients of the matrix in (2.1) are given by:
$$\begin{eqnarray}\displaystyle \!\!\! & \left.\begin{array}{@{}rcl@{}}M_{11}\ & =\ & \displaystyle \unicode[STIX]{x0394}M_{2}+M_{1}(M_{1}-\unicode[STIX]{x0394}M_{2}-2)+\unicode[STIX]{x1D6E9}-\cdots \\ \ & \ & \displaystyle \text{e}^{-s\unicode[STIX]{x1D70F}_{1}}R_{1}(-\unicode[STIX]{x0394}M_{2}+M1(M_{1}-\unicode[STIX]{x0394}M_{2}+2)+\unicode[STIX]{x1D6E9})\\ M_{12}\ & =\ & \displaystyle \unicode[STIX]{x1D6E9}\left[1+M_{2}-\text{e}^{-s\unicode[STIX]{x1D70F}_{2}}R_{2}(M_{2}-1)\right]\\ M_{21}\ & =\ & \displaystyle \frac{1}{2(\unicode[STIX]{x1D6FE}-1)}\left[(M_{1}-1)\left(2+M_{1}(M_{1}-2)(\unicode[STIX]{x1D6FE}-1)+\unicode[STIX]{x0394}^{2}M_{2}^{2}(1-\unicode[STIX]{x1D6FE})\right)-\right.\cdots \\ \ & \ & \displaystyle \left.\text{e}^{-s\unicode[STIX]{x1D70F}_{1}}R_{1}\left((M_{1}+1)(2+M_{1}(M_{1}+2)(\unicode[STIX]{x1D6FE}-1)+\unicode[STIX]{x0394}^{2}M_{2}^{2}(1-\unicode[STIX]{x1D6FE}))\right)\right]\\ M_{22}\ & =\ & \displaystyle \frac{\unicode[STIX]{x1D6E5}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x1D6FE}-1}\left[(M_{2}+1)(1+M_{2}(\unicode[STIX]{x1D6FE}-1))+\text{e}^{-s\unicode[STIX]{x1D70F}_{2}}R_{2}(M_{2}^{2}-1)(1-\unicode[STIX]{x1D6FE})\right],\end{array}\right\}\qquad & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FE}$
is the specific heat ratio of air and
$M_{1,2}\equiv \overline{u}_{1,2}/c_{1,2}$
is the Mach number.
Appendix B. Nonlinear expansion terms
In the following all the terms obtained by expanding (3.1) to fifth order in
$\unicode[STIX]{x1D716}$
are listed and classified by their physical origin.
-
(i) Fast time scale (
$\unicode[STIX]{x2202}\boldsymbol{x}/\unicode[STIX]{x2202}t_{0}$
): (B 1)
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{1}}{\unicode[STIX]{x2202}t_{0}}+\unicode[STIX]{x1D716}^{2}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{2}}{\unicode[STIX]{x2202}t_{0}}+\unicode[STIX]{x1D716}^{3}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{3}}{\unicode[STIX]{x2202}t_{0}}+\unicode[STIX]{x1D716}^{4}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{4}}{\unicode[STIX]{x2202}t_{0}}+\unicode[STIX]{x1D716}^{5}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{5}}{\unicode[STIX]{x2202}t_{0}}. & \displaystyle\end{eqnarray}$$
-
(ii) First slow time scale (
$\unicode[STIX]{x2202}\boldsymbol{x}/\unicode[STIX]{x2202}t_{2}$
): (B 2)
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}^{3}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{1}}{\unicode[STIX]{x2202}t_{2}}+\unicode[STIX]{x1D716}^{4}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{2}}{\unicode[STIX]{x2202}t_{2}}+\unicode[STIX]{x1D716}^{5}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{3}}{\unicode[STIX]{x2202}t_{2}}. & \displaystyle\end{eqnarray}$$
-
(iii) Second slow time scale (
$\unicode[STIX]{x2202}\boldsymbol{x}/\unicode[STIX]{x2202}t_{4}$
): (B 3)
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}^{5}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{1}}{\unicode[STIX]{x2202}t_{4}}. & \displaystyle\end{eqnarray}$$
-
(iv) Linear acoustics (
$\unicode[STIX]{x1D63C}\boldsymbol{x}$
): (B 4)
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}\unicode[STIX]{x1D63C}\boldsymbol{x}_{1}+\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D63C}\boldsymbol{x}_{2}+\unicode[STIX]{x1D716}^{3}\unicode[STIX]{x1D63C}\boldsymbol{x}_{3}+\unicode[STIX]{x1D716}^{4}\unicode[STIX]{x1D63C}\boldsymbol{x}_{4}+\unicode[STIX]{x1D716}^{5}\unicode[STIX]{x1D63C}\boldsymbol{x}_{5}. & \displaystyle\end{eqnarray}$$
-
(v) Linear heat release (
$\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}(t-\unicode[STIX]{x1D70F}))$
): (B 5)
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D716}\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))+\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{2}(t-\unicode[STIX]{x1D70F}))+\unicode[STIX]{x1D716}^{3}\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{3}(t-\unicode[STIX]{x1D70F}))\cdots \nonumber\\ \displaystyle & & \displaystyle \quad +\,\unicode[STIX]{x1D716}^{4}\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{4}(t-\unicode[STIX]{x1D70F}))+\unicode[STIX]{x1D716}^{5}\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{5}(t-\unicode[STIX]{x1D70F})).\end{eqnarray}$$
-
(vi) Quadratic heat release (
$\unicode[STIX]{x1D6FC}_{2}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}(t-\unicode[STIX]{x1D70F}))^{2}$
): (B 6)
$$\begin{eqnarray}\displaystyle & & \displaystyle \!\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D6FC}_{2}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))^{2}+\unicode[STIX]{x1D716}^{3}2\unicode[STIX]{x1D6FC}_{2}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))(\boldsymbol{C}\boldsymbol{x}_{2}(t-\unicode[STIX]{x1D70F}))\cdots \nonumber\\ \displaystyle & & \displaystyle \!\quad +\,\unicode[STIX]{x1D716}^{4}2\unicode[STIX]{x1D6FC}_{2}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))(\boldsymbol{C}\boldsymbol{x}_{3}(t-\unicode[STIX]{x1D70F}))+\unicode[STIX]{x1D716}^{5}2\unicode[STIX]{x1D6FC}_{2}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))(\boldsymbol{C}\boldsymbol{x}_{4}(t-\unicode[STIX]{x1D70F}))\cdots \nonumber\\ \displaystyle & & \displaystyle \!\quad +\,\unicode[STIX]{x1D716}^{4}\unicode[STIX]{x1D6FC}_{2}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{2}(t-\unicode[STIX]{x1D70F}))^{2}+\unicode[STIX]{x1D716}^{5}2\unicode[STIX]{x1D6FC}_{2}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{2}(t-\unicode[STIX]{x1D70F}))(\boldsymbol{C}\boldsymbol{x}_{3}(t-\unicode[STIX]{x1D70F})).\quad\end{eqnarray}$$
-
(vii) Linear heat release and
$\unicode[STIX]{x0394}K$
coupling (
$\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D6FF}_{2}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}(t-\unicode[STIX]{x1D70F}))$
): (B 7)
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}^{3}\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D6FF}_{2}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))+\unicode[STIX]{x1D716}^{4}\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D6FF}_{2}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{2}(t-\unicode[STIX]{x1D70F}))+\unicode[STIX]{x1D716}^{5}\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D6FF}_{2}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{3}(t-\unicode[STIX]{x1D70F})). & \displaystyle\end{eqnarray}$$
-
(viii) Cubic heat release (
$\unicode[STIX]{x1D6FC}_{3}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}(t-\unicode[STIX]{x1D70F}))^{3}$
): (B 8)
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D716}^{3}\unicode[STIX]{x1D6FC}_{3}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))^{3}+\unicode[STIX]{x1D716}^{4}3\unicode[STIX]{x1D6FC}_{3}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))^{2}(\boldsymbol{C}\boldsymbol{x}_{2}(t-\unicode[STIX]{x1D70F}))\cdots \nonumber\\ \displaystyle & & \displaystyle \quad +\,\unicode[STIX]{x1D716}^{5}3\unicode[STIX]{x1D6FC}_{3}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))^{2}(\boldsymbol{C}\boldsymbol{x}_{3}(t-\unicode[STIX]{x1D70F}))+\unicode[STIX]{x1D716}^{5}3\unicode[STIX]{x1D6FC}_{3}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))(\boldsymbol{C}\boldsymbol{x}_{2}(t-\unicode[STIX]{x1D70F}))^{2}.\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
-
(ix) Quadratic heat release and
$\unicode[STIX]{x0394}K$
coupling (
$\unicode[STIX]{x1D6FC}_{2}\unicode[STIX]{x1D6FF}_{2}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}(t-\unicode[STIX]{x1D70F}))^{2}$
): (B 9)
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}^{4}\unicode[STIX]{x1D6FC}_{2}\unicode[STIX]{x1D6FF}_{2}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))^{2}+\unicode[STIX]{x1D716}^{5}2\unicode[STIX]{x1D6FC}_{2}\unicode[STIX]{x1D6FF}_{2}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))(\boldsymbol{C}\boldsymbol{x}_{2}(t-\unicode[STIX]{x1D70F})). & \displaystyle\end{eqnarray}$$
-
(x) Quartic heat release (
$\unicode[STIX]{x1D6FC}_{4}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}(t-\unicode[STIX]{x1D70F}))^{4}$
): (B 10)
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}^{4}\unicode[STIX]{x1D6FC}_{4}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))^{4}+\unicode[STIX]{x1D716}^{5}4\unicode[STIX]{x1D6FC}_{4}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))^{3}(\boldsymbol{C}\boldsymbol{x}_{2}(t-\unicode[STIX]{x1D70F})). & \displaystyle\end{eqnarray}$$
-
(xi) Cubic heat release and
$\unicode[STIX]{x0394}K$
coupling (
$\unicode[STIX]{x1D6FC}_{3}\unicode[STIX]{x1D6FF}_{2}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}(t-\unicode[STIX]{x1D70F}))^{3}$
): (B 11)
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}^{5}\unicode[STIX]{x1D6FC}_{3}\unicode[STIX]{x1D6FF}_{2}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))^{3}. & \displaystyle\end{eqnarray}$$
-
(xii) Quintic heat release (
$\unicode[STIX]{x1D6FC}_{5}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}(t-\unicode[STIX]{x1D70F}))^{5}$
): (B 12)
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}^{5}\unicode[STIX]{x1D6FC}_{5}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))^{5}. & \displaystyle\end{eqnarray}$$
-
(xiii) Time delay (proportional to
$\unicode[STIX]{x1D70F}$
,
$\unicode[STIX]{x1D70F}^{2}$
) (B 13)
$$\begin{eqnarray}\displaystyle & & \displaystyle -\unicode[STIX]{x1D70F}\boldsymbol{B}\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{C}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{1}}{\unicode[STIX]{x2202}t_{2}}(t-\unicode[STIX]{x1D70F})\unicode[STIX]{x1D716}^{3}\cdots \nonumber\\ \displaystyle & & \displaystyle \quad +\,\unicode[STIX]{x1D70F}\boldsymbol{B}\left(-2\unicode[STIX]{x1D6FC}_{2}K^{c}\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F})\boldsymbol{C}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{1}}{\unicode[STIX]{x2202}t_{2}}(t-\unicode[STIX]{x1D70F})-\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{C}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{2}}{\unicode[STIX]{x2202}t_{2}}(t-\unicode[STIX]{x1D70F})\right)\unicode[STIX]{x1D716}^{4}\cdots \nonumber\\ \displaystyle & & \displaystyle \quad +\,\unicode[STIX]{x1D70F}\boldsymbol{B}\left[-\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{C}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{1}}{\unicode[STIX]{x2202}t_{4}}(t-\unicode[STIX]{x1D70F})-\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D6FF}_{2}\boldsymbol{C}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{1}}{\unicode[STIX]{x2202}t_{2}}(t-\unicode[STIX]{x1D70F})\right.\cdots \nonumber\\ \displaystyle & & \displaystyle \quad -\left.3\unicode[STIX]{x1D6FC}_{3}K^{c}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))^{2}\boldsymbol{C}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{1}}{\unicode[STIX]{x2202}t_{2}}(t-\unicode[STIX]{x1D70F})-\,2\unicode[STIX]{x1D6FC}_{2}K^{c}\boldsymbol{C}\boldsymbol{x}_{2}(t-\unicode[STIX]{x1D70F})\boldsymbol{C}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{1}}{\unicode[STIX]{x2202}t_{2}}(t-\unicode[STIX]{x1D70F})\right.\cdots \nonumber\\ \displaystyle & & \displaystyle \quad -\left.2\unicode[STIX]{x1D6FC}_{2}K^{c}\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F})\boldsymbol{C}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{2}}{\unicode[STIX]{x2202}t_{2}}(t-\unicode[STIX]{x1D70F})-\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{C}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{3}}{\unicode[STIX]{x2202}t_{2}}(t-\unicode[STIX]{x1D70F})\cdots \,\right.\nonumber\\ \displaystyle & & \displaystyle \quad +\left.\frac{1}{2}\unicode[STIX]{x1D70F}\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{C}\frac{\unicode[STIX]{x2202}^{2}\boldsymbol{x}_{1}}{\unicode[STIX]{x2202}t_{2}^{2}}(t-\unicode[STIX]{x1D70F})\right]\unicode[STIX]{x1D716}^{5}.\end{eqnarray}$$
Appendix C. Forcing terms
C.1
$O(\unicode[STIX]{x1D716}^{3})$
The list of forcing terms is:
$$\begin{eqnarray}\displaystyle & \displaystyle 2\unicode[STIX]{x1D6FC}_{2}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))(\boldsymbol{C}\boldsymbol{x}_{2}(t-\unicode[STIX]{x1D70F}))+\unicode[STIX]{x1D6FC}_{3}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))^{3}+\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D6FF}_{2}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F})).\quad & \displaystyle\end{eqnarray}$$
Using the solutions (4.8), (4.11), we can expand these terms and classify them by their amplitude dependence as:
$$\begin{eqnarray}\boldsymbol{F}_{3}^{W}\equiv \unicode[STIX]{x1D6FC}_{1}\boldsymbol{B}\unicode[STIX]{x1D6FF}_{2}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{c}\unicode[STIX]{x1D70F}}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \boldsymbol{F}_{3}^{|W|^{2}W} & \equiv & \displaystyle \boldsymbol{B}K^{c}\big[3\unicode[STIX]{x1D6FC}_{3}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})|(\boldsymbol{C}\boldsymbol{x}_{1}^{W})|^{2}\cdots \nonumber\\ \displaystyle & & \displaystyle +\,2\unicode[STIX]{x1D6FC}_{2}\big((\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})(\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{\ast }+(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{2}^{|W|^{2}})\big)\big]\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{c}\unicode[STIX]{x1D70F}}\end{eqnarray}$$
$$\begin{eqnarray}\boldsymbol{F}_{3}^{W^{3}}\equiv \boldsymbol{B}K^{c}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})\big(\unicode[STIX]{x1D6FC}_{3}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{2}+2\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})\big)\text{e}^{-3\text{i}\unicode[STIX]{x1D714}^{c}\unicode[STIX]{x1D70F}}.\end{eqnarray}$$
C.2
$O(\unicode[STIX]{x1D716}^{4})$
The list of forcing terms is:
$$\begin{eqnarray}\displaystyle & & \displaystyle -\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{2}}{\unicode[STIX]{x2202}t_{2}}+\unicode[STIX]{x1D70F}\boldsymbol{B}\left(-2\unicode[STIX]{x1D6FC}_{2}K^{c}\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F})\boldsymbol{C}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{1}}{\unicode[STIX]{x2202}t_{2}}(t-\unicode[STIX]{x1D70F})-\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{C}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{2}}{\unicode[STIX]{x2202}t_{2}}(t-\unicode[STIX]{x1D70F})\right)\cdots \nonumber\\ \displaystyle & & \displaystyle \quad +\,2\unicode[STIX]{x1D6FC}_{2}K\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))(\boldsymbol{C}\boldsymbol{x}_{3}(t-\unicode[STIX]{x1D70F}))+\unicode[STIX]{x1D6FC}_{2}K\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{2}(t-\unicode[STIX]{x1D70F}))^{2}\cdots \nonumber\\ \displaystyle & & \displaystyle \quad +\,\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D6FF}_{2}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{2}(t-\unicode[STIX]{x1D70F}))+3\unicode[STIX]{x1D6FC}_{3}K\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))^{2}(\boldsymbol{C}\boldsymbol{x}_{2}(t-\unicode[STIX]{x1D70F}))\cdots \nonumber\\ \displaystyle & & \displaystyle \quad +\,\unicode[STIX]{x1D6FC}_{2}\unicode[STIX]{x1D6FF}_{2}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))^{2}+\unicode[STIX]{x1D6FC}_{4}K\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))^{4}.\end{eqnarray}$$
The definitions of the forcings in (4.26) read:
$$\begin{eqnarray}\displaystyle & & \displaystyle \boldsymbol{F}_{4}^{W^{4}}\equiv \boldsymbol{B}K\big(\unicode[STIX]{x1D6FC}_{4}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{4}+3\unicode[STIX]{x1D6FC}_{3}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{2}(\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})+2\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{3}^{W^{3}})+\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})^{2}\big)\qquad\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \boldsymbol{F}_{4}^{|W|^{4}} & \equiv & \displaystyle \boldsymbol{B}K^{c}\big(6\unicode[STIX]{x1D6FC}_{4}|\boldsymbol{C}\boldsymbol{x}_{1}^{W}|^{4}+(\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})^{\ast }\big(3\unicode[STIX]{x1D6FC}_{3}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{2}+2\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})\big)\cdots \nonumber\\ \displaystyle & & \displaystyle +\left.3\unicode[STIX]{x1D6FC}_{3}(\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})\left((\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{\ast }\right)^{2}+2(\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{\ast }(3\unicode[STIX]{x1D6FC}_{3}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{2}^{|W|^{2}})+\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{3}^{|W|^{2}W}))\cdots \,\right.\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}_{2}\big(2(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{3}^{|W|^{2}W})^{\ast }+(\boldsymbol{C}\boldsymbol{x}_{2}^{|W|^{2}})^{2}\big)\big)-2\boldsymbol{B}K^{c}\unicode[STIX]{x1D70F}\,\text{Re}(\unicode[STIX]{x1D708}_{3})\big(\unicode[STIX]{x1D6FC}_{1}(\boldsymbol{C}\boldsymbol{x}_{2}^{|W|^{2}})\cdots \nonumber\\ \displaystyle & & \displaystyle +\,2\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{\ast }\big)-2\boldsymbol{x}_{2}^{|W|^{2}}\,\text{Re}(\unicode[STIX]{x1D708}_{3})\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \boldsymbol{F}_{4}^{|W|^{2}W^{2}} & \equiv & \displaystyle \boldsymbol{B}K^{c}\big(2(\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{\ast }\big((\boldsymbol{C}\boldsymbol{x}_{1}^{W})\big(2\unicode[STIX]{x1D6FC}_{4}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{2}+3\unicode[STIX]{x1D6FC}_{3}(\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})\big)+\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{3}^{W^{3}})\big)\cdots \nonumber\\ \displaystyle & & \displaystyle +\,(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(3\unicode[STIX]{x1D6FC}_{3}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{2}^{|W|^{2}})+2\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{3}^{|W|^{2}W}))+2\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})(\boldsymbol{C}\boldsymbol{x}_{2}^{|W|^{2}})\big)\text{e}^{-2\text{i}\unicode[STIX]{x1D714}^{c}\unicode[STIX]{x1D70F}}\cdots \nonumber\\ \displaystyle & & \displaystyle -\,2\boldsymbol{B}K^{c}\unicode[STIX]{x1D708}_{3}\unicode[STIX]{x1D70F}\big(\unicode[STIX]{x1D6FC}_{1}(\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})+\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{2}\big)\text{e}^{-2\text{i}\unicode[STIX]{x1D714}^{c}\unicode[STIX]{x1D70F}}-2\unicode[STIX]{x1D708}_{3}\boldsymbol{x}_{2}^{W^{2}}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \boldsymbol{F}_{4}^{|W|^{2}} & \equiv & \displaystyle \boldsymbol{B}\big(2\unicode[STIX]{x1D6FC}_{2}\big((\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{\ast }(K^{c}(\boldsymbol{C}\boldsymbol{x}_{3}^{W})+\unicode[STIX]{x1D6FF}_{2}(\boldsymbol{C}\boldsymbol{x}_{1}^{W}))+K^{c}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{3}^{W})^{\ast }\big)\cdots \nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D6FF}_{2}(\boldsymbol{C}\boldsymbol{x}_{2}^{|W|^{2}})\big)-2\boldsymbol{B}K^{c}\unicode[STIX]{x1D70F}\,\text{Re}(\unicode[STIX]{x1D706}_{3})\big(\unicode[STIX]{x1D6FC}_{1}(\boldsymbol{C}\boldsymbol{x}_{2}^{|W|^{2}})+2\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{\ast }\big)\cdots \nonumber\\ \displaystyle & & \displaystyle -\,2\boldsymbol{x}_{2}^{|W|^{2}}\,\text{Re}(\unicode[STIX]{x1D706}_{3})\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \boldsymbol{F}_{4}^{W^{2}} & \equiv & \displaystyle \boldsymbol{B}\big(\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(2K(\boldsymbol{C}\boldsymbol{x}_{3}^{W})+\unicode[STIX]{x1D6FF}_{2}(\boldsymbol{C}\boldsymbol{x}_{1}^{W}))+\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D6FF}_{2}(\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})\big)\text{e}^{-2\text{i}\unicode[STIX]{x1D714}^{c}\unicode[STIX]{x1D70F}}\cdots \nonumber\\ \displaystyle & & \displaystyle -\,2\boldsymbol{B}K^{c}\unicode[STIX]{x1D706}_{3}\unicode[STIX]{x1D70F}\big(\unicode[STIX]{x1D6FC}_{1}(\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})+\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{2}\big)\text{e}^{-2\text{i}\unicode[STIX]{x1D714}^{c}\unicode[STIX]{x1D70F}}-2\unicode[STIX]{x1D706}_{3}\boldsymbol{x}_{2}^{W^{2}}.\end{eqnarray}$$
C.3
$O(\unicode[STIX]{x1D716}^{5})$
The list of forcing terms is:
$$\begin{eqnarray}\displaystyle & & \displaystyle -\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{3}}{\unicode[STIX]{x2202}t_{2}}+\unicode[STIX]{x1D70F}\boldsymbol{B}\left[-\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D6FF}_{2}\boldsymbol{C}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{1}}{\unicode[STIX]{x2202}t_{2}}(t-\unicode[STIX]{x1D70F})-3\unicode[STIX]{x1D6FC}_{3}K^{c}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))^{2}\boldsymbol{C}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{1}}{\unicode[STIX]{x2202}t_{2}}(t-\unicode[STIX]{x1D70F})\right.\cdots \nonumber\\ \displaystyle & & \displaystyle \quad -\,2\unicode[STIX]{x1D6FC}_{2}K^{c}\boldsymbol{C}\boldsymbol{x}_{2}(t-\unicode[STIX]{x1D70F})\boldsymbol{C}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{1}}{\unicode[STIX]{x2202}t_{2}}(t-\unicode[STIX]{x1D70F})-2\unicode[STIX]{x1D6FC}_{2}K^{c}\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F})\boldsymbol{C}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{2}}{\unicode[STIX]{x2202}t_{2}}(t-\unicode[STIX]{x1D70F})\cdots \nonumber\\ \displaystyle & & \displaystyle \quad -\left.\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{C}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{3}}{\unicode[STIX]{x2202}t_{2}}(t-\unicode[STIX]{x1D70F})+\frac{1}{2}\unicode[STIX]{x1D70F}\unicode[STIX]{x1D6FC}_{1}K^{c}\boldsymbol{C}\frac{\unicode[STIX]{x2202}^{2}\boldsymbol{x}_{1}}{\unicode[STIX]{x2202}t_{2}^{2}}(t-\unicode[STIX]{x1D70F})\right]\cdots \nonumber\\ \displaystyle & & \displaystyle \quad +\,2\unicode[STIX]{x1D6FC}_{2}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))(\boldsymbol{C}\boldsymbol{x}_{4}(t-\unicode[STIX]{x1D70F}))+2\unicode[STIX]{x1D6FC}_{2}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{2}(t-\unicode[STIX]{x1D70F}))(\boldsymbol{C}\boldsymbol{x}_{3}(t-\unicode[STIX]{x1D70F}))\cdots \nonumber\\ \displaystyle & & \displaystyle \quad +\,\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D6FF}_{2}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{3}(t-\unicode[STIX]{x1D70F}))+3\unicode[STIX]{x1D6FC}_{3}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))^{2}(\boldsymbol{C}\boldsymbol{x}_{3}(t-\unicode[STIX]{x1D70F}))\cdots \nonumber\\ \displaystyle & & \displaystyle \quad +\,3\unicode[STIX]{x1D6FC}_{3}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))(\boldsymbol{C}\boldsymbol{x}_{2}(t-\unicode[STIX]{x1D70F}))^{2}+2\unicode[STIX]{x1D6FC}_{2}\unicode[STIX]{x1D6FF}_{2}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))(\boldsymbol{C}\boldsymbol{x}_{2}(t-\unicode[STIX]{x1D70F}))\cdots \nonumber\\ \displaystyle & & \displaystyle \quad +\,4\unicode[STIX]{x1D6FC}_{4}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))^{3}(\boldsymbol{C}\boldsymbol{x}_{2}(t-\unicode[STIX]{x1D70F}))+\unicode[STIX]{x1D6FC}_{3}\unicode[STIX]{x1D6FF}_{2}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))^{3}+\unicode[STIX]{x1D6FC}_{5}K^{c}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{1}(t-\unicode[STIX]{x1D70F}))^{5}.\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
For the purpose of this study, we are interested only in the explicit expressions of the resonant forcing terms in (4.29), which read:
$$\begin{eqnarray}\displaystyle \boldsymbol{F}_{5}^{W} & \equiv & \displaystyle \unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D6FF}_{2}\boldsymbol{B}(\boldsymbol{C}\boldsymbol{x}_{3}^{W})\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{c}\unicode[STIX]{x1D70F}}\nonumber\\ \displaystyle & & \displaystyle -\,\boldsymbol{B}\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D706}_{3}\unicode[STIX]{x1D70F}\left(\unicode[STIX]{x1D6FF}_{2}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})+K^{c}(\boldsymbol{C}\boldsymbol{x}_{3}^{W})\cdots -\,{\textstyle \frac{1}{2}}K^{c}\unicode[STIX]{x1D706}_{3}\unicode[STIX]{x1D70F}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})\right)\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{c}\unicode[STIX]{x1D70F}}-\unicode[STIX]{x1D706}_{3}\boldsymbol{x}_{3}^{W}\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \boldsymbol{F}_{5}^{|W|^{2}W}\! & \equiv & \displaystyle \!\boldsymbol{B}\big(3\unicode[STIX]{x1D6FC}_{3}\unicode[STIX]{x1D6FF}_{2}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})|\boldsymbol{C}\boldsymbol{x}_{1}^{W}|^{2}+2(\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{\ast }\big(3\unicode[STIX]{x1D6FC}_{3}K^{c}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{3}^{W})\cdots \nonumber\\ \displaystyle & & \displaystyle \!+\,\left.\unicode[STIX]{x1D6FC}_{2}K^{c}(\boldsymbol{C}\boldsymbol{x}_{4}^{W^{2}})+\unicode[STIX]{x1D6FC}_{2}\unicode[STIX]{x1D6FF}_{2}(\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})\big)+K^{c}(\boldsymbol{C}\boldsymbol{x}_{3}^{W})^{\ast }\big(3\unicode[STIX]{x1D6FC}_{3}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{2}+2\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})\big)\right.\cdots \nonumber\\ \displaystyle & & \displaystyle \!+\,2\unicode[STIX]{x1D6FC}_{2}\big(K^{c}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{4}^{|W|^{2}})+K^{c}(\boldsymbol{C}\boldsymbol{x}_{2}^{|W|^{2}})(\boldsymbol{C}\boldsymbol{x}_{3}^{W})+\unicode[STIX]{x1D6FF}_{2}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{2}^{|W|^{2}})\big)\cdots \nonumber\\ \displaystyle & & \displaystyle \!+\,\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D6FF}_{2}(\boldsymbol{C}\boldsymbol{x}_{3}^{|W|^{2}W})\big)\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{c}\unicode[STIX]{x1D70F}}-\boldsymbol{B}\unicode[STIX]{x1D70F}\big(\unicode[STIX]{x1D6FF}_{2}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D708}_{3}+K^{c}\big((\boldsymbol{C}\boldsymbol{x}_{3}^{|W|^{2}W})\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D706}_{3}\cdots \nonumber\\ \displaystyle & & \displaystyle \!+\,2(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{2}^{|W|^{2}})\unicode[STIX]{x1D6FC}_{2}\unicode[STIX]{x1D706}_{3}+\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D708}_{3}((\boldsymbol{C}\boldsymbol{x}_{3}^{W})-(\boldsymbol{C}\boldsymbol{x}_{1}^{W})\unicode[STIX]{x1D706}_{3}\unicode[STIX]{x1D70F})\big)+K^{c}\big(2(\boldsymbol{C}\boldsymbol{x}_{3}^{|W|^{2}W})\unicode[STIX]{x1D6FC}_{1}\cdots \nonumber\\ \displaystyle & & \displaystyle \!+\,4(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{2}^{|W|^{2}}\!\unicode[STIX]{x1D6FC}_{2}-(\boldsymbol{C}\boldsymbol{x}_{1}^{W})\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D708}_{3}\unicode[STIX]{x1D70F})\,\text{Re}(\unicode[STIX]{x1D706}_{3})3K^{c}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})\unicode[STIX]{x1D6FC}_{3}|\boldsymbol{C}\boldsymbol{x}_{1}^{W}|^{2}(\unicode[STIX]{x1D706}_{3}+2\,\text{Re}(\unicode[STIX]{x1D706}_{3}))\cdots \nonumber\\ \displaystyle & & \displaystyle \!+\,2K^{c}(\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{\ast }(\unicode[STIX]{x1D706}_{3}+2\,\text{Re}(\unicode[STIX]{x1D706}_{3}))\big)\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{c}\unicode[STIX]{x1D70F}}\cdots \nonumber\\ \displaystyle & & \displaystyle \!-\,2\boldsymbol{x}_{3}^{|W|^{2}W}\,\text{Re}(\unicode[STIX]{x1D706}_{3})-\unicode[STIX]{x1D708}_{3}\boldsymbol{x}_{3}^{W}-\unicode[STIX]{x1D706}_{3}\boldsymbol{x}_{3}^{|W|^{2}W}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \boldsymbol{F}_{5}^{|W|^{4}W} & \equiv & \displaystyle \boldsymbol{B}K^{c}\big(12\unicode[STIX]{x1D6FC}_{4}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{2}^{|W|^{2}})|\boldsymbol{C}\boldsymbol{x}_{1}^{W}|^{2}+2(\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})^{\ast }\big(2\unicode[STIX]{x1D6FC}_{4}C^{2}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{3}\cdots \nonumber\\ \displaystyle & & \displaystyle +\,3\unicode[STIX]{x1D6FC}_{3}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})+\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{3}^{W^{3}})\big)+(\boldsymbol{C}\boldsymbol{x}_{3}^{|W|^{2}W})^{\ast }\big(3\unicode[STIX]{x1D6FC}_{3}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{+}2\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})\big)\cdots \nonumber\\ \displaystyle & & \displaystyle +\,2(\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{\ast }(3\unicode[STIX]{x1D6FC}_{3}((\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{3}^{|W|^{2}W})+(\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})(\boldsymbol{C}\boldsymbol{x}_{2}^{|W|^{2}}))+\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{4}^{|W|^{2}W^{2}}))\cdots \nonumber\\ \displaystyle & & \displaystyle +\,3((\boldsymbol{C}\boldsymbol{x}_{1}^{W})^{\ast })^{2}(4\unicode[STIX]{x1D6FC}_{4}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})+\unicode[STIX]{x1D6FC}_{3}(\boldsymbol{C}\boldsymbol{x}_{3}^{W^{3}}))\big)+\boldsymbol{B}K^{c}\big(3\unicode[STIX]{x1D6FC}_{3}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{2}^{|W|^{2}})^{2}\cdots \nonumber\\ \displaystyle & & \displaystyle +\,2\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{4}^{|W|^{4}})+2\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{2}^{|W|^{2}})(\boldsymbol{C}\boldsymbol{x}_{3}^{|W|^{2}W})\big)\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{c}\unicode[STIX]{x1D70F}}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{1}{2}\boldsymbol{B}K^{c}\unicode[STIX]{x1D70F}(\unicode[STIX]{x1D708}_{3}+2\,\text{Re}(\unicode[STIX]{x1D708}_{3}))\big(2\unicode[STIX]{x1D6FC}_{1}(\boldsymbol{C}\boldsymbol{x}_{3}^{|W|^{2}W})+4\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{2}^{|W|^{2}})\cdots \nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D708}_{3}\unicode[STIX]{x1D70F}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})+4\unicode[STIX]{x1D6FC}_{2}(\boldsymbol{C}\boldsymbol{x}_{2}^{W^{2}})(\boldsymbol{C}(\boldsymbol{x}_{1}^{W})^{\ast })+6\unicode[STIX]{x1D6FC}_{3}(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}\boldsymbol{x}_{1}^{W})(\boldsymbol{C}(\boldsymbol{x}_{1}^{W})^{\ast })\big)\text{e}^{-\text{i}\unicode[STIX]{x1D714}^{c}\unicode[STIX]{x1D70F}}\nonumber\\ \displaystyle & & \displaystyle -\,2\boldsymbol{x}_{3}^{|W|^{2}W}\,\text{Re}(\unicode[STIX]{x1D708}_{3})-\unicode[STIX]{x1D708}_{3}\boldsymbol{x}_{3}^{|W|^{2}W}.\end{eqnarray}$$
