Experiments and modelling of rate-dependent transition delay in a stochastic subcritical bifurcation.

Complex systems exhibiting critical transitions when one of their governing parameters varies are ubiquitous in nature and in engineering applications. Despite a vast literature focusing on this topic, there are few studies dealing with the effect of the rate of change of the bifurcation parameter on the tipping points. In this work, we consider a subcritical stochastic Hopf bifurcation under two scenarios: the bifurcation parameter is first changed in a quasi-steady manner and then, with a finite ramping rate. In the latter case, a rate-dependent bifurcation delay is observed and exemplified experimentally using a thermoacoustic instability in a combustion chamber. This delay increases with the rate of change. This leads to a state transition of larger amplitude compared with the one that would be experienced by the system with a quasi-steady change of the parameter. We also bring experimental evidence of a dynamic hysteresis caused by the bifurcation delay when the parameter is ramped back. A surrogate model is derived in order to predict the statistic of these delays and to scrutinize the underlying stochastic dynamics. Our study highlights the dramatic influence of a finite rate of change of bifurcation parameters upon tipping points, and it pinpoints the crucial need of considering this effect when investigating critical transitions.

Introduction

Many systems exhibit abrupt changes, or tipping, e.g. population extinction [1,2], emergence of infectious diseases [3], financial systems crisis [4], compression buckling of mechanical structures [5] and climate transitions [6-8].

Tipping is dangerous if some states of the system are associated with extreme or catastrophic events, and this explains the interest this subject has received in the last decades. Recently, different studies demonstrated that economical or environmental disasters can be modelled as dynamical systems incurring a tipping. Therefore, the development of tipping forecasting techniques with early indicators is an active research area [9-11].

A key aspect in this context is the distinction between three types of tipping, rooted in different causes [12].

B-tipping is induced by a bifurcation where the system state changes drastically for a small change of a control parameter. In this case, the tipping can be often predicted with techniques that rely on the popular concept of critical slowing down [13-17], or that make use of other properties of the attractor [18,19].

In N-tipping, dynamic noise induces jumps between several coexisting attractors (e.g. [20-23]); in this case, the analysis of the time-series statistic can help in detecting precursor of critical transitions [24,25].

R-tipping is induced by the rate at which a control parameter is varied, if several possible attractors are present in the range of parameter variation. Inertial effects play a central role in R-tipping. In the case of standard R-tipping, the system starts from an attractor but, if the parameter rate of change is larger than a critical value, it cannot follow this attractor and tips to another one [26-29]. In the case of ‘preconditioned R-tipping’, the system starts from an unstable condition and, depending on the rate of change, it evolves towards one of the possible attractors [30].

Inertial effects can also delay the bifurcation, moving the tipping point to higher/lower values of the bifurcation parameter as observed, for example, by Baer & Gaekel [31] for the FitzHugh–Nagumo model. This delay is in general a function of the parameter rate of change. Therefore, we will refer to this effect as ‘rate-delayed tipping’.

All those mechanisms can manifest independently, or, like in the present study, simultaneously. In this case, the evolution of the system results from the interplay of different time scales set by the ramp rate, the noise intensity and the system relaxation time [32]. Several examples can be found in the recent literature. Ashwin et al. [33] study the regimes of transition and the escape time in a network of bistable nodes as a function of the coupling and noise strengths. Sun et al. [34] assess the possibility of tipping for a Duffing–Van der Pol oscillator with time-delayed feedback, as a function of forcing noise intensity, feedback time delay and feedback intensity. The work from Clements & Ozgul [35] deals with two stochastic models for population dynamics, and studies the effect of the rate of change of one governing parameter on the system dynamics and on the predictability of tipping. Berglund & Gentz [36] provide theoretical and numerical analyses for rate-delayed tipping in the presence of noise in supercritical pitchfork bifurcations. An analogous study is carried out by Ritchie & Sieber [37] for a rate-dependent tipping in a saddle-node bifurcation. Kwasniok [38] introduces a method to predict a fold and a Hopf bifurcation in the presence of noise. Kuehn [39] studies the delay in a Hopf bifurcation with a random initial condition.

