1. Introduction
Canonical configurations of premixed laminar flames attract significant research interest because their unsteady dynamics is rich and representative of many practical configurations, yet relatively simple in comparison with turbulent flames (Lieuwen Reference Lieuwen2003; Schuller, Poinsot & Candel Reference Schuller, Poinsot and Candel2020). A V-flame, the configuration considered in the present study, is an inverted conical flame anchored on the centrebody of an annular burner. In industrial applications, V-flames are often enhanced with swirl to stabilise the flame and prevent blow-off (Candel et al. Reference Candel, Durox, Schuller, Bourgouin and Moeck2014). However, even in the absence of swirl, the V-flame dynamics remains intricate. Flow fields corresponding to a non-swirling V-flame, computed on a two-dimensional axisymmetric grid at two different Reynolds numbers, are presented in figure 1. The flame surface region, characterised by high-magnitude heat release rates, is anchored to the centrebody. A strong free shear layer is shed from the outer corner of the jet nozzle, as illustrated by the vorticity distribution. Unlike the cylinder-anchored flame investigated in our previous study (Wang, Lesshafft & Oberleithner Reference Wang, Lesshafft and Oberleithner2022b), the superposed streamlines in the present V-flame indicate no prominent recirculation region behind the annular centrebody. As a result, the V-flame dynamics more closely resembles that of amplifier flows such as jets, resulting in increased sensitivity to acoustic perturbations (Schuller Reference Schuller2003; Schuller, Durox & Candel Reference Schuller, Durox and Candel2003). Numerous studies have been dedicated to the exploration of the linear and nonlinear dynamics of V-flames under various forcing scenarios, both in confined and unconfined arrangements (Durox, Schuller & Candel Reference Durox, Schuller and Candel2005; Birbaud et al. Reference Birbaud, Durox, Ducruix and Candel2007; Durox et al. Reference Durox, Schuller, Noiray and Candel2009). Self-excited oscillations have been observed in confined V-flames, which are furthermore prone to a chaotic dynamics (Vishnu, Sujith & Aghalayam Reference Vishnu, Sujith and Aghalayam2015). In the case of unconfined configurations, although self-excited oscillations were observed in the experiments conducted by Durox et al. (Reference Durox, Schuller and Candel2005), a comprehensive understanding of the critical conditions and the dynamics associated with instability onset remains elusive. This present work aims to elucidate the onset of axisymmetric instability in a non-swirling annular V-flame through a combined approach of linear analysis and nonlinear time-domain simulation.
Figure 1. Steady base states of the calculated V-flame at (a) 
$Re=1674$ and (b) 
$Re=2891$. Top of each plot: azimuthal vorticity 
$\varOmega _z$. Bottom of each plot: heat release rate 
$\dot {\omega }_T$. Both variables are normalised by their maximum. Streamlines are superposed.
Unsteadiness and instability in flames is a fascinating field of application for linear instability analysis because of the rich dynamics that can result from the interplay of vortical, thermal, chemical and acoustic elements. However, reactive flows are associated with steep gradients, stiff reaction terms in the governing equations and additional unknown variables for each chemical species, posing challenges for global linear stability analysis computations. As a result, many pioneering works invoked simplifications involving parallel flow and decoupled chemistry assumptions. For example, Emerson et al. (Reference Emerson, O'Connor, Juniper and Lieuwen2012) modelled the reacting wake flow behind a bluff body as an incompressible flow using discontinuous local profiles of velocity and density in the spanwise direction. Similar local analyses have also been performed on profiles measured experimentally in swirling jet flames (Oberleithner et al. Reference Oberleithner, Stöhr, Im, Arndt and Steinberg2015; Douglas et al. Reference Douglas, Emerson, Hemchandra and Lieuwen2021a), backward-facing step flames (Manoharan & Hemchandra Reference Manoharan and Hemchandra2015) and bluff-body wake flames (Emerson, Lieuwen & Juniper Reference Emerson, Lieuwen and Juniper2016), among others. However, reacting flow fields are generally strongly non-parallel in the streamwise direction and flames interact with the flow through more than just density gradients. In recent years, there has been a gradual integration of more comprehensive global linear stability analysis methods, which preserve the streamwise variations of the base state, into canonical premixed laminar flame configurations. For example, Qadri, Chandler & Juniper (Reference Qadri, Chandler and Juniper2015) investigated the self-sustained oscillations in lifted jet diffusion flames with a heat release source term included in the governing flow equations. Using compressible reacting flow equations, Avdonin, Meindl & Polifke (Reference Avdonin, Meindl and Polifke2019) calculated the eigenvalues and the associated eigenmode structures of a premixed slot flame, revealing significant fluctuations of streamwise velocity and chemical heat release rate along the flame surface region. In previous work, we investigated self-excited oscillations of a premixed laminar flame stabilised by a square cylinder in a channel (Wang et al. Reference Wang, Lesshafft and Oberleithner2022b). The study involved the examination of global eigenmodes associated with the steady base state of the reacting cylinder wake as well as nonlinear simulations, employing a one-step methane–air reaction in the low Mach number limit. The critical Reynolds number, corresponding to zero temporal growth rate, was identified, marking the onset of limit-cycle oscillations through a supercritical Hopf bifurcation. Endogeneity analysis (Marquet & Lesshafft Reference Marquet and Lesshafft2015) further indicated that this global instability was driven by momentum feedback in the wake recirculation zone, with only marginal contributions from additional feedback in the flame region, characterising the flame oscillations as a passive effect of the essentially hydrodynamic wake instability.
Other studies using global linear analysis focus on flame responses to external perturbations, including acoustic forcing, inertial waves and optimal forcing identified through resolvent analysis. Novel reduced-order models and predictive tools for the dynamic response of flames have emerged by building upon these methods. In the case of premixed slot flames, flame transfer functions have been computed using linearised reacting flow equations with one-step (Avdonin et al. Reference Avdonin, Meindl and Polifke2019; Brokof, Douglas & Polifke Reference Brokof, Douglas and Polifke2024) and two-step chemistry schemes (Meindl, Silva & Polifke Reference Meindl, Silva and Polifke2021; Wang et al. Reference Wang, Kaiser, Meindl, Oberleithner, Polifke and Lesshafft2022a), and their quantitative validation against reference results has been achieved. Resolvent analysis (Wang et al. Reference Wang, Kaiser, Meindl, Oberleithner, Polifke and Lesshafft2022a) revealed that the identified optimal excitation frequency corresponds to an intrinsic thermoacoustic mode (Silva Reference Silva2023). In the context of M-flames, computations were carried out to determine the linear response to impulsive and harmonic perturbations (Blanchard Reference Blanchard2015; Blanchard et al. Reference Blanchard, Schuller, Sipp and Schmid2015). An adjoint analysis (Skene & Schmid Reference Skene and Schmid2019) was also employed to assess the sensitivity of the M-flames’ optimal response gain to swirl. For a swirling V-flame, impulse response calculations were conducted to investigate the modulation of flame fronts by inertial waves (Albayrak, Bezgin & Polifke Reference Albayrak, Bezgin and Polifke2018). Recent studies have also calculated the linear responses of turbulent jet flames (Casel et al. Reference Casel, Oberleithner, Zhang, Zirwes, Bockhorn, Trimis and Kaiser2022; Kaiser et al. Reference Kaiser, Varillon, Polifke, Zhang, Zirwes, Bockhorn and Oberleithner2023) and a reacting jet in cross-flow (Sayadi & Schmid Reference Sayadi and Schmid2021).
The global linear stability analysis procedure employed to study the present V-flame closely follows our previous work (Wang et al. Reference Wang, Lesshafft and Oberleithner2022b) concerning a flame stabilised by a square cylinder. However, based on local stability intuition, the unsteady dynamics in these two configurations are expected to be distinct: wake recirculation as a source of local absolute instability is almost absent in the V-flame, whereas a strong convective instability resides in the jet shear layer. Moreover, the characteristic global eigenspectrum of a jet exhibits marked differences from that of a wake. In both non-reacting and reacting wakes, the global dynamics is dominated by an isolated linear eigenmode, which leads to the nonlinear shedding of counter-rotating vortices (Noack & Eckelmann Reference Noack and Eckelmann1994; Barkley Reference Barkley2006; Wang et al. Reference Wang, Lesshafft and Oberleithner2022b). In the context of round jets, such an isolated mode only exists in the presence of strong density contrast (Lesshafft et al. Reference Lesshafft, Huerre, Sagaut and Terracol2006; Coenen et al. Reference Coenen, Lesshafft, Garnaud and Sevilla2017; Chakravarthy, Lesshafft & Huerre Reference Chakravarthy, Lesshafft and Huerre2018), where self-excited oscillations have been observed (Monkewitz et al. Reference Monkewitz, Bechert, Barsikow and Lehmann1990; Kyle & Sreenivasan Reference Kyle and Sreenivasan1993; Hallberg & Strykowski Reference Hallberg and Strykowski2006). In isothermal jets, the numerical eigenspectrum is typically dominated by artificial modes that arise from spurious pressure feedback between the boundaries (Garnaud et al. Reference Garnaud, Lesshafft, Schmid and Huerre2013a; Coenen et al. Reference Coenen, Lesshafft, Garnaud and Sevilla2017). Such artificial modes do not converge with respect to domain size, and they are highly affected when a sponge layer is employed to artificially dampen pressure feedback (Cerqueira & Sipp Reference Cerqueira and Sipp2014; Lesshafft Reference Lesshafft2018).
Recent findings indicate that the onset of instability in jet flows can be subcritical in many circumstances. Zhu, Gupta & Li (Reference Zhu, Gupta and Li2017) identified a hysteretic bistable region when incrementally adjusting the jet velocity in their helium-air jet experiment. This observation suggests a subcritical Hopf bifurcation, which may be modelled through a truncated Landau equation. Demange et al. (Reference Demange, Qadri, Juniper and Pinna2022) investigated the self-sustained oscillation of a heated jet with real gas effects through numerical simulation and stability analysis. The critical temperature ratio for the Hopf bifurcation, predicted by the global stability analysis around the steady base flow, was found to be considerably lower than in nonlinear simulation, indicating subcriticality. Other recent work has identified subcritical instability dynamics in constant-density, swirling circular jets (Douglas, Emerson & Lieuwen Reference Douglas, Emerson and Lieuwen2021b) as well as swirling and non-swirling annular jets (Douglas, Emerson & Lieuwen Reference Douglas, Emerson and Lieuwen2022) using nonlinear branch tracing.
For the present V-flame configuration, we first conduct a global eigenmode analysis of the laminar reacting base flows shown in figure 1, obtained as steady solutions of the nonlinear reacting flow equations at various Reynolds numbers. A prominent branch of ‘flame-tip’ modes, arising from non-local feedback between the nozzle and the flame tip, is discussed in particular, and non-modal behaviour is briefly characterised in the framework of resolvent analysis. Nonlinear time stepping reveals a subcritical onset of instability, leading to periodic, quasi-periodic and chaotic regimes as the Reynolds number is varied. It is argued that, although the flow does not follow a simple supercritical Hopf bifurcation scenario, the limit cycle is still underpinned by a similar non-local feedback mechanism as the leading linear eigenmode. Finally, a linear eigenmode analysis of the time-averaged mean flow is attempted, to assess its capability to predict the nonlinear global frequency.
The remainder of this paper is structured as follows. In § 2, we present the V-flame configuration and describe the governing equations. Section 3 focuses on the global linear analysis around the steady base state, where various categories of eigenmodes are characterised. In § 4, nonlinear time stepping and analysis of the nonlinear flame dynamics are carried out. Section 5 presents a linear analysis around time-averaged mean flows and discusses the ambiguity of linearised chemical reaction terms in the limit-cycle regime. Conclusions are provided in § 6.
2. Calculation of steady V-flames
The burner geometry aligns with the numerical study by Birbaud et al. (Reference Birbaud, Ducruix, Durox and Candel2008) on an inverted dihedral flame, but the burner considered here is axisymmetric, formulated in the cylindrical coordinates 
$\boldsymbol {x}=(x, r)$. Here, 
$x$ represents the streamwise direction, and 
$r$ represents the radial direction. Figure 2(a) illustrates the entire calculation domain. The streamwise position of the outer corner of the nozzle is defined as 
$x=0$, and the symmetry axis is designated as 
$r=0$. Lean premixed methane–air reactant is injected from an annular nozzle measuring 30 mm in length and 
$D=11$ mm in diameter. The V-flame is anchored on a centrebody with a diameter of 3 mm, protruding 2 mm from the nozzle. The domain extends to 
$x_{{max}}=200$ mm downstream and 
$r_{{max}}=50$ mm radially. An axisymmetric condition is applied at the centreline to confine the calculation to a two-dimensional setting.
Figure 2. The standard mesh with 180 481 elements used in the study. Colours represent the normalised base flow heat release. (a) The complete numerical domain. (b) Zoom on a flame surface region.
The governing equations for the reacting flow, considering an ideal gas in the low Mach number limit, are the same as in Wang et al. (Reference Wang, Lesshafft and Oberleithner2022b), except that they are formulated here in axisymmetric coordinates. Primitive variables (
$\rho, \boldsymbol {u}, h, Y_{\mathrm {CH}_{4}}, p$) are chosen, where 
$\boldsymbol {u}=(u_x,u_r)$ represents the streamwise and radial velocity components, 
$\rho$ denotes density, 
$h$ stands for sensible enthalpy, 
$Y_{\mathrm {CH}_{4}}$ is the mass fraction of methane and 
$p$ is the hydrodynamic pressure. The governing equations are expressed as
$$\begin{gather} \frac{\partial \rho}{\partial t} ={-}\boldsymbol{\nabla}\boldsymbol{\cdot}\left (\rho \boldsymbol{u} \right ), \end{gather}$$
$$\begin{gather}\rho\frac{\partial \boldsymbol{u}}{\partial t}={-} \rho \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} - \boldsymbol{\nabla} p + \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\tau}, \end{gather}$$
$$\begin{gather}\rho\frac{\partial Y_{\mathrm{CH}_{4}}}{\partial t} ={-}\rho \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} Y_{\mathrm{CH}_{4}}+ \boldsymbol{\nabla} \boldsymbol{\cdot}\left (D_s \boldsymbol{\nabla}Y_{\mathrm{CH}_{4}} \right ) + \dot{\omega}_{\mathrm{CH}_{4}}, \end{gather}$$
$$\begin{gather}\rho\frac{\partial h}{\partial t} ={-} \rho \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} h + \boldsymbol{\nabla}\boldsymbol{\cdot} \left (D_h \boldsymbol{\nabla} h \right ) + \dot{\omega}_T, \end{gather}$$
$$\begin{gather}p_0=R_s \rho T. \end{gather}$$
The molecular stress tensor 
$\boldsymbol{\tau}$ is defined as
\begin{equation} \boldsymbol{\tau} =\mu \left ( \boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} \boldsymbol{u}^{\mathrm{T}} - \tfrac{2}{3} \left (\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} \right ) \boldsymbol{I} \right ), \end{equation}
where the molecular viscosity follows Sutherland's law 
$\mu =A_s T^{1/2} /(1+T_s/T)$ with constants 
$A_s=1.672\times 10^{-6}\ \mathrm {kg}\ (\mathrm {m}\ \mathrm {s}\ \mathrm {K}^{1/2})^{-1}$ and 
$T_s=170.7$ K. The species and enthalpy diffusive transport coefficients, 
$D_s$ and 
$D_h$, respectively, are assumed to be proportional to the molecular viscosity, with constant values of Schmidt number 
${Sc}=\mu /D_s=0.7$ and Prandtl number 
${Pr}=\mu /D_h=0.7$. Enthalpy is expressed as 
$h=C_pT$, where the specific heat capacity 
$C_p=1.3\ \mathrm {kJ}\ (\mathrm {kg}\ \mathrm {K})^{-1}$ is assumed constant. Following the original low Mach number expansion of McMurtry et al. (Reference McMurtry, Jou, Riley and Metcalfe1986), the thermodynamic pressure 
$p_0$ in the ideal gas law (2.5) is the zeroth-order pressure component, prescribed as 
$p_0=101.3$ kPa, while the unknown hydrodynamic pressure 
$p$ in the momentum equation (2.2) is its first-order complement in the squared Mach number. Also, 
$R_s=264.6\ \mathrm {J}\ (\mathrm {kg}\ \mathrm {K})^{-1}$ in the ideal gas law represents the specific gas constant. The one-step reaction scheme 
$\mathrm {1S\_CH4\_MP1}$ from CERFACS (2017) is used, as it is sufficient to accurately reproduce the laminar burning speed of lean premixed methane–air mixtures with the current set of parameters (Birbaud et al. Reference Birbaud, Ducruix, Durox and Candel2008). This scheme was also employed in our previous calculations of a bluff-body flame, successfully reproducing the length of the recirculation bubble compared with reference calculations (Wang et al. Reference Wang, Lesshafft and Oberleithner2022b). The reaction rate 
$\mathcal {Q}$ is governed by an Arrhenius law
\begin{equation} \mathcal{Q}=A_r[X_{\mathrm{CH}_{4}}]^{n_\mathrm{CH_4}}[X_\mathrm{O_{2}}]^{n_\mathrm{O_2}} \exp \left ( -\frac{T_{a}}{T} \right ). \end{equation}
For the lean global methane–air reaction modelled here, the values of the reaction exponents are 
$n_\mathrm {CH_4}=1$ and 
$n_\mathrm {O_2}=1/2$. With these exponents, the values of the Arrhenius pre-exponential factor and activation temperature are, respectively, 
$A_r=1.1\times 10^{7} \mathrm {m}^{3/2}\ ({\rm s}\ {\rm mol}^{1/2})^{-1}$ and 
$T_a=1.007\times 10^4$ K (CERFACS 2017). The molar concentrations of 
$\mathrm {CH_4}$ and 
$\mathrm {O_2}$ are given by 
$[X_{\mathrm {CH_4}}]=\rho ({Y_{\mathrm {CH_4}}}/{W_{\mathrm {CH_4}}})$ and 
$[X_{\mathrm {O_2}}]=\rho ( {Y_{\mathrm {O_2}}}/{W_{\mathrm {O_2}}})$, respectively, where 
$W_\mathrm {CH_4}=16.0\ {\rm g}\ {\rm mol}^{-1}$ and 
$W_\mathrm {O_2}=32.0\ {\rm g}\ {\rm mol}^{-1}$ represent their molecular masses. We assume a very lean mixture such that 
$Y_{\mathrm {O_2}}=0.2128$ is constant. The reaction rate in the species equation is denoted by 
$\dot {\omega }_\mathrm {CH_4}=-W_{\mathrm {CH_4}}\mathcal {Q}$, and the heat release due to combustion in the enthalpy equation is expressed as 
$\dot {\omega }_T=-\Delta h_{f}^o \mathcal {Q}$, where 
$\Delta h_{f}^o=-804.1\ \mathrm {kJ}\ \mathrm {mol}^{-1}$ is the standard enthalpy of reaction.
Figure 3 shows the base flow streamwise velocity and temperature fields at a bulk velocity of 
$U_0=3.0\ \mathrm {m}\ \mathrm {s}^{{-1}}$. For this condition, the corresponding inflow Reynolds number defined as 
$Re = \rho _0 U_0 D/\mu _0$ is equal to 2282. A simple parabolic profile is prescribed as the inlet velocity. The inflow conditions for temperature and fuel mass fraction are set to 
$T_0=300$ K and 
$Y_{\mathrm {CH_4}}=0.04256$. No-slip and adiabatic conditions are imposed on the inflow channel walls. The dump-plane wall (
$x=0$) connected to the outer corner of the annular nozzle is assumed to be no slip and isothermal at 
$T_0=300$ K. Note that the use of a slow co-flow, as employed for stabilisation purposes in Birbaud et al. (Reference Birbaud, Ducruix, Durox and Candel2008), is not considered here. No-slip and isothermal conditions are also imposed on the annular centrebody, with the wall temperature fixed at 
$T_r=700$ K, as in the reference. At the lateral boundary far from the flame (
$r=r_{{max}}$) a free-slip condition is imposed on the velocity, and the temperature and fuel mass fraction are set to 
$T_0=300$ K and 
$Y_{\mathrm {CH_4}}=0$, respectively. A symmetry condition is imposed along the centreline. Finally, a traction-free condition is employed at the downstream boundary (
$x=x_{{max}}$), with Neumann conditions on the temperature and fuel mass fraction.
Figure 3. Streamwise velocity and temperature associated with the steady base state at 
$Re=2282$.
In the following, we investigate V-flames in a range of bulk velocities from 
$U_0=2.2\ \mathrm {m}\ \mathrm {s}^{-1}$ to 
$U_0=3.8\ \mathrm {m}\ \mathrm {s}^{-1}$. The prescribed conditions result in 
$Re=\rho _0 U_0 D/\mu _0$ ranging from 1674 to 2891, where 
$\rho _0$ and 
$\mu _0$ denote the density and molecular viscosity at the inflow. Figure 1 shows the base flame shape and vorticity distribution at both ends of the Reynolds number range investigated. The flame is elongated downstream as the velocity is increased.
Sponge layers are introduced at the downstream (
$x=x_{{max}}$) and lateral boundaries (
$r=r_{{max}}$) to prevent spurious back scattering (Lesshafft Reference Lesshafft2018). The sponge layer implementation from Meliga, Sipp & Chomaz (Reference Meliga, Sipp and Chomaz2010) is employed, wherein viscosity is artificially increased in the sponge regions. The molecular viscosity 
$\mu$ is divided by the sponge layer expression 
$sg(r,x)$, formulated as
\begin{equation} \left.\begin{gathered} sg(r,x) = 1\quad \text{if } r\leqslant r_{{sg}}\quad \text{and}\quad x\leqslant x_{{sg}},\\ sg(r,x) = 1+(\alpha-1) \zeta(x,x_{{sg}})\quad \text{if } r\leqslant r_{{sg}}\quad \text{and}\quad x > x_{{sg}},\\ sg(r,x) = sg(r_{{sg}},x)+\left[\alpha-sg(r_{{sg}})\right]\zeta(r,r_{{sg}})\quad \text{if } r > r_{{sg}},\end{gathered}\right\} \end{equation}
where 
$\alpha ^{-1}$ is listed in table 1. Here, 
$x_{{sg}}$ and 
$r_{{sg}}$ denote the starting position of the sponge layers in the streamwise and radial directions. The function 
$\zeta$, defined by
\begin{equation} \zeta(a,b)=\frac{1}{2}+\frac{1}{2}\tanh \left\{\gamma \tan \left(-\frac{\rm \pi}{2}+{\rm \pi} \frac{|a-b|}{l} \right) \right\}, \end{equation}
yields a smooth transition of 
$\mu$ from the inner region to the truncated boundaries. In this expression, the constant 
$\gamma$ is set to 1 and 
$l$ denotes the sponge layer thickness, equal to 
$x_{{max}}-x_{{sg}}$ or 
$r_{{max}}-r_{{sg}}$ in the streamwise and radial directions, respectively. It is important to note that the transport coefficients 
$D_s$ and 
$D_h$ are also increased in the sponge layer as they are proportional to 
$\mu$. This implementation is applied both in the computation of the base flow and in the subsequent linear analysis.
Table 1. Sponge layer parameters and number of mesh elements for eigenmode calculations in figure 4.
Figure 4. Eigenspectra at 
$Re=2282$. (a) Convergence with respect to parameters of sponge layers. (b) Convergence with respect to grid resolution. Legends are given in table 1. Different categories of eigenmodes are marked with colours: flame-column modes (green ‘
’) in figures 5(a)–5(c), plume modes (blue ‘
’) in figures 5(d)–5(f), flame-tips modes (red ‘
’) in figure 6 and flame-tip-downstream modes (brown ‘
’) in figure 5(g–i).
Figure 5. Three families of eigenmodes at 
$Re=2282$. Radial velocity fluctuations are shown. (a–c) Flame-column modes: first three modes with green markers ‘
’ in figure 4. (d–f) Plume modes: first three modes with blue markers ‘
’. (g–i) flame-tip-column modes: first three modes with brown markers ‘
’. Black contour: an isocontour of base state vorticity illustrating the free shear layer position; yellow contour: an isocontour of base state heat release rate illustrating the flame front position.
Figure 6. Branch of flame-tip modes at 
$Re=2282$. First three modes (a–c) with red markers ‘
’ in figure 4. Top of each plot: heat release rate fluctuation. Bottom of each plot: radial velocity fluctuation. Black contour: an isocontour of base state vorticity illustrating the free shear layer position; yellow contour: an isocontour of base state heat release rate illustrating the flame front position.
The nonlinear governing equations are discretised on the unstructured mesh presented in figure 2 using 180 481 triangular elements. High grid resolution is allocated at the flame surface, as shown in figure 2(b), with a maximum mesh resolution of 
$\Delta x=0.08$ mm. A continuous Galerkin method is employed using mixed Taylor–Hood finite elements, of quadratic order for the velocity and linear order for other flow variables. The base states are calculated by Newton iteration. Readers may refer to Appendix B of Wang (Reference Wang2022) for more details on the numerical methods.
3. Global linear analysis
3.1. Survey of the global eigenspectrum
To analyse the onset and behaviour of self-sustained oscillations in V-flames, a global linear stability analysis is conducted by computing the eigenmodes associated with the Jacobian matrix of the governing equations, which are linearised around the steady base state. Equations (2.1)–(2.4) are linearised with respect to the primitive variables (
$\rho, \boldsymbol {u}, h, Y_\mathrm {CH_4}, p$) at each grid point of the entire numerical domain. The flow fluctuations 
$\boldsymbol {q'}(\boldsymbol {x},t)$ can be expanded in the basis of eigenmodes 
$\boldsymbol {q}'_j(\boldsymbol {x},t)=\boldsymbol {\phi }_j(\boldsymbol {x})\exp (\lambda _j t)$ obtained from the generalised eigenvalue problem
\begin{equation} \lambda_j \boldsymbol{\mathsf{B}}\boldsymbol{\phi}_j =\boldsymbol{\mathsf{A}}\boldsymbol{\phi}_j. \end{equation}
The matrices 
$\boldsymbol{\mathsf{B}}$ and 
$\boldsymbol{\mathsf{A}}$ are constructed from the linearisation of (2.1)–(2.4). The eigenvalues 
$\lambda _j$ and associated eigenvectors 
$\boldsymbol {\phi }_j$ are computed using the eigs() function in Matlab. The frequency is defined as 
$2{\rm \pi} f_j=\textrm {Im}[\lambda _j]$, and the growth rate is defined as 
$2{\rm \pi} \sigma _j=\textrm {Re}[\lambda _j]$, so that 
$\lambda _j/(2{\rm \pi} )=\sigma _j+\textrm {i}f_j$. Furthermore, their non-dimensional counterparts are defined as the Strouhal number 
${St}=fD/{U_0}$ and 
${\tilde {\sigma }}={\sigma D}/{U_0}$. Instabilities associated with non-zero azimuthal numbers are not considered in this study, although they are known to arise spontaneously in isothermal annular jets (Douglas et al. Reference Douglas, Emerson and Lieuwen2022) or in conical flames (Douglas, Polifke & Lesshafft Reference Douglas, Polifke and Lesshafft2023).
The linear eigenmodes associated with steady base states at Reynolds numbers ranging from 
$Re=1674$ to 
$Re=2891$ with an increment of 
$\Delta Re=76$ are calculated (corresponding to an inflow velocity increment of 
$\Delta U_0=0.1\ \mathrm {m}\ \mathrm {s}^{-1}$). Different families of eigenmodes are identified based on their dependence on boundary conditions, their eigenmode structures, and their trends with respect to the Reynolds number. Specifically, their convergence with respect to parameters of the sponge layers is examined, by using the different values given in rows (a–c) of table 1. From case (
$a$) to (
$c$), the sponge layers are increasingly large and viscous. The sponge layer in case (
$c$) corresponds to around one half of the whole domain length, with a sixtyfold increased molecular viscosity at the boundaries. The results are presented in figure 4(a), where the eigenmodes that overlap with the three different markers are considered converged.
The eigenfunctions associated with the identified categories are shown in figures 5 and 6, marked with different cross colours in figure 4(
$a$), using the parameters of case (
$a$) in table 1. In each figure, the position of the flame surface is illustrated by a yellow isocontour corresponding to 10 % of the maximum heat release rate, while the shear is illustrated by a black isocontour corresponding to 33 % of the maximum azimuthal vorticity. The eigenmode spectra at Reynolds numbers equal to 1978, 2282, 2586 and 2891 are presented in figure 7. Four main families of eigenmodes are identified:
(i) Flame-column modes (green crosses ‘
’ in figure 4, eigenmodes shown in 5a–c). These modes are characterised by their low frequencies; they are the least stable family of eigenmodes for $Re<2586$. These modes do not undergo destabilisation as the Reynolds number increases but are notably influenced by the parameters of the sponge layer. The eigenmode structures exhibit spatial growth primarily in the streamwise direction. Specifically, the mode associated with the lowest frequency, illustrated in figure 5(a), displays maximum amplitude at the downstream end of the computational domain. This fluctuation expands radially to both sides of the shear layer, with the inner side extending close to the centreline. At higher frequencies, the peak fluctuation moves upstream, and spreads closer to both sides of the shear layer. Such modes, characterised by an amplitude maximum near the domain's downstream end, are commonly encountered in incompressible and compressible jets, where they are attributed to the stable advection of nearly neutral structures within the shear layer (Garnaud et al. Reference Garnaud, Lesshafft, Schmid and Huerre2013a). In general, these modes do not converge with larger domain sizes and are influenced by the sponge layers at the boundary (Lesshafft Reference Lesshafft2018). Given that the fluctuations fill the column of the plume, we denote these modes as flame-column modes, akin to jet-column modes.(ii) Plume modes (blue crosses ‘
’ in figure 4, eigenmodes associated with the three lowest frequencies shown in 5d–f). Along this branch, identified at relatively low frequencies, the temporal growth rates exhibit an increasing trend with frequency. Remarkably, these modes remain unaffected by the sponge layer, suggesting that they form a family distinct from the flame-column modes. Although the branch experiences a slight destabilisation with growing Reynolds numbers, all modes within this range of $Re$ remain stable. The fluctuation structures are primarily localised slightly downstream of the flame surface and within the jet shear region, corresponding to areas of heightened base vorticity. Notably, the maximum perturbation is consistently observed around $x/D=4$ for the three leading eigenmodes, positioned closely along the shear in the radial direction. This specific location corresponds to the point where the base flow velocity and temperature profiles develop to be parallel in the streamwise direction, creating a base flow reminiscent of a non-reacting hot jet. Similar eigenmode structures were previously identified in the mixing layer of a plume, as illustrated in figure 4(b) of Chakravarthy et al. (Reference Chakravarthy, Lesshafft and Huerre2018). In their study of plumes, the maximum fluctuation amplitude is found very close to the inflow, corresponding to the region with the maximum density gradient in the base flow and confined within the mixing layer. Despite the absence of buoyancy effects in the present governing equations, we refer to this branch as plume modes.(iii) Flame-tip modes (red crosses ‘
’ in figure 4, leading eigenmodes shown in figure 6, red crosses ‘’ in figure 7). This branch constitutes a prominent feature in V-flames, notably separated from other more stable eigenvalues in the spectra. The eigenmodes along this branch remain unaffected by the sponge layer. However, with an increase in Reynolds number, certain members of this branch become unstable. Specifically, the mode associated with the lowest frequency becomes unstable at $Re=2815$, suggesting a Hopf bifurcation. Subsequently, we refer to this mode as the ‘leading flame-tip mode’. The fluctuation of heat release rate and velocity is confined to the inner mixing layer, extending closely downstream from the intersection of the flame surface and the jet shear. The maximum fluctuation is identified in the flame extinction zone at around $x/D=3$ for the leading flame-tip mode. For higher-frequency flame-tip modes, the maximum is located further upstream, closer to the flame surface, and the associated fluctuations exhibit higher wavenumbers.(iv) Flame-tip-column modes (brown crosses ‘
’ in figure 4, eigenmodes shown in 5g–i). This branch encompasses stable eigenmodes across a broad range of frequencies. The eigenmodes exhibit fluctuations at the flame extinction region and downstream close to the centreline. Notably, the fluctuation amplitudes at the downstream end become more pronounced with higher frequency. However, these modes are sensitive to the sponge layer parameters and remain unaffected by changes in Reynolds number.
’ in 
’ in 
’ in 
’ in 
’ in 
Figure 7. Eigenspectra at different Reynolds numbers. Branches of flame-tip modes are marked as red. The dimensional temporal growth rates and frequencies are used for comparison: (a) 
$Re=1978$; (b) 
$Re=2282$; (c) 
$Re=2586$; (d) 
$Re=2891$.
The convergence of eigenmodes with respect to the mesh refinement is examined at 
$Re=2282$ using meshes with 180 481, 221 019 and 251 757 elements. Figure 4(b) shows that convergence has been achieved on the reference mesh with 180 481 elements for all the eigenmodes described above except certain plume modes. The poor convergence of the plume modes seems to result from the fact that the maximum fluctuation of those plume modes are located at around 
$x/D=4$ where the change of local refinement is relatively abrupt. Nonetheless, as the least stable plume mode is converged with respect to refinement on the reference mesh and its associated growth rate is negative, we conclude that further refinement is not necessary for our global instability study.
To summarise, the flame-column modes and the flame-tip-column modes exhibit intense oscillations at the end of the numerical domain, with their associated eigenvalues significantly influenced by sponge layer parameters. This suggests that these modes arise from spurious pressure feedback from the outflow, rather than from physical instabilities. In contrast, the plume and the flame-tip modes remain unaffected by the sponge layers, indicating a more physical nature. The plume modes appear to result from the extended shear layer located downstream from 
$x/D=4$, sharing similar base flow and eigenmode structures with a non-reacting plume. Conversely, the flame-tip modes, characterised by strong oscillations at the flame extinctions around 
$x/D=3$, become considerably more destabilised as the Reynolds number increases, leading to a positive growth rate at 
$Re=2891$. The leading flame-tip mode, i.e. the least stable mode along the flame-tip branch, demonstrates an increased frequency with the Reynolds number, as shown in figure 7.
3.2. Analysis of the flame-tip mode mechanisms
We now aim to characterise the physical mechanisms governing the behaviour of the leading flame-tip modes, drawing inspiration from earlier experimental work by Schuller (Reference Schuller2003) and Durox et al. (Reference Durox, Schuller and Candel2005). In their studies of perturbed laminar V-flames, the flame oscillations were interpreted as a consequence of vortex structures, advected from the burner exit along the shear layer, modulating the flame shape. Building upon this interpretation, the time delay between the heat release rate and an imposed velocity perturbation at the inlet can be linked to the convective velocity of the vortices and the convection distance. Experimental measurements of this time delay were performed, and the convective velocity was estimated as one half of the maximum streamwise velocity at the burner exit. The convection distance was then determined by multiplying the time delay by the convective velocity, revealing a correspondence with a region where the roll-up of the flame surface was notably pronounced.
In contrast to these acoustically forced experimental investigations, our focus is specifically toward the frequency of the leading flame-tip eigenmode, which is representative of the intrinsic flame dynamics in the absence of external forcing. When the leading eigenmode is unstable, self-excited flame oscillations are expected. In our approach, we posit that the eigenfrequency of a flame-tip mode in a V-flame is linked to the advection time from the nozzle exit to the flame surface along the outer shear layer. The assumption of linearity underlines that the base flow remains unaffected by perturbations, ensuring steady streamlines and a steady flame shape in the base state. This allows for a more precise calculation of the convective velocity and distance compared with the finite-amplitude (i.e. nonlinear) oscillating flame scenarios explored experimentally by Schuller (Reference Schuller2003) and Durox et al. (Reference Durox, Schuller and Candel2005). However, the calculation is not without challenges, primarily due to the following factors. Firstly, the start- and endpoints of the convective distance are not clearly defined. Secondly, the shear layer exhibits a continuum of advection velocities across its finite thickness, and the precise streamline along which the vortices are advected remains to be chosen. Finally, hydrodynamic perturbations are generally dispersive, and are not advected by the base flow in a completely passive manner.
To address these challenges, we rely on two key assumptions in our calculation. First, we presume that the eigenmode oscillation stems from a non-local interaction between known up- and downstream points. The upstream interaction initiates at the nozzle exit (
$x=0$), while the downstream interaction concludes at the point where the streamline intersects with the flame front. Second, we assume that the perturbations in the shear layer traverse through the region marked by the highest base flow vorticity without dispersion. Based on these assumptions, the advection time is computed from the base flow in three sequential steps:
(i) We identify the points corresponding to the strongest vorticity in the shear layer, denoted by ‘
$\triangle$’ in figure 8(a).(ii) A backward integration is executed from ‘
$\triangle$’ based on the base flow velocity field. The starting point where the obtained streamline intersects with the burner exit at $x=0$ is then determined and marked as ‘$\blacksquare$’.(iii) A forward integration is initiated from ‘
$\blacksquare$’, through ‘$\triangle$’. Along this trajectory, we record the Lagrangian heat release rate, as illustrated in figure 8(b). The intersection between the streamline and the flame surface is defined by the peak of the Lagrangian heat release rate. This point is interpreted as the downstream terminus of the non-local interaction and is marked as ‘$\bullet$’. The advection time is subsequently determined as the time lag along the streamline between ‘$\blacksquare$’ and ‘$\bullet$’.
Figure 8. Calculation of Lagrangian advection time and comparison with leading flame-tip mode frequencies. (a) Steady base states at 
$Re=2282$. Top: vorticity 
$\varOmega$. Bottom: heat release rate 
$\dot {\omega }_T$. Both variables are normalised with their maximum. Marker ‘
$\triangle$’ represents the point associated with maximum vorticity. Marker ‘
$\blacksquare$’ represents the intersection of backward integration and 
$x=0$. Marker ‘
$\bullet$’ represents the peak of Lagrangian heat release rate along the streamline stemming from ‘
$\blacksquare$’. (b) Lagrangian heat release as a function of time 
$t$ normalised by 
$D/U_0$. The advection time is obtained as the time lag between ‘
$\blacksquare$’ and ‘
$\bullet$’. (c) Evolution of the Strouhal number with Reynolds number. Marker ‘
’ denotes the leading flame-tip modes. Marker 
$\circ$ represents the Strouhal number calculated from the average advection time over all starting grid points where the vorticity exceeds 99 % of its global maximum magnitude. The associated error bars are shown.
This process is repeated for all points in the shear layer where the vorticity exceeds 99 % of its global maximum magnitude. Consequently, a collection of differing advection times is obtained for each base flow at a specified Reynolds number. In figure 8(c), the black dots depict the averaged values of the Strouhal number, calculated as the inverse of the advection time, with error bars indicating the maximum and minimum within each set. The red dots represent the Strouhal number associated with the leading flame-tip eigenmode. Remarkably, over the extensive range of Reynolds numbers investigated, the frequency of the leading flame-tip mode closely aligns with the Lagrangian advection time. It is crucial to again note that the calculation does not consider the dispersion of perturbations, and instead assumes that perturbations travel at the base flow velocity. This assumption, along with spatial variations in the effective centre of the up- and downstream interaction regions, may explain the small errors observed at certain Reynolds number ranges.
This outcome suggests that the identified flame-tip modes indeed originate from a non-local feedback between the nozzle exit and the flame surface with a frequency characterised by the advection time of downstream-travelling hydrodynamic perturbations along the outer shear layer. It can be hypothesised that the non-local feedback loop is closed by upstream-travelling pressure waves generated by the fluctuations at the flame tip. In the low Mach number limit, this pressure impulse is felt instantaneously at the nozzle outlet.
3.3. Resolvent analysis
When the flame is subjected to external perturbations, its linear response may be dominated by pseudo-resonance (Trefethen & Embree Reference Trefethen and Embree2005), owing to the non-normality of the operator. This possibility is assessed through resolvent analysis, performed on the base flow at 
$Re=2282$ and 
$Re=2586$, with volume forcing in the axial and the radial momentum equations. The response is calculated using the same boundary conditions as in the global eigenspectrum calculations, and the standard 2-norm of velocity is employed both for the forcing and response norms (Garnaud et al. Reference Garnaud, Lesshafft, Schmid and Huerre2013b). The optimal gain 
$\varGamma ^2$ relating the forcing and response norms at different Strouhal numbers 
${St}$ are presented in figure 9, where the frequencies corresponding to the flame-tip eigenmodes in figures 7(b) and 7(c) are illustrated as the red dashed lines. At both Reynolds numbers, the dominant resolvent gain peaks align well with flame-tip mode frequencies. Figure 10 presents the optimal response structures at the local peak frequencies 
${St}=0.31$ and 
${St}=0.68$ related to 
$Re=2586$, which are characterised by structures that are clearly similar to the flame-tip eigenmodes displayed in figure 6. These results indicate that the linear flame dynamics is dominated by modal mechanisms, even in the presence of flow forcing.
Figure 9. Optimal gain curves at 
$Re=2282$ (a) and 
$Re=2586$ (b) obtained by the resolvent analysis. The red dashed lines represent the frequency of the flame-tip eigenmodes in figure 7.
Figure 10. Resolvent response structures at 
$Re=2586$ for 
${St}=0.31$ (a) and 
${St}=0.68$ (b). Top of each plot: heat release rate fluctuation. Bottom of each plot: radial velocity fluctuation. Black contour: an isocontour of base state vorticity illustrating the free shear layer position; yellow contour: an isocontour of base state heat release rate illustrating the flame front position.
4. Nonlinear analysis
4.1. Time-series analysis
The outcomes derived from the linear analysis indicate a Hopf bifurcation at 
$Re=2815$, as evidenced by the temporal growth rate evolution associated with the flame-tip modes in figure 11. In this section, we delve into the nonlinear dynamics of the V-flame by superimposing finite-amplitude perturbations on the steady base states at each 
$Re$ and conducting nonlinear time integration using a first-order implicit scheme. The results depicted in figure 12 correspond to a linearly stable base state at 
$Re=1978$. We denote the unperturbed base flow velocity 
$\boldsymbol {u}_b (\boldsymbol {x})$. The superposed perturbations are prescribed as the velocity components of the leading flame-tip eigenmode, represented by 
$\boldsymbol {u}_p(\boldsymbol {x})$. Both components of the velocity perturbation are normalised such that their maximum is the same maximum amplitude associated with the streamwise base flow velocity. The perturbed initial velocity field 
$\boldsymbol {u}_0 (\boldsymbol {x})$ can be expressed as
\begin{equation} \boldsymbol{u}_0 (\boldsymbol{x})=\boldsymbol{u}_b (\boldsymbol{x})+\beta\boldsymbol{u}_p(\boldsymbol{x}), \end{equation}
where 
$\beta$ is a ratio constant to be varied. We consider two perturbation amplitude ratios of 
$\beta =20\,\%$ and 
$\beta =50\,\%$. In each scenario, the temporal signal of the streamwise velocity is recorded at 
$(x/D,r/D) = (2.7,1.2)$ in a region proximal to the flame extinction zone. Figure 12 also presents snapshots of the flame surface at the time steps marked with red dots in the time signal.
Figure 11. Normalised temporal growth rate of the leading flame-tip mode with Reynolds number.
Figure 12. Illustration of subcritical behaviour at 
$Re=1978$. Streamwise velocity signal recorded at 
$(x/D,y/D) = (2.7,1.2)$ with initial velocity perturbations of (a) 
$\beta =20\,\%$ and (b) 
$\beta =50\,\%$ magnitude. The snapshots of flame movement are shown in (a i–a iv) and (b i–b iv), corresponding to the time steps marked with red dots ‘
’. A zoom-in of the time series is given in (b) over one period of the limit-cycle oscillation.
At the initial perturbation amplitude of 
$\beta =20\,\%$, the velocity perturbation undergoes transient growth before it enters modal decay. Only a small section of the flame surface near the downstream edge (
$x/D>2$) displays visible oscillations during the transient growth, after which the flame surface converges to the steady base state. This transient growth is an indicator of the strong non-normality of the system, which may promote bypass to other sub- or non-critical attractors, as in classical Poiseuille and Couette flows (Schmid & Henningson Reference Schmid and Henningson2001). Indeed, at a larger initial amplitude of 
$\beta =50\,\%$, the perturbation undergoes growth before settling into a limit-cycle oscillation. Four snapshots in an oscillation cycle are presented in figure 12(b i–b iv). A considerable portion of the flame surface (
$x/D>1$) exhibits oscillations, resulting in pronounced wrinkling along the length of the flame as well as dramatic flapping of the flame tip, resembling the oscillating flame surfaces observed in Durox et al. (Reference Durox, Schuller and Candel2005). The results indicate that the system has at least two distinct nonlinear attractors at 
$Re=1978$.
The bifurcation diagram with respect to the Reynolds number is depicted in figure 13, revealing the presence of a bi-stable region between 
$Re=1750$ and 
$Re=2815$. Both forward and backward continuation paths, characterised by the streamwise velocity perturbation 
$u'_x$, are displayed in the figure, and the methods for tracking both paths are distinct.
Figure 13. Bifurcation diagram. The backward continuation path is represented with black dots (
$\boldsymbol {{\cdot }}$), corresponding to the local maxima of streamwise velocity perturbation signal measured. The stable range of forward continuation path is represented by red squares (
). The Hopf bifurcation point of the steady state is represented by a blue diamond (
).
Forward path: the forward continuation path along the steady base state is traced by introducing small-amplitude velocity fluctuations of the leading flame-tip mode onto the steady states at each investigated Reynolds number. The initial perturbation is set to the velocity of the leading eigenmode with a small amplitude of 1 % of the maximum base flow streamwise velocity. For 
$Re<2815$, all flames disturbed in this way reconverge to the original steady states, akin to the time series in figure 12(a). These steady states are represented by the red squares (
) in figure 13. At 
$Re=2815$, marked by the blue diamond (
), the disturbed flame undergoes temporal growth before entering an unsteady fluctuating state. A similar dynamic process is observed for the disturbed steady flame at 
$Re=2891$. At both Reynolds numbers, the simulation duration corresponds to around six flow-through times of the computational domain including the sponge layer. This long duration was necessary to ensure subsidence of the initial transient, such that the fluctuations converge to a statistically stationary state. The resulting loss of stability of the steady state apparent from the time-domain simulations aligns with the critical 
$Re$ for instability identified via global linear analysis in § 3.
Backward path: the backward continuation path is identified by gradually reducing the Reynolds number with 
$\Delta Re=76$, starting with the unsteady fluctuating state at the next-highest 
$Re$. (For example, the 
$Re=2815$ case is initialised with the final state from the 
$Re=2891$ case). The black dots (
$\boldsymbol {{\cdot }}$) in figure 13 represent any local maxima identified in the velocity fluctuation signal at the probe location 
$(x/D,r/D) = (2.7,1.2)$. The cutoff time horizon before recording the local maxima corresponding to three flow-through times for the Reynolds number 
$2130 \leqslant Re \leqslant 2815$. The values of the local maxima are normalised by the temporal average of the velocity time series, representing the relative fluctuation amplitudes. The results show that fluctuations persist even when 
$Re$ is decreased below the critical Reynolds number to 
$Re=2663$. This indicates that the system exhibits subcritical dynamics and possesses at least two attractors below the critical point. The distribution of local maxima markedly narrows when the Reynolds number is reduced to 2510, indicating a change in the dynamic state. The distribution then further narrows when reducing the Reynolds number from 2510 to 2130. At 
$Re=2130$, the distribution of local maxima collapses essentially to a single value, with the local maxima varying by less than 1 % of the sliding average of ten successive local maximum samples. The presence of a single local maximum value is consistent with the limit-cycle oscillation behaviour observed in figure 14(a). Once the Reynolds number is reduced below 1750, the self-sustained oscillations vanish, and the system converges again to the steady state, as shown by the collapse of red squares and black dots at 
$Re=1598$ and 
$Re=1674$, the lower limit of our investigation.
Figure 14. (a,d,g) Time signal of streamwise velocity perturbation measured at 
$(x/D,r/D) = (2.7,1.2)$. (b,e,k) Three-dimensional Poincaré sphere. (c,f,i) Poincaré section. Panels show (a–c) 
$Re=1978$, (d–f) 
$Re=2358$ and (g–i) 
$Re=2815$.
The dynamic states associated with the oscillating flame at 
$Re=1978$, 
$Re=2358$ and 
$Re=2815$ are further characterised by plotting the temporal signal of 
$u'_x(t)$ at 
$(x/D,r/D) = (2.7,1.2)$, the associated three-dimensional Poincaré trajectory computed with a time delay of 
$\Delta t=5\times 10^{-3}$ s and the corresponding two-dimensional Poincaré section in figure 14. Figures 14(b) and 14(c) illustrate that the associated phase space at 
$Re=1978$ is an enclosed trajectory with two points depicted in the intersection plane, indicating a limit-cycle state. At 
$Re=2358$, the temporal signal displays a nearly enclosed trajectory and plane intersection points scattered in clusters, a pattern representative of a quasi-periodic dynamics. Finally, unlike the organised state-space nature identified for the lower Reynolds numbers, the temporal signal for 
$Re=2815$ reveals erratic and intermittent behaviours, corresponding to the local maxima in figure 13 being densely distributed from zero to the maximum amplitude. In particular, the Poincaré section of figure 14(f) exhibits scattered points that suggest a chaotic nature of the flame fluctuations. Overall, this progression is consistent with a Reulle–Takens–Newhouse scenario for the onset of chaos in the annular V-flame.
Power spectral density (PSD) contours, computed along the backward continuation path, are depicted in figure 15, tracking the axial velocity perturbation at 
$(x/D, r/D) = (2.7,1.2)$ and the global heat release rate. The abscissa ranges from 
$Re=1750$, the lower extent of the identified bi-stable region, to the upper limit represented by the critical 
$Re$ at the white dashed line. The linear flame-tip mode frequencies obtained from the imaginary part of the leading linear eigenvalue are overlaid on the PSD. Tracking the point velocity evolution reveals a single fundamental frequency peak for 
$Re<2586$, distinct from the flame-tip modes. At 
$Re=2586$, a second incommensurate frequency peak, closely aligned with the linear flame-tip mode frequency, emerges. This suggests that a 2-torus is born in the phase space via a Neimark–Sacker bifurcation, and the system has transitioned from periodicity now to two-frequency quasi-periodicity. In contrast to the two-frequency quasi-periodicity presented in a previous forced synchronisation study (Guan et al. Reference Guan, Gupta, Kashinath and Li2019), where quasi-periodicity arises from the competition between a self-excited mode and a forced mode, this one arises from the competition between two incommensurate self-excited modes. Further increases in Reynolds number result in more continuously distributed frequency spectra, indicative of the non-periodic state observed at 
$Re>2586$. Examining the global heat release rate unveils dominant frequency components at twice the frequency observed when tracking point velocity. This indicates nonlinear harmonic interactions – a notable feature of the flame dynamics that cannot be captured by linear analysis.
Figure 15. Power density spectra along the backward continuation path by tracking (a) local streamwise velocity and (b) the global heat release rate. The frequencies of the leading linear flame-tip mode are superposed as white dots. The white dashed lines represent the Reynolds number associated with Hopf bifurcation. The spectra are displayed in the Reynolds number range 
$1750 \leqslant Re \leqslant 2891$. No flame oscillation is observed for 
$Re \leqslant 1674$ : (a) 
$\hat{u}_{x}$; (b) 
$\hat{\dot {\omega }}_{T}$.
In summary, increasing 
$Re$ along the forward continuation path confirms that even small-amplitude perturbations lead to unsteady oscillations above the critical Reynolds number identified through linear analysis. Subcritical dynamics corresponding to hysteresis is identified along the backward continuation path when gradually reducing the Reynolds number, resulting in progressive transitions among unsteady dynamic states before the end of the hysteresis interval. At 
$Re=2510$, close to the critical Reynolds number at 
$Re=2815$, the dominant nonlinear oscillation frequency peak aligns with the linear flame-tip mode frequency, although other frequencies remain present. The nonlinear frequency associated with the subcritical oscillation below 
$Re=2510$ is apparently unrelated to the flame-tip mode frequencies identified in the linearly stable subcritical base flows. A qualitative interpretation of the subcritical dynamics can be hypothesised from the flame shape snapshots in figure 12: nonlinear oscillations are self-sustaining when a sufficiently large portion of the downstream flame edge is perturbed with sufficiently large amplitude and vice versa.
4.2. Analysis of the nonlinear non-local interaction
Having already explored the linear dynamics in § 3.2, we here investigate the physical mechanics of the subcritical limit-cycle oscillation at 
$Re=1978$ in the nonlinear case. As in figure 15, the system exhibits a strong limit-cycle oscillation in this regime with a fundamental frequency of 
$f_0=110.0$ Hz at 
$Re=1978$. Informed by the linear non-local mechanism explored above, we hypothesise a nonlinear non-local scenario where vortex structures advected from the burner exit interact with the flame surface. Thus, as before, we posit that the limit-cycle frequency can be estimated by the convection time from the burner exit to the flame surface. However, unlike the linear case, in the nonlinear scenario, the unsteady flame surface exhibits substantial movement. Hence, the convective time delay between when a vortex structure departs from the burner exit and when it arrives at the flame surface depends strongly on the starting phase, denoted as 
$\phi _{{start}}=2{\rm \pi} {}t_{{start}}/T_0$, where 
$T_0=1/f_0$ is the oscillation period and 
$t_{{start}}$ is the starting time instant. The value of convective time delay from the burner exit at 
$x=0$ to the downstream flame surface also depends on the exact starting radial position of a measured trajectory, denoted as 
$r_{{start}}$.
We conduct a Monte Carlo simulation to generate pathlines initiated at different 
$t_{{start}}$ and 
$r_{{start}}$. A numerical non-inertial particle tracer solver is implemented and coupled with the nonlinear flame solver. The tracer trajectories are computed using the same value of discretisation time as the flame solver, and the particle is advanced using a fourth-order Runge–Kutta scheme. The starting time 
$t_{{start}}$ is considered at each time step over one cycle period 
$T_0$, resulting in 92 points of phase 
$\phi _{{start}}$. Regarding 
$r_{{start}}$, we prescribe a distribution centred at 
$r=0.436D$, corresponding to the mean radial position of trajectories experiencing the strongest vorticity magnitudes identified in the linear regime. The prescribed 
$r_{{start}}$ distribution covers the range from 
$0.409D$ to 
$0.464D$ with 31 points. Hence, the overall number of sampling points for the Monte Carlo simulation is 
$N_t\times N_r = 92\times 31=2852$.
The obtained 
$N_t\times N_r$ passive tracer trajectories are analysed. The convective time delay from the burner exit to the flame surface, denoted as 
$\Delta t_{{conv}}$, is identified as the difference between the peak instant of Lagrangian heat release rate and the starting time in the same manner as in § 3.2. To identify the most relevant tracer trajectories, a selection process is conducted, based on an assumption that only the trajectories experiencing sufficient magnitudes of vorticity and heat release rate along the trajectories are essential to the global flame dynamics. The vorticity magnitude criterion is designed to select trajectories corresponding to the advection of vortex structures, while the heat release rate criterion serves to select trajectories that actually cross the flame surface. A vorticity magnitude threshold 
$\varOmega _{{th}}$ and a heat release rate threshold 
$\dot {\omega }_{T,{th}}$ are defined, respectively, to select the trajectories and estimate the oscillation frequency 
$f_e$. A formula for the estimation of 
$f_e$ is proposed as
\begin{equation} f_e=\frac{N_t N_r}{\sum_{i=1}^{N_t} \sum_{j=1}^{N_r} \delta_{i,j}}, \end{equation}where
\begin{equation} \delta_{i,j}=\begin{cases} \Delta t_{{conv},i,j}, & \text{if } \max(\dot{\omega}_{T,i,j})>\dot{\omega}_{T,{th}} \text{ and } \max( \varOmega_{i,j})>\varOmega_{{th}}, \\ 0, & \text{else.} \end{cases} \end{equation}
The formula implies that 
$f_e$ is simply estimated as the inverse of the mean convective time delay of the selected trajectories. The indices 
$i,j$ represent the variables of a trajectory of the sampling space with 
$i \in [1,N_t]$ and 
$j \in [1,N_r]$. Also, 
$\max (\dot {\omega }_{T,i,j})$ refers to the maximum heat release rate along the entire trajectory, whereas 
$\max (\varOmega _{i,j})$ refers to the maximum vorticity magnitude along the trajectory portion 
$x/D>0.3$. (The latter restriction prevents large vorticity magnitudes localised at the burner lip from interfering with the identification of advecting vortex structures.)
Figure 16(a) illustrates the selected sampling points corresponding to a prescribed set of threshold 
$\varOmega _{{th}}=0.89 \varOmega _{max}$ and 
$\dot {\omega }_{T,{th}}=0.94 \dot {\omega }_{T,max}$, where 
$\varOmega _{max}$ and 
$\dot {\omega }_{T,max}$ represent the maximum vorticity magnitude and heat release rate recorded among all the trajectories of the selected sampling points. Each point is associated with a set of (
$t_{start}$, 
$r_{start}$), coloured by the value of convective time delay 
$\Delta t_{conv}$. The trajectories at smaller 
$r_{start}$ (i.e. close to the annular jet core) experience smaller convective time delays before reaching the flame surface, and vice versa.
Figure 16. (a) Selected trajectories using a set of threshold 
$\varOmega _{{th}}=0.89 \varOmega _{max}$ and 
$\dot {\omega }_{T,{th}}=0.94 \dot {\omega }_{T,max}$. Each trajectory is associated with a set of starting instant and starting radial position (
$t_{start}$, 
$r_{start}$), coloured by the value of convective time delay 
$\Delta t_{conv}$. (b) Differences between the estimated frequency 
$f_e$ and the actual nonlinear frequency 
$f_0$ using different sets of selection thresholds (
$\dot {\omega }_{T,{th}}$, 
$\varOmega _{{th}}$). Differences 
$|\,f_e-f_0|$ larger than 10 Hz are represented by the same degree of red or blue colours. The black circle marker corresponds to the selection criterion used in (a). The dashed rectangle in the top right corner represents the trajectories filtered by the highest values of 
$\dot {\omega }_{T,{th}}$ and 
$\varOmega _{{th}}$ where 
$|\,f_e-f_0|<5$ Hz.
The estimated frequency 
$f_e$ obtained by (4.2) on different sets of threshold values (
$\dot {\omega }_{T,{th}}$, 
$\varOmega _{{th}}$) is presented in figure 16(b). The contour colours indicate the difference between 
$f_e$ and the reference frequency 
$f_0=110.0$ Hz extracted from the nonlinear simulation. The bottom right of this figure corresponds to high values of 
$\dot {\omega }_{T,{th}}$ and low values of 
$\varOmega _{{th}}$. Hence, this region selects trajectories closer to the jet core, where tracers are more likely to reach the flame surface but are also likely to experience smaller vorticity. As noted above, due to the higher fluid velocities near the jet core, these trajectories also experience smaller 
$\Delta t_{conv}$ and therefore larger values of 
$f_e$. Conversely, the top left part of figure 16(b) corresponds to low values of 
$\dot {\omega }_{T,{th}}$ and high values of 
$\varOmega _{{th}}$. This region tends to contain trajectories starting near the burner wall with higher vorticity but also with less likelihood to intersect the flame surface at a point where the heat release is large. Likewise, since these trajectories are initialised in the slow-moving fluid near the wall, they feature longer convective time delays and lower values of 
$f_e$. The trajectories with high values of both 
$\dot {\omega }_{T,{th}}$ and 
$\varOmega _{{th}}$ lead to a relatively small difference between 
$f_e$ and 
$f_0$. More specifically, the absolute difference 
$|\,f_e-f_0|$ obtained over the range 
$\dot {\omega }_{T,{th}}/\dot {\omega }_{T,{max}} \in [0.79,0.94]$ and 
$\varOmega _{{th}}/\varOmega _{{max}} \in [0.86,0.90]$ is less than 5 Hz within the dashed rectangle on the top right corner of figure 16(b). (Further increase in the direction of either 
$\dot {\omega }_{T,{th}}$ or 
$\varOmega _{{th}}$ leads to a sharp decline in the number of available sampling points, leading to deviated values of 
$f_e$.) Hence, the trajectories that yield the best frequency estimates according to (4.2) both (i) exhibit prominent vortex structures and (ii) intersect the flame surface where and when the local instantaneous heat release rate is large. This supports our hypothesis that these processes are essential to the flame dynamics.
The sampling points illustrated in figure 16(a) correspond to the black circle marker in figure 16(b) at high threshold values of 
$\dot {\omega }_{T,{th}}$ and 
$\varOmega _{{th}}$. Using these points as input for (4.2), 
$f_e$ is identified as 109.8 Hz, very close to the reference value of 
$f_0=110.0$ Hz. The time instant when a trajectory intersects the flame surface can be calculated as 
$t_{int}=t_{start}+\Delta {}t_{conv}$. We find that the standard deviation of 
$t_{int}$ calculated from the different trajectories in figure 16(a) is only 
$2.4\,\%$ of a cycle period, indicating that these trajectories reach the flame surface at nearly the same time instant. Figure 17 presents the snapshot of heat release rate and vorticity magnitude at the time instant associated with the mean value of 
$t_{int}$. Comparing this snapshot against the snapshots in figure 12, we note that the point where the pathlines associated with the strongest vorticity intersect the flame surface is situated at a position relatively close to the upstream limit of the flame flapping motion. Further, a structure of significant vorticity magnitude can be visualised close to the flame tip in figure 17(b). The trajectory at 
$t_{{start}}/T_0=0.30$ and 
$r_{{start}}/D=0.45$, corresponding to 
$\Delta {}t_{{conv}}/T_0=0.99$, is superposed in figure 17. Note that this pathline trajectory does not correspond to any streamline of the instantaneous flow field, because the flow is unsteady. Compared with the base flow streamline identified within linear regime in figure 8(a), the axial position where the trajectory intersects the flame surface is significantly reduced. The result aligns with the higher oscillation frequency 
$f_0=110.0$ Hz observed in the nonlinear regime in comparison with the linear flame-tip mode identified at 87.8 Hz.
Figure 17. Snapshots of heat release rate (a) and vorticity magnitude (b) at the time instant of mean intersection time 
$t_{int}$ corresponding to the selected trajectories in figure 16(a). The trajectory at 
$2{\rm \pi} {}t_{{start}}/T_0=0.30$ and 
$r_{{start}}/D=0.45$ is superposed: (a) 
$\dot {\omega }_{T}$; (b) 
$|\varOmega|$.
5. Global linear analysis of the nonlinear mean flow
Finally, we compute the global eigenmodes of the time-averaged mean flow to investigate the potential recovery of the nonlinear oscillation frequency and/or structure by these mean flow eigenmodes. From the outset, this procedure is complicated by the non-unique definition of a mean flow, exposed by Karban et al. (Reference Karban, Bugeat, Martini, Towne, Cavalieri, Lesshafft, Agarwal, Jordan and Colonius2020), that arises from nonlinearity. Those authors demonstrated the issue by comparing resolvent analysis results for a compressible jet, obtained from averaging the same large-eddy simulation data in either primitive or conservative variables. The present reacting flow case provides an even more compelling illustration of mean-flow ambiguity. We will consider two equally plausible definitions of the mean reaction rate 
$\mathcal {Q}$, noting that many more are possible. Representing the time average by an overbar, the first definition is the average of 
$\mathcal {Q}$ itself
\begin{equation} \bar{\mathcal{Q}}(X_i, T)=\overline{A_r[X_\mathrm{CH_{4}}]^{n_\mathrm{CH_4}}[X_\mathrm{O_{2}}]^{n_\mathrm{O_2}} \exp \left ( -\frac{T_{a}}{T} \right )}. \end{equation}
The second definition inserts the mean-flow variables into the definition of 
$\mathcal {Q}$
\begin{equation} \mathcal{Q}(\bar{X}_i, \bar{T})=A_r[\bar{X}_\mathrm{CH_{4}}]^{n_\mathrm{CH_4}}[\bar{X}_\mathrm{O_{2}}]^{n_\mathrm{O_2}} \exp \left ( -\frac{T_{a}}{\bar{T}} \right ). \end{equation}
The mean reaction rate, assessed by both definitions, is depicted in figure 18. Here, 
$\bar {\mathcal {Q}}(X_i, T)$ reveals pronounced oscillations of the flame surfaces, leading to the progress rate being distributed in a region around flame extinction. Conversely, 
$\mathcal {Q}(\bar {X}_i, \bar {T})$ displays only a thin flame surface and fails to represent the unsteady oscillations. It is worthwhile to note that, in turbulent reacting flows, the difference between two reaction rates is often employed to assess the turbulence–chemistry interaction, which is important in turbulent reaction modelling (Poinsot & Veynante Reference Poinsot and Veynante2005; Duan & Martín Reference Duan and Martín2011; Di Renzo & Urzay Reference Di Renzo and Urzay2021). Conventionally, the second definition 
$\mathcal {Q}(\bar {X}_i, \bar {T})$, called the laminar chemistry or laminar reaction rate model, only takes the frozen mean-flow quantity into account, as if in a steady laminar flame. Conversely, the first definition 
$\bar {\mathcal {Q}}(X_i, T)$, called the turbulent reaction rate, includes the species and temperature fluctuations modulated by the turbulent flow field. Figure 18 shows that even in a laminar flame, unsteadiness can lead to significant differences between the mean reaction rates evaluated by these definitions. Thus, the conventional notion of a laminar chemistry model requires caution.
Figure 18. Time-averaged reaction rate 
$\mathcal {Q}(Y_{\mathrm {CH_4}},T)$ (top) and reaction rate of time-averaged variables 
$\mathcal {Q}(\bar {Y}_{\mathrm {CH_4}},\bar {T})$ (bottom). The variables are normalised by their maximum.
In the linearised reacting flow equations, the expression for the linearised reaction rate, 
$\mathcal {Q}'$, is derived as
\begin{equation} \mathcal{Q}'=\bar{\mathcal{Q}}\left ( (n_\mathrm{CH_4}+n_\mathrm{O_2})\frac{\rho'}{\bar{\rho}}+n_\mathrm{CH_4}\frac{Y'_\mathrm{CH_{4}}}{\bar{Y}_\mathrm{CH_{4}}}+ n_\mathrm{O_2}\frac{Y'_\mathrm{O_{2}}}{\bar{Y}_\mathrm{O_{2}}} +\frac{T_a T'}{\bar{T}^2} \right ), \end{equation}
where the mean reaction rate 
$\bar {\mathcal {Q}}$ is factored out and provided by the second expression 
$\mathcal {Q}(\bar {X}_i, \bar {T})$. Note that, in accordance with our fuel-limited reaction model, 
$Y'_{\mathrm {O_2}}=0$ by assumption.
A mean-flow analysis is conducted at 
$Re=1978$, where the saturated unsteady dynamics corresponds to a subcritical limit-cycle oscillation. The Fourier mode is extracted at the limit-cycle fundamental frequency of 110.0 Hz, as identified in figures 14(a)–14(c) and 15. Figure 19(a) presents the associated Fourier mode of progress rate fluctuation, serving as a reference for fluctuation structures. These structures manifest as progress rate wrinkles convecting along the distributed mean reaction zone, as assessed through 
$\bar {\mathcal {Q}}(Y_{\mathrm {CH_4}},T)$. Concurrently, figure 21(a) illustrates the Fourier mode of radial velocity, revealing oscillations in both the flame region and the extended shear layer. The identified Fourier mode structures underscore the distinctive characteristics of the subcritical nonlinear dynamics in the V-flame compared with the linear flame-tip modes. Note that the extracted frequency 110.0 Hz is not the dominant frequency of the heat release rate, which appears at the first harmonic (cf. figure 15b).
Figure 19. Reaction rate fluctuation 
$\mathcal {Q}'$ at 
$Re=1978$. (a) Fourier modes at the limit-cycle fundamental frequency of 110.0 Hz, obtained from nonlinear simulation. (b) A priori test at 110.0 Hz. (c) Computed mean-flow eigenmode corresponding to the eigenvalue at 
$(108.8-45.8\mathrm {i})$ Hz. (d) Computed mean-flow eigenmode corresponding to the eigenvalue at 
$(159.1-0.7\mathrm {i})$ Hz.
The outcomes of the a priori assessments for 
$\mathcal {Q}'$ utilising 
$\bar {\mathcal {Q}}=\mathcal {Q}(\bar {Y}_{\mathrm {CH_4}},\bar {T})$ in its formulation are displayed. Fourier modes of the flow variables 
$Y'_\mathrm {CH_{4}}$ and 
$T'$ are inserted into (5.3), following the methodology employed for a priori tests in turbulent reaction models (Kaiser et al. Reference Kaiser, Varillon, Polifke, Zhang, Zirwes, Bockhorn and Oberleithner2023). The structure of 
$\mathcal {Q}'$ acquired using 
$\bar {\mathcal {Q}}=\mathcal {Q}_(\bar {Y}_{\mathrm {CH_4}},\bar {T})$ is depicted in figure 19(b), showing travelling waves along the thin flame surface. These structures are akin to the shape of mean-flow chemical progress rates; however, they deviate markedly from the reference Fourier mode structures.
The resulting mean-flow spectrum presented in figure 20 reveals a separated branch of eigenvalues, akin to the flame-tip modes identified for the base flow in § 2. The two leading eigenvalues on this branch exhibit temporal growth rates close to zero, but their associated frequencies are distinct from the 110.0 Hz fundamental frequency of the nonlinear limit-cycle state, as indicated by the dashed blue line. The mean-flow eigenmode structures associated with the eigenvalues marked with red are also displayed. Of these, one occurs at 
$(108.8 -45.8\mathrm {i})$ Hz, a heavily damped eigenvalue with a frequency close to the fundamental tone of the nonlinear oscillation. Another is the leading, marginally damped eigenvalue on the separated branch at 
$(159.1 -0.7\mathrm {i})$ Hz. The mean-flow modes of progress rate fluctuation in figures 19(c)–19(d) and those of radial velocity fluctuation in figures 21(b)–21(c) were found not to align with the reference Fourier modes. The strong velocity fluctuations downstream in figure 21(b) indicates that the eigenvalue 
$(108.8 -45.8\mathrm {i})$ Hz is similar to a flame-column mode, and it is not relevant to the fundamental frequency of the nonlinear oscillation. The structure of velocity fluctuation of the leading eigenvalue 
$(159.1 -0.7\mathrm {i})$ Hz in figure 21(c) shares a certain resemblance to the Fourier mode in figure 21(a), although the identified eigenvalue frequency is notably different.
Figure 20. Eigenspectrum of the time-averaged mean flow at 
$Re=1978$. Blue dashed line indicates the limit-cycle frequency of the velocity perturbation at 110.0 Hz. Eigenvalues marked in red are 
$(108.8 -45.8\mathrm {i})$ Hz and 
$(159.1 -0.7\mathrm {i})$ Hz, for which the corresponding eigenmode structures are given in figures 19 and 21, respectively.
Figure 21. Radial velocity fluctuation 
$u'_r$ at 
$Re=1978$. (a) Fourier modes at the limit-cycle fundamental frequency of 110.0 Hz, obtained from nonlinear simulation. (b) Computed mean-flow eigenmode corresponding to the eigenvalue at 
$(108.8 -45.8\mathrm {i})$ Hz. (c) Computed mean-flow eigenmode corresponding to the eigenvalue at 
$(159.1-0.7\mathrm {i})$ Hz.
At this point, we emphasise that linear eigenmode analysis of the time-averaged mean flow is not expected, in principle, to predict the nonlinear limit-cycle frequency. For example, Sipp & Lebedev (Reference Sipp and Lebedev2007) have shown through weakly nonlinear analysis that its success depends on whether nonlinear interactions are strongly dominated by mean-flow distortion effects, provided these effects are stabilising (i.e. that the bifurcation is supercritical). Building on this work and others (e.g. Karban et al. Reference Karban, Bugeat, Martini, Towne, Cavalieri, Lesshafft, Agarwal, Jordan and Colonius2020; Beneddine et al. Reference Beneddine, Sipp, Arnault, Dandois and Lesshafft2016; Tammisola & Juniper Reference Tammisola and Juniper2016), the failure of mean-flow analysis in the present V-flame can be attributed to its particular instability dynamics which has been identified and analysed in the previous sections.
Foremost, when the relevant bifurcation is subcritical, there should be no expectation of validity for mean-flow analysis even if the Sipp & Lebedev (Reference Sipp and Lebedev2007) criterion on the Stuart–Landau coefficients is satisfied. Unlike the supercritical case, subcritical dynamics indicates that harmonic interactions induce a fundamentally destabilising effect in the vicinity of the bifurcation (i.e. for sufficiently small amplitudes). This may be because mean-flow distortion is destabilising or unsteady harmonic feedback is destabilising, or both. In subcritical systems, saturation does not become possible until the amplitude reaches some finite values, indicating that harmonic interactions play a different role on the limit cycle in comparison with the base state. The ‘real-zero imaginary-frequency’ property of a mean flow (Turton, Tuckerman & Barkley Reference Turton, Tuckerman and Barkley2015) results from the growth of an unstable linear eigenmode that gradually distorts the mean flow until it can grow no further; clearly, in a subcritical situation where linear growth is impossible, that scenario cannot play out.
Additionally, the oscillations identified in the limit-cycle regime are strongly dichromatic, and the non-local instability mechanism crucially depends on interactions between features associated with both frequencies. While there is no reason why a global mean-flow analysis could not deal with non-local interactions on their own, the unsteady heat release rate in this case cannot be accurately represented at the fundamental frequency, as evidenced by the a priori test in figure 19(b). Indeed, as is common of premixed flames exhibiting strong flame–flow interactions (see, for example, the seminal observations of Joos & Vortmeyer Reference Joos and Vortmeyer1986; Lang Reference Lang1991), the global heat release rate oscillates predominantly at the first harmonic of the fundamental velocity oscillation tone (see figure 15). This harmonic generation phenomenon can be understood intuitively, since the global heat release is proportional to the flame surface area, which increases twice per cycle – at both extremes of a flow oscillation Lieuwen (Reference Lieuwen2005). Conversely, as discussed by Tammisola & Juniper (Reference Tammisola and Juniper2016), linear analysis assumes that all fields (including the heat release rate) oscillate at a single global frequency. Hence, the harmonic interaction between the heat release rate and flow velocity in oscillating premixed flames cannot be captured by linear mean flow analysis – this interaction is fundamentally nonlinear.
A strong intrinsic nonlinearity of the mean reaction rate was also encountered in our prior work on flames anchored in the wake of a two-dimensional square cylinder in a channel (Wang et al. Reference Wang, Lesshafft and Oberleithner2022b). However, the mean-flow analysis in that work accurately captured the oscillation frequency. This difference can be mainly attributed to the distinct instability mechanisms in a cylinder wake flame and an annular V-flame. The dominant dynamics in the former case arises from a supercritical hydrodynamic shear instability (the Bénard–von Kármán instability) localised in the wake recirculation zone, with only secondary influences from the spatially separated flame front region. Conversely, the investigated V-flame dynamics is subcritical and there is a significant non-local flame–flow interaction involving two distinct frequencies, as discussed in § 4.2.
6. Conclusions
This study computationally investigates the self-excited axisymmetric oscillations of a lean premixed V-flame in a laminar annular jet. The reactive flow is simulated using an irreversible single-step global chemistry model representing a lean premixed methane–air reaction coupled to the low Mach number compressible Navier–Stokes equations. Following the identification of steady states of the linearised reacting flow equations, we conduct a detailed survey of the axisymmetric global eigenmodes computed around these base states. For sufficiently high Reynolds number, destabilisation occurs for an eigenmode on an ‘arc branch’ separated from other families of more stable eigenmodes. These arc branch modes, which we term ‘flame-tip’ modes, are characterised by strong fluctuations near the flame tip and are independent of numerical domain truncation. A detailed, quantitative description of the linear feedback mechanism driving their destabilisation is provided by associating the frequency of the leading flame-tip mode with the Lagrangian advection time along the outer shear layer from the nozzle exit to the flame tip. A non-modal resolvent analysis demonstrates that the receptivity of the flame to forcing of the linear operator is largely consistent with simple resonance of the linear eigenmodes. Significant pseudo-resonance is not observed.
Upon assessing these linear results against nonlinear time-domain simulations, however, a more complex picture emerged. For small initial perturbations, the onset of sustained unsteadiness corresponds to the critical Reynolds number identified by linear analysis. Conversely, for sufficiently large perturbations, self-sustained oscillations occur even at Reynolds number values where the flame is linearly stable, revealing a substantial interval of hysteresis. Continuation analysis along this branch of unsteady solutions reveals an ordered sequence of state transitions in the subcritical regime. Along most of the unsteady branch, the unsteady flow settles into a limit-cycle state with a periodicity that does not match any linear eigenmodes of the base flow along the steady solution branch. However, as the Reynolds number approaches the critical value for linear instability of the steady state, the dynamics becomes enriched by an ordered sequence of increasingly high-dimensional features, including apparent quasi-periodicity and chaos. Interestingly, the frequency associated with the leading (stable) eigenmode of the base state becomes prominently represented in the power density spectrum during this process, suggesting a Neimark–Sacker bifurcation arising from two competitive modes. This dynamics is consistent with a Reulle–Takens–Newhouse scenario for the onset of chaos in the V-flame. Together, these findings shed new light on the nonlinear dynamical elements underpinning the observed V-flame behaviour.
Building from the analysis of linear flame-tip mode, we carry out an analysis of the non-local interaction in the subcritical oscillation in the nonlinear regime. A Monte Carlo simulation of passive flow tracers is conducted, with tracers departing from the burner exit at various phase time and radial positions. Based on the hypothesis that advected vortex structures interact with the flame surface, we design criteria to test this hypothesis by selecting the conforming trajectories and estimating the nonlinear oscillation frequency based on the mean values of the convective time delay corresponding to each conforming trajectory. In agreement with our physical reasoning, the estimated frequency is found to be close to the reference nonlinear oscillation frequency when similarly stringent selection thresholds are employed simultaneously for the vorticity and heat release rate. The resulting trajectories intersect with the flame surface at phase instances where the flapping flame tip is significantly further upstream than in the steady base states.
Finally, we assess the capacity of linear methods to predict basic features of the V-flame dynamics, as is commonly attempted in the reduced-order modelling literature. As neither modal nor non-modal analysis of the steady base flow provides any hint about the observed subcritical limit-cycle oscillations, we attempt an eigenmode analysis of the time-averaged mean flow. This approach is hampered by the strong nonlinearity of the system, particularly visible in the reaction rate term. Both a priori assessment and the computed mean-flow eigenmodes reveal notable disparities in the mean-flow eigenvalues and eigenmodes when compared with the reference Fourier modes associated with the nonlinear frequency of the limit-cycle fundamental. This result highlights known limitations of linear mean-flow analysis, namely its failure to provide valid predictions for subcritical instabilities and in the presence of strong non-monochromatic interactions. Such issues may indeed be quite common in premixed laminar flame oscillations driven by flame–flow feedback due to the well-known interaction between convective velocity perturbations at a fundamental tone and the unsteady heat release rate at the first harmonic (Lang Reference Lang1991; Lieuwen Reference Lieuwen2005). The present failure of mean-flow analysis in the V-flame configuration may be contrasted with our successful earlier efforts in cylinder-stabilised flames (Wang et al. Reference Wang, Lesshafft and Oberleithner2022b), which featured supercritical bifurcation behaviour dominated by monochromatic hydrodynamic amplification mechanisms localised to the wake region without significant feedback from downstream heat release fluctuations.
Acknowledgements
This study was inspired by the chapter ‘Modal analysis of a V-shaped flame’ of C.-H.W.'s PhD thesis. For the original work in that chapter, C.-H.W. and L.L. worked in collaboration with G. Varillon and W. Polifke (TU Munich), who computed base states and conducted time-resolved simulation of V-flames with OpenFOAM. The authors would like to express their sincere gratitude for the scientific inputs from G. Varillon and W. Polifke in the original work that inspired the present study.
Funding
C.-H.W. was supported through a PhD scholarship from École Polytechnique and by the Shuimu postdoc fellowship from Tsinghua University. C.-H.W. acknowledges the China Postdoctoral Science Foundation (CPSF) and the National Natural Science Foundation of China (NSFC) for funding this work (Grant No. 2022TQ0181, 12388101). C.D. acknowledges funding for this project from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska–Curie Grant Agreement No. 899987. Y.G. was supported by the PolyU Start-up Fund (Project No. P0043562) and the NSFC (Grant No. 52306166).
Declaration of interests
The authors report no conflict of interest.