In this study, we show experimental evidence of simultaneous B-, N- and rate-delayed tipping mechanisms at a Hopf subcritical bifurcation, in a laboratory-scale combustor subject to thermoacoustic instabilities in the presence of turbulence-induced noise.

The four panels in figure 1 illustrate how the three types of tipping combine in our system. In these diagrams, the amplitude A is plotted as a function of the bifurcation parameter ν. In the absence of noise and for a quasi-steady change of the bifurcation parameter (figure 1a), the system state evolves on the deterministic attractor, leading to B-tipping and hysteresis (blue and red for increasing and decreasing ν). This quasi-steady picture changes if the bifurcation parameter varies at a finite rate (figure 1b): bifurcation delay occurs, and it is a function of the rate (e.g. [40-42]). For a quasi-steady variation of ν in the presence of stochastic forcing (figure 1c), the hysteresis is suppressed in a statistical sense. For each value of the bifurcation parameter, the state is defined by a probability density distribution. In this case, N-tipping occurs in the bistable region (e.g. [6,43]). Finally, when the bifurcation parameter is varied at a finite rate in the presence of stochastic forcing (figure 1d), the highest probability of state transition is delayed. This is the case discussed in this work. Our scenario, therefore, results from the combination of a finite-rate ramping through a stochastic subcritical bifurcation.

Figure 1.

Illustration of the various types of tipping encountered in the vicinity of the bistable region of a subcritical bifurcation according to the classification proposed in [12]. Solid and dashed black lines: deterministic attractor and repeller, respectively. Light to dark hues: low to high probability density. (a–c) B-tipping, rate-delayed tipping and N-tipping. (d) Present work where B-, N- and rate-delayed tipping mechanisms occur simultaneously (see figure 4 for the experimental data).

Figure 4.

Experimental evidence of the bifurcation delay and of the dynamic hysteresis in the ramped ϕ experiment. The panels are divided in ramp up (top row) and ramp down (bottom row). The stationary probability density function at seven equivalence ratios ϕ (grey) is given as a reference and compared to the evolution in time of the ensemble PDF P(A;t), where t is the time at which ϕ(t)=ϕ with ϕ(t)=ϕ0+Rt for the ramp up and ϕ(t)=ϕE−Rt for the ramp down.

This paper is organized as follows. In §2, we introduce the physical problem of thermoacoustic instabilities. In §§3.1 and 4.1, we show experimental results where the average tipping point is delayed when the control parameter is ramped at a finite rate. In §§3.2 and 4.2, we develop a low-order stochastic model of the system and demonstrate with a quantitative first-passage time analysis how the bifurcation delay statistic varies with the ramping rate. Finally, in §5 we consider a situation where a control parameter is ramped up and, if tipping is detected, ramped down in order to come back to the initial safe state. In this situation, the system may suffer irreversible damage if the ramp up is too fast, which applies to many industrial applications or to natural systems like, for instance, climate transitions.

Thermoacoustic instabilities

Thermoacoustic coupling is a phenomenon that has fascinated scientists for over two centuries. In 1777, Dr William Higgins reported, with surprise, on a hydrogen flame emitting ‘sweet tones’ if placed inside a glass tube [44]. In 1894, Lord Rayleigh provided an explanation to this observation: the gas in the tube resonates if the flame (or any other source) provides heat at the moment of maximum gas compression [45].

Many years after, during the Cold War, thermoacoustic instabilities became a very critical issue for one of the most challenging projects ever realized by humankind: the Apollo programme to take man to the Moon. As detailed in [46], the F-1 engines propelling the Saturn V had destructive combustion instabilities that required 2000 full-scale tests, with empirical modifications of the chamber geometry before the rocket was ready for take-off.

More recently, thermoacoustic instabilities became a recurrent issue in the development phase of heavy-duty gas turbines for power generation and turbofans for air transportation. This is because the resulting intense acoustic fields induce high-cycle fatigue of the combustion chambers [47]. For heavy-duty gas turbines, the pressing demand for machines with high power density and ultra-low emissions, which are capable of compensating the production intermittency of the wind and solar sources, led to the use of lean premixed flames. Unfortunately, these flames are more prone to thermoacoustic instabilities. In the case of aeroplane turbofans, these instabilities constitute an increasingly serious obstacle to the development of new aeroengines complying with more stringent emission regulations [48].

The suppression of these instabilities is very challenging due to the uniqueness and complexity of engines in real-life application [49]. Despite the achievements attained over the past decades in terms of passive mitigation implementation, development engineers cannot predict if a combustor prototype will have a sufficiently large pulsation-free operating window, over which the acoustic level is acceptable for the mechanical integrity of the components.

Figure 2a shows a schematic of our laboratory-scale combustor.[1] The air premixed with methane enters the plenum, a volume that, in practice, evens the flow delivered by the compressor and guides it to the inlet of the burner. Then, the mixture passes through the swirler, a set of curved blades that rotate the flow. This rotational motion is essential to achieve a stable anchoring of the flame. Then, the flow enters the combustion chamber, where combustion takes place. At any operating point, the fluctuating component q of the heat release rate acts as a source term in the wave equation: where p is the acoustic pressure, c the speed of sound and γ the specific heat ratio. In practice, the unsteady heat release of the flame q is influenced by the acoustic field p, via, for instance, acoustically triggered fuel supply modulation or coherent vortex shedding, which leads to a thermoacoustic feedback loop [50]. As illustrated in figure 2b, the geometry of the combustor and the temperature distribution define a set of acoustic modes in the chamber. Each mode is characterized by a shape and an eigenvalue. The latter determines whether the thermoacoustic mode is linearly stable or unstable. The system stability depends on several operating parameters, such as the mass flows of fuel or air . The transition from linearly stable to linearly unstable regime occurs at Hopf bifurcations, where the sign of the growth rate of the mode changes. If unstable, the thermoacoustic dynamics is characterized by a limit cycle, with amplitudes prms and qrms being defined by the natural acoustic damping of the chamber, and by the linear and nonlinear components of the flame response to acoustic perturbations [50,51]. The non-coherent component of the heat release rate fluctuations, which is induced by turbulence, acts as a broadband forcing on this self-sustained thermoacoustic oscillation.

Figure 2.

Thermoacoustic instabilities occur in combustion chambers for aeronautic and power-generation applications. (a) Schematic of our laboratory-scale swirled combustor. (b) Illustration of this unstable coupling. Under a certain phase difference relationship, well known as the Rayleigh criterion, a constructive feedback establishes between the unsteady heat release rate q(,t) and the acoustic field p(,t). (c) Thermoacoustic limit cycle in the laboratory-scale combustor used for this work. In the left loop, four snapshots of the flame showing the coherent motion of the flame leading to q and originating from the thermoacoustic feedback. Time traces of acoustic pressure and heat release rate are shown on the right.

A typical operating condition for which we observe a thermoacoustic limit cycle is presented in figure 2c (see also the movie in the electronic supplementary material). The four panels in the loop show instantaneous flame pictures and the corresponding phase-averaged flame shapes. The right plot displays the time traces of the acoustic pressure signal p (in red) and the heat release rate q (in blue) (note the symbols on the time trace corresponding to the four flame snapshots in the left loop). The flame exhibits a periodic motion at the frequency of the first acoustic mode (150 Hz), with sound intensity at the anti-nodes exceeding 150 dB, which is considerable for a burner operated at atmospheric pressure. This dynamic state would not be acceptable in a commercial aeronautical engine or in a heavy-duty gas turbine, because the acoustic loading, which scales with the engine operating pressure, would be highly detrimental for the mechanical components.

In this work, we focus on the transient thermoacoustic dynamics associated with the passage through the Hopf bifurcation when one of the critical operating parameters—the equivalence ratio—is ramped. We show experimental evidence of a bifurcation delay and explain the phenomenon using a surrogate low-order model. This is particularly relevant for the development of new aeronautical and land-based gas turbines, which require fast loading or deloading, and which may be at risk due to such rate-delayed tipping points.

Subcritical bifurcation

This section presents two main results. In the first part, the results of the experimental mapping of the combustor dynamics as a function of the equivalence ratio are shown. In the second part, a low-order model of the system is derived.

Stationary experiment

The combustor was operated selecting one equivalence ratio ϕ at a time. The stationary operation was reached and a long acoustic pressure signal p(t) was recorded using a microphone placed inside the chamber. The oscillation amplitude A(t) was then extracted by applying the Hilbert transform to p(t). The procedure was repeated for different equivalence ratios ϕ in the range [0.580; 0.635]. The results for five selected ϕ are presented in figure 3a. On the left, the measured acoustic pressure and amplitude signals are plotted, together with their probability density functions (PDFs) P(p) and P(A). On the right, the joint PDFs for three of the presented operative points show the statistic of the phase portraits.

Figure 3.

(a) Experimental records of the thermoacoustic subcritical Hopf bifurcation investigated in this work. According to the methane/air mixture equivalence ratio , the acoustic pressure recorded in the combustor has three different signatures, reflected in the different shapes of the PDFs P(p) and P(A). From top to bottom (increasing ϕ): small amplitude acoustic pressure resulting from the forcing of the linearly stable thermoacoustic mode by turbulence-induced noise; bistable thermoacoustic dynamics with two intermittently visited attractors; high amplitude limit cycle. These three possible regimes are also presented by the joint PDFs of the oscillation phase portrait at three exemplary ϕ. (b) Surrogate oscillator model (3.7) that mimics the subcritical Hopf bifurcation when the parameter ν is increased. In the three-dimensional plot, and three cuts, resembling the experimental . (c) The stationary PDF for the slow-varying oscillation amplitude A, obtained with the transformation of variables . On top of it, the deterministic pitchfork and saddle-node bifurcation diagram (in blue), and the stationary PDF with the corresponding potential V (A) for an overdamped particle at the three selected ν.

These results demonstrate how the system undergoes a subcritical Hopf bifurcation when the control parameter is varied. For low equivalence ratio ϕ, the system state is attracted towards zero. The small fluctuations of the acoustic signal envelope are due to the stochastic forcing exerted by the intense turbulence in the combustor. For intermediate values of ϕ, two states are possible: small amplitude acoustic pressure and high amplitude limit cycle. The intermittency between the two states is triggered by the turbulence-induced stochastic forcing (N-tipping, as in figure 1c). For higher equivalence ratio ϕ, the stochastically forced limit cycle is the only stable state. The reader can refer to the electronic supplementary material for the movies of the three regimes.

Nonlinear oscillator model

The thermoacoustic behaviour described in the previous section can be reproduced by a low-order model derived from first principles. The Helmholtz equation (2.1) (hereafter rewritten in Laplace space) defines the acoustic pressure field in the combustor, given an unsteady source of heat in the volume and impedance conditions at the boundaries: and where s is the Laplace variable, and are the transforms of acoustic pressure and velocity fluctuations, respectively, x the position, c the local speed of sound, γ the specific heat ratio, the heat release rate source term, n the outward normal to the boundary and Z the acoustic impedance. This equation is valid under low Mach number conditions.

Although nonlinear coupling among different thermoacoustic modes can occur in some practical configurations, we focus on situations where, like in the present case, one mode is dominant in the thermoacoustic dynamics. Therefore, it is possible to project the acoustic field on an orthogonal basis Ψ and approximate the system’s dynamics with the one of the dominant mode only, which will be denoted by ψ [52,53]. This yields the approximation , being the mode amplitude: where ρ is the gas density and Λ the mode normalization coefficient. This equation can be rewritten as with and Therefore, the system dynamics can be approximated by a forced damped harmonic oscillator (3.4) of resonance frequency ω0. The term α>0 represents the damping mechanisms, and it is assumed to be constant, since the impedance at the boundaries is generally a smooth function of s and therefore is not expected to vary significantly around ω0. The term is the result of the weighting on the mode shape of the volumetric heat release rate and can be decomposed into two contributions: . The first component represents the non-coherent part of the heat release rate oscillations, induced by the intense turbulence present in practical combustors. The term refers to the coherent heat release rate fluctuations, which result from a feedback interaction with the acoustic field established in the combustor. Hence, this term can be expressed as a nonlinear function of the modal amplitude η. It is customary to simplify this function by replacing it with its truncated Taylor expansion [52,53]. The linear term coefficient β of this expansion defines, together with the linear damping α, the linear stability of the system. Absorbing in the constants the mode shape ψ(x) at the pressure probe location x and considering only the odd terms of the series expansion up to the fifth order leads to the following oscillator model for the thermoacoustic system: where ν=(β−α)/2 is the oscillation linear growth rate in rad s−1 and κ and μ the two positive constants that define the nonlinear component of the oscillator response. The term ξ(t) is a white noise forcing of intensity Γ that models the non-coherent turbulence-induced heat release rate fluctuations. In figure 3b, the plot shows the stochastic Hopf bifurcation featured by this oscillator, as a function of the bifurcation parameter ν. This three-dimensional representation of the stationary joint-probability density is depicted together with three orthogonal cuts resembling the ones obtained from the experiments and showing that the bifurcation parameter ν of the surrogate model (3.7) corresponds to the equivalence ratio ϕ in the experiments.

It is convenient to describe the system evolution in terms of the slowly varying amplitude A and phase φ, with . By inserting this ansatz for p into the second-order stochastic differential equation (3.7) and by performing deterministic and stochastic averaging (e.g. [54]), one gets first-order Langevin equations for A and φ. The equation for A is , where and ζ is a white noise forcing of intensity . The deterministic dynamics derives from a potential with , and the equation does not depend on φ, which leads to the corresponding Fokker–Planck equation (FPE) for the variation in time of the amplitude PDF P(A;t): Setting ∂P/∂t=0, one obtains the stationary PDF , plotted in figure 3c as a function of the linear growth rate ν, in a pitchfork bifurcation diagram fashion. To provide a visual reference, the bifurcation diagram of the deterministic case (i.e. in the absence of noise, Γ=0) is superimposed in blue. This diagram shows the subcritical pitchfork and the saddle-node bifurcations governing the system. In the bottom insets, the PDFs for three selected values of the bifurcation parameter ν are presented. In the upper insets, the corresponding potentials are plotted. The linearly stable and stable limit cycle conditions feature a single potential well at low or high amplitude, while the bistable case has two potential wells. The stochastic forcing causes the jumps from one basin of attraction to the other, and hence the intermittency between low-amplitude noisy fluctuations and high-amplitude limit cycle oscillations.

Ramping

In this section, the dynamics of the system under transient operation is analysed. In the first part, experimental results obtained by ramping the bifurcation parameter are provided. They highlight the presence in the system dynamics of B- and N-tipping mechanisms combined with inertial and hysteresis effects. In the second part, the model introduced in §3.2 is used to study the influence of the ramp rate on the system dynamics.

Ramp experiment

The following test was performed on the test rig to highlight the peculiar dynamics of this combustor. The methane and air mass flows were controlled such that the equivalence ratio ϕ repeated 100 times the following four-step cycle: (1) linear increase for 4 s from ϕ0=0.580 to ϕE=0.635; (2) idle for 10s at ϕE; (3) linear decrease for 4 s back to ϕ0; (4) idle for 10 s at the lowest equivalence ratio. Figure 4 presents the results of this experiment. The panels are grouped in two rows: the upper row corresponds to the statistic of the 100 ramps up, the lower row to the one of the 100 ramps down. Each column corresponds to an equivalence ratio. The PDFs of this ramp experiment were obtained via a kernel density estimation (KDE) applied over the 100 realizations, and they are plotted in colour (blue for the ramp up, red for the ramp down). In all the panels, the stationary PDF for the corresponding ϕ (no ramping, already presented in figure 3a) is given in grey as a reference.

The system experiences dynamic hysteresis: in the bistable region, even though the stationary PDF features two maxima, the system stays in the low-amplitude (respectively, high-amplitude) range when ϕ is ramped up (respectively, down). Another feature is the delay in transition, easily observable in the bottom row: the dynamic PDF peak is at an amplitude that is higher than the one of the stationary PDF at the same ϕ. This means that the system experiences inertial effects, remaining close to the initial state longer: a bifurcation delay is observed. This observation corresponds to the case depicted in figure 1d.

Rate-dependent bifurcation delay

The ramp rate, together with the ramp profile, is expected to influence the bifurcation delay, as shown in [31] for a deterministic system. We therefore used the surrogate oscillator model to investigate this aspect in more detail. The parameter ν was varied linearly in time between two values ν0 and νE at different rates R. Two approaches were used. The first is to simulate (3.7) in Simulink, varying the initial condition and running different realizations of the process. Extracting the envelope for each realization and considering the ensemble statistic, it is possible to draw maps of the evolution in time of the amplitude PDF P(A;t). The other approach is to integrate numerically the FPE (3.8) and obtain directly P(A;t). The two methods closely agree, as shown in appendix A.2. In figure 5, the results of the FPE integration are presented. A ramp up/idle/ramp down/idle cycle is solved, for two different ramp rates R=50 rad s−2 and R=10 rad s−2. The dynamic stochastic bifurcation delay is captured and it is observed that a faster ramp leads to a more pronounced delayed transition from one stable point to the other.

Figure 5.

Ramping cycle of ν between the values ν0=−4.5 and νE=5.5 for the oscillator model (3.7) for two different ramp rates R. The contour plot represents the PDF P(A;t) computed by time-marching the FPE (3.8). The time is normalized with the ramp time tramp=(νE−ν0)/R. In the insets, the ramp cycle in dimensional time. The other parameters of the model are: ω0/2π= 120 s−1, κ=8 s−1, μ=2 s−1, .

An important aspect of the phenomenon depicted in this figure is that the state transitions are delayed with respect to the bifurcation point, but not time delayed (the horizontal axis in these figures is normalized by the physical duration of the cycle). In other words, a higher rate of change of the time-varying potential induces a faster transition into the neighbouring basin of attraction, but the transition occurs for a delayed value of the bifurcation parameter compared to the quasi-steady picture of the system.

First passage analysis

In this section, we imagine that a tipping point is feared due to the monotonous change of a key parameter of the system, and that one wants to ramp back this parameter sufficiently early to avoid the critical transition. In that situation, the underlying time-varying potential landscape is unknown and a controller monitors the state of the system while the parameter varies. As is usually done when new prototype engines are tested, we will take the current state of the system to feed the controller. Indeed, gas turbines and aeronautical engine combustors are equipped with a controller that constantly monitors the acoustic pressure level in the chamber. In case the measured acoustic pressure is too high, the control system intervenes, either changing the parameters to bring the operating condition back to a safe point, or in extreme cases, shutting off the flame by closing the fuel supply valve.

If the combustor features a subcritical bifurcation on the varied parameter (e.g. ϕ), the system inertia is a factor that has to be taken into account. In this case, the transition from low to high amplitudes happens suddenly and, if the bifurcation delay is long, the reached acoustic pressure level can be considerable. In this situation, the control system detects the danger late and might be ineffective in avoiding damages to the system. A way to estimate the hazard represented by the delayed bifurcation is to compute, using the surrogate model, the statistic of the time tfp needed to reach a certain danger level. This is similar to the classical problem of first passage time, often addressed in the context of bifurcation theory for stochastic dynamics in steady double-well potential [3,55-58]. A major difference in the present situation is that the potential evolves with time. Ramp rate and noise intensity are expected to influence this escape problem as theoretically shown for other types of bifurcation in [59] or [60]. The statistic of the first passage time can be computed either performing an ensemble average over many time-domain simulations of the process, or solving the unsteady Fokker–Planck equation and imposing an absorbing boundary condition at that threshold level. Details about the two methods, with results in close agreement, are provided in appendix A.3. The value ν(tth)=νth of the control parameter ν(t) at the first passage time tth is of particular interest: this quantity is proportional to the danger of the delayed transition, as it determines the limit cycle amplitude when the transition occurs. This νth statistic can be determined as νth=ν0+Rtth. The results are presented in figure 6. The contour levels represent the probability density of νth as a function of the ramp rate R. The mean value of νth (plotted in blue, 〈νth〉) increases with the ramp rate R, while the time needed to reach the danger level is shorter (see the iso-time lines). This finding indicates that a fast ramp of the control parameter is dangerous if a subcritical bifurcation is present, as exemplified in the two test cases presented in figure 7. Here the process was simulated in Simulink: the parameter ν was ramped up at two different rates R (10 and 50 rad s−2) and when the danger level was reached, ramped back down at the maximum rate R=−50 rad s−2. In figure 7a,b, many realizations of this process are presented. As a function of the initial condition and of the random excitation, each realization has a different evolution and, therefore, a different first passage time. The two extreme realizations (shortest and longest first passage times) are highlighted with thick lines. The respective deterministic bifurcation diagrams are superimposed to provide a visual reference. The PDFs obtained with a KDE over the realizations are plotted in figure 7c,d. The control system effectively brings the oscillations back to a safe level in both cases. However, the combined action of the finite ramp-down rate, dynamic hysteresis and inertia causes the system to stay in the danger zone for a certain time. The faster case R=50 m s−2 is more critical: as discussed before, the crossing of the threshold level happens on average when the target ν is already high. As a result, the system abruptly reaches high-amplitude oscillations and has to travel a long distance on the bifurcation diagram upper branch before reaching the safety zone. This effect can be gauged by comparing two quantities for the two cases R=10 and 50 rad s−2: in the latter case, the mean residence time over the safety threshold Δtth is twice larger and the mean released energy is nine times larger.

Figure 6.

P(νth;R) is the probability density of the instantaneous linear growth rate ν at the first passage over the threshold amplitude Ath, as function of the ramp rate R. It is obtained from simulations of the unsteady FPE with absorbing boundary at A=Ath. In blue 〈νth(R)〉, the linear growth rate of the system at the mean first passage time.

Figure 7.

Two exemplary cases (square R=10, circle R=50) are simulated in Simulink, implementing a control system that ramps down ν if the danger level is reached. Row (a,b) different realizations of the process (thin lines, grey and red in the safety and danger zones, respectively) and the two extreme realizations in terms of first passage time (thick grey lines) with their associated quasi-steady deterministic bifurcation diagrams (blue lines). (c,d) KDE of the PDFs. The mean residence and mean released energy in the danger zone are also indicated.

Conclusion

A subcritical Hopf bifurcation of a thermoacoustic system was investigated in this work. A laboratory-scale combustor was operated under different values of methane/air equivalence ratio, which serves as bifurcation parameter: depending on its value, acoustic pressure amplitude in the chamber is either damped, intermittently switching between low and high amplitudes, or attracted towards high-amplitude, which corresponds to a stable limit cycle. The main focus of the work was on the transient dynamics: the equivalence ratio was ramped in time and dynamic hysteresis and delayed bifurcation were observed. A nonlinear oscillator surrogate model was used to investigate the effect of the ramp rate on the bifurcation delay. It was shown that when the control parameter is ramped faster, the transition from the damped regime to the limit cycle occurs for higher values of the bifurcation parameter. The corresponding first passage problem in a time-varying potential was solved with the unsteady Fokker–Planck equation and with Monte Carlo simulations of the process. This study primarily addresses a major problem of practical combustion systems. Operating conditions of gas turbines are often varied in time, for matching power grid requirement, and similar rapid changes of the combustion regimes also occur in aeronautical engines, especially at take-off. If a subcritical thermoacoustic bifurcation is present, a delayed bifurcation results in a sudden and unexpected acoustic pressure level rise, which is detrimental to the machine integrity. Therefore, a slow variation of the machine parameters is advisable, especially when mapping the operating points of a new combustor. More broadly, this study is relevant for the countless systems which exhibit critical transitions. This work highlights the importance of carefully considering the rate of change of the bifurcation parameter, when investigating tipping points.