Hopf Bifurcation Analysis of the Combustion Instability in a Liquid Rocket Engine

: The bifurcation process of self-sustained combustion instability pressure perturbations in a liquid rocket combustor is investigated based on the Helmholtz equations and a pressure dependent ﬂame describing function. The modal frequency and growth rates are numerically resolved by the commercial software COMSOL multiphysics. Validation of the numerical approach is ﬁrstly conducted on a Rijke tube combustor, and a supercritical bifurcation for the ﬁrst longitudinal mode is observed. The bifurcation diagrams for the ﬁrst transverse mode for different time delays and gain index of the ﬂame describing function are analyzed. Only the supercritical bifurcation presents for this conﬁguration. The trajectory of Hopf points and the bifurcation diagram feature period motions with increasing the time delay. The effect of ﬂame length distributions on the bifurcation diagrams is analyzed by considering a non-uniform ﬂame length distribution model. Results show that the distribution has a large impact on the bifurcation process, e.g., the ﬁrst transverse mode is more unstable for the non-uniform distribution. Finally, a subcritical bifurcation is found when a more complicated ﬂame describing function is considered; the bistable region presents and the condition for this is discussed. G.W. and J.L.; data curation, X.L.; writing—original draft preparation, X.L.; writing—review and editing, L.Y., G.W. and J.L.; visualization, X.L.; supervision, L.Y.; project administration, L.Y.; funding acquisition, L.Y. and J.L. All authors have read and agreed to the published version of the manuscript.


Introduction
Combustion instability (CI) has been considered as one of the most difficult and challenging problems in aeronautics [1], astronautics [2,3] and other burners. Due to the severe operating conditions in the combustion chamber (high temperature, extremely high pressure and large energy density), CI more easily to occurs in liquid rocket engines (LRE) [4]. CI may lead to severe pressure oscillation in the combustion chamber, resulting in damages to injectors, the combustion chamber or other components and even causing an explosion. However, due to the highly complex multi-process coupling among injecting, vaporizing, combustion, acoustic characteristic, boundary condition and other factors, it is still much difficult to predict and eliminate the CI [5][6][7][8][9].
In 1859, Rijke found that unstable acoustic oscillations occured when the heat source was placed in a vertical tube, which provided a simple configuration to study the complicated CI phenomenon [10]. Based on the Rijke tube, many complex CI phenomena presenting in real geometries could be reproduced [11][12][13]. These studies could be extended to more complicated models by thoughtful deduction and appropriate extension [12,14,15]. Sujith et al. investigated the non-normality and nonlinearity in the CI interaction in a Rijke tube [16]. Bifurcation conduct has been identified in Rijke tube experimentally [17,18]. An online monitoring and optimization algorithm was developed based on the least mean square method to mitigate the CI in a Rijke tube [19,20]. Although these simplified models have captured the major features of the CI and facilitate the understanding of CI [11,21], the structure of it is relatively simple compared to full-scale combustion devices. At the same time, the flow in Rijke tubes is essentially regarded as laminar condition, which would not be sufficient to explain the dynamics of practical combustors with turbulent flow. Moreover, since only longitudinal modes are suited for studying in Rijke tubes, they still cannot reveal more factors of the true combustion chambers.
The BKD model combustor, which is a reduced-scale rocket engine developed at the DLR Lampoldshausen, has been considered as a benchmark test rig for the investigations of LRE combustion instabilities, and a great many researches have been conducted based on this configuration [22][23][24]. The linearized Euler code PIANO considering the mean flow and response function was applied to perform a complete stability analysis [25]. Gröning et al. [26] analysed the effect of hydrogen temperature on the high frequency CI in the BKD combustor. The compressible LES study of self sustained transversal CI and flame response subjected to transversal perturbations in the BKD were carried out by Urbano and her colleagues [27,28], and more practical results, e.g., mode triggering, were obtained. Machine learning has been applied for thermoacoustic instabilities in the BKD combustor, performance evaluated [29].
Due to the nonlinear responses of the flame and even the acoustic waves (these typically only occur in rocket engines), the CI features more complex behaviours. During past several decades, many researches have been conducted to develop linear and nonlinear theories [30][31][32][33]. Linear analysis is simple but is limited. Only linear stability can be obtained according to the linear analysis, and temporary status and the developing process of fluctuation typically cannot be determined [12]. Dowling has developed and summarized the low order models for the predictions of the CI in a simple duct combustor [11] and annular combustion chambers [14]. In order to capture the mode triggering and the limit cycle characteristics of the CI, nonlinear analysis is needed and a great many researches have been conducted to develop the nonlinear theories. Lieuwen studied the limit cycle phenomenon in a gas turbine combustor experimentally [34]. A typical self-exited CI of lean premixed combustion was examined by Sun [35], and the nonlinear bifurcation was captured. A comprehensive nonlinear combustion instability model has been proposed by Flandro [36]. In order to obtain the bifurcation diagram, systematic variation of parameters and tracking direct time integration were conducted by [37,38]. Besides, the nonlinear behavior of the solid rocket motor was analyzed to investigate the bistable region for the different responses [39].
As in many combustors, the flame is the major source of nonlinearity. Thus many researches have been conducted to develop the nonlinear flame response model, and it is also called the flame describing function (FDF) for the weakly nonlinear analysis. A lot of researches have worked out the nonlinear flame response for certain simple flames, e.g., the laminar premixed conical flame, based on the G-equation model [33,[40][41][42][43][44]. These models reveal the nonlinear dynamic flame responses and are much closer to real flames. However, they are limited to simple flame geometries, and are not suitable for highly complex turbulent flames. In the Rijke tube fed by a hot wire, the King's law is widely used to relate the nonlinear heat and the oncoming flow velocities, and this model has been coupled in the Galerkin theory to capture the nonlinear properties and even the subcritical and supercritical bifurcation [45]. The heat source of course can also be replaced by a laminar premixed conical flame model based on the G-equation method, and the results are more meaningful to the combustion research [13,46].
Other researches try to put forward simple nonlinear models which capture the main features of flame subjected to oncoming disturbances. Dowling has proposed a saturation model of the normalised heat release rate based on the fact that the heat release rate perturbation cannot be larger than the time-averaged heat release rate [31]. This model has been successfully used to capture the nonlinear interaction between longitudinal modes and circumferential modes in an annular combustion chamber [47], and has been further used to capture the nonlinear saturation in a Rijke tube combustor [48]. Li and Morgans had improved this model, considering a nonlinear factor contributing to both amplitude and phase of the FDF, and applied it to reconstruct different nonlinear CI behaviours in a Rijke tube [12]. Other researchers established polynomial relations between the heat release rate disturbance and oncoming perturbation, which is much easier to be used in the low order time domain simulation scheme [49,50]. Noiray has coupled a three order flame response model as a function of the pressure perturbation into a Galerkin approach to analyse the nonlinear circumferential mode in gas turbines [49]. Campa and Juniper have coupled the three orders and five orders flame models in the low order network code-LOTAN, to simulate the nonlinear bifurcations [51]. The nature of the Hopf bifurcation was analysed for different nonlinear flame models. These two models have been further coupled in the Helmholtz equations to simulate the nonlinear CI in more complex geometries [52,53], which used the usual velocity dependent flame transfer function. As for most combustors, the major nonlinearities come from the nonlinear flame responses; this facilitates the analysis of combustion instabilities based on the Helmholtz equation, with the nonlinear flame response as the source.
The objective of the present work is to couple the nonlinear polynomial pressure dependent flame describing function into a Helmholtz equation model of the BKD combustor to study the Hopf bifurcation of combustion instability in the LRE. The combustion instability in a Rijke tube was studied firstly and then that of the BKD combustor with different distributed flame lengths. The remainder of the paper is organized as follows. The simulation models, geometry models (Rijke tube and BKD combustor), nonlinear flame model and distributed flame lengths model are introduced in Section 2. The nonlinear analysis method was introduced in Section 3. Section 4 presents the results and the discussion of the bifurcation diagram. Finally, the conclusion is shown in Section 5.

Governing Equations and the Model of Flame Describing Function
The fluid in the target combustion chamber is considered as ideal gas. The effects of viscosity, heat diffusion and heat conduction can be neglected. The pressure and the velocity within the chamber are supposed to be uniform. Because the mean flow velocity is far smaller than the speed of sound, it is reasonable to assume a zero mean flow, even for LRE combustor [54,55]. Under these assumptions and the small perturbation assumption, it is easy to get the inhomogeneous wave equation as a function of the pressure perturbation p [11,14]: where q represents the heat release rate per unit volume and c is the speed of sound. ρ and γ denote the density and the specific heat ratio, respectively.( ) means an average value and ( ) indicates a perturbation value. With the assumption of perturbations on the basic frequency, a (t) =â exp(iωt), where a could be p or q. ω is the complex angular frequency; its real part R(ω) = 2π f represents the angular frequency, and its imaginary part I(ω) = −α corresponds to the modal growth rate. When α > 0, the acoustic mode is unstable, and the disturbances grow with time. When α < 0, the acoustic mode is stable, and the disturbances decay with time. When α = 0, the acoustic mode is marginal stable or a limit cycle for a perturbation is established. By substituting the above expression into it, the wave equation can be changed to the expression in the frequency domain: A finite element method (FEM) based the commercial software COMSOL multiphysics is used to solve the Helmholtz equation. This numerical scheme is based on the ARPACK numerical routine for large-scale eigenvalue problems [56]. The numerical procedure is iterated until the error is smaller than 10 −4 Hz. In previous researches, the flame response to oncoming disturbances is typically based on the relation between the heat release rate and the axial flow velocity perturbations [53]. However, in the rocket combustion chamber, transverse CI typically occurs and the model based on the axial velocity perturbation is limited. One thus prefers the flame response model to local pressure perturbations.
In order to account for the nonlinearity of the heat release rate, a flame response model adopted the classical n − τ model [57] and the nonlinearity is introduced by a third order term as used in [51,53], and this model can be expressed in the time domain as: where n and τ are the gain index and the time delay of the n − τ model, respectively. It is noted that the presence of the third order term leads to the nonlinear flame response. µ is a coefficient related to the strength of the third order term compared the first order term. In the present paper, the coefficient µ equals 2 for the nonlinear analyses and equals 0 for the linear analyses. In order to couple the flame response model in the frequency domain analysis, a weakly nonlinear flame model, also called the FDF, can be obtained by accounting for the Laplace transform of the fundamental term of heat release rate perturbation from Equation (3), and the FDF can be expressed as: where β = |p/p|. It is noted that, when normalised pressure perturbation amplitude β remains small, the third order term is much smaller than the first order term, and the flame describing function becomes a linear relation; with increasing β, the gain of the flame describing function obviously decreases, that captures the saturation behaviour of the flame dynamic response [58]. As the time delays in the first order and the third order terms are the same, the gain of the FDF decreases monotonously [50]. The nonlinearity term of a fifth order polynomial was similarly introduced, where µ 3 and µ 5 are the strength of the third order term and the fifth order compared to the first order term, respectively. In the present paper, the coefficient of the third order µ 3 equals −5 and the coefficient of the fifth order 6. The corresponding FDF was obtained according to the procedure of the third FDF: at the moment, the amplitude of the pressure disturbance increases with the disturbance developing, resulting in the non-monotonic change of the nonlinear term, which causes more complicated nonlinear conduct. By substituting Equations (4) and (6) into the Helmholtz Equation (2), one gets ODE equations as purely functions ofp, and this equation thus can now be numerically solved.

Rijke Tube Model Combustor
Following previous research [53], the proposed models and solver have been firstly applied to the simplest geometry-Rijke tube model combustor, as shown in Figure 1. The length of the tube is 3 m and the radius of the tube is 0.1 m. Both ends of the tube considered open to the atmosphere, i.e., p = 0. The flame is placed at one quarter of the tube (as represented by the red block in Figure 1), and has a width of 0.02 m. The mean temperature before and after the flame aret 1 = 300 K andt 2 = 700 K, respectively. The flow in the two sections is assumed as air, and the thermal properties of air are used in the calculation. The specific ratio is assumed constant and equals γ = 1.4. Figure 1 also shows the mesh for the simulation.

BKD Combustor Model
Further analyses are conducted on the reduced-scale rocket engine model BKD combustor. The combustion chamber of BKD is 8 cm in diameter and 20 cm in length, fed by 42 injectors. There is an injection system before the chamber and a convergent-divergent nozzle after. More parameters are described in [59]. The first transverse (1T) mode showed an unstable mode of this combustor under case LP4, presented in [60]. The pressure of combustion chamber is 80.04 bar. The operation parameters are shown as Table 1 in terms of propellant mass flow ratesṁ O 2 /H 2 , ratio of oxidizer to fuel mass flow (ROF=ṁ O 2 /ṁ H 2 ), velocity at injector exit u O 2 /H 2 ,inj and the average pressure of the combustion chamberp cc . The acoustic model of BKD is built upon its real geometry, whose main body is a cylinder. The FEM model is presented in Figure 2. Walls of the injectors and combustion chamber are treated as rigid walls. Hydrogen tubes were modeled by the impedance proposed by Urbano [27] (Z = −1.160 − i0.255 for the first order transverse mode) instead of modeling true hydrogen tubes to simplify the simulation. For the same purpose, the nozzle is replaced by an acoustic boundary condition, whose pressure reflection coefficient was calculated by the Magnus expansion of the governing linearised Euler equations of the nozzle [61,62]. The resulting acoustic reflection coefficient is shown in Figure 3, where Ω means the reduced frequency and Ω = f L n /c 0 , L n is the nozzle length, c 0 is the speed of sound at the inlet of the nozzle. It is obvious that the reflection coefficient features the shape of a low pass filter. It is noted that the mean properties vary along the combustion chamber. In order to resolve the Helmholtz equation (Equation (2)), the mean properties, especiallyρ andc, are obtained by a RANS (Reynolds Average Navier-Stokes) of a single injector connected to a combustion chamber with an equivalent radius [63]. Figure 4 shows the mean properties along the axis from RANS and thermal calculation software REFPROP considering the supercritical state. The mean temperature arises with increasing the axial location and reaches around 3000 K when x = 0.2 m. The mean speed of sound also increases along the combustion chamber. Due to the presence of hydrogen, the speed of sound exceeds 1500 m/s when x = 0.2 m. Based on the numerical simulation results [27], the flame region is assumed to be extended to the half of the combustion chamber to construct the combustor model and is represented by the red block in Figure 2. Furthermore, in order to study the effect of the distribution of flame length along the radius direction on the CI, three kinds of flame length distributions are accounted for, and are presented in the following section.

Models of Flame Distribution in the BKD Combustor
However, in the true liquid rocket engine, the flame distribution is usually not a standard cylinder, e.g., the simulation results of the BKD combustor as shown in Figure 11b in [27]. Previous researches [64,65] about real liquid rocket engines also showed that the flame region (or high temperature region) in the combustion chamber was distributed along the radius direction. One thus also accounts for another two flame distribution models, a convex shape and a concave shape. The variation of flame length along the radius obeys a Gaussian distribution, and the flame length distributions for these two models can be expressed as: where L f ,+ (r) and L f ,− (r) represent the flame lengths as a function of radius r for the convex and concave shapes, respectively. L cc represents the length of the entire combustion chamber. σ equals 1.58 herein, a and b equal 5.7737 and 1.1934, respectively. The difference between the flame length at the center r = 0 and that near the lateral wall r = R is 0.05 m. The uniform flame in the BKD combustor is considered half of the length of the combustion chamber, as shown in Figure 2. The flame length distribution for the three types of envisaged flames are shown as Figure 5. The variables are nondimensionalized to describe the shape of the heat resource. The effect of flame length distribution on the nonlinear CI will be described in details in Section 4.3. Figure 5. Flame lengths as a function of the radius r for the three kinds of distribution shapes.

Nonlinear Analysis
The objective of this work is to investigate the Hopf bifurcation of thermoacoustic mode when gradually increasing or decreasing a parameter, e.g., the index n, of the flame describing function. This can be achieved by looking into the growth rate α path of the envisaged mode with increasing the normalised pressure perturbation amplitude β for a given value of n for instance [53]. By searching the value that α = 0, one can get the normalised pressure perturbation amplitude β * when the limit cycle is established. It should be noted that, for certain values of n, the system remains stable for all β, and one assigns 0 to β * . It is then possible to consider the value of n when β * changes from 0 to a positive value as the Hopf point or stability margin, which means that the system changes from linear stable to unstable. The bifurcation diagram can then be established by plotting β * versus n. Note that the above procedure can be simplified as the gain of FDF decreases monotonously with increasing β. The mode thus is typically mostly unstable when β = 0. Thus, one can firstly calculate the linear stability for all n and search the Hopf point when the modal growth rate α = 0. Considering more complex FDF, e.g., the fifth order FDF as used in [53], more complex nonlinear combustion instability behaviours, e.g., the hysteresis phenomenon, occur with increasing β. The bifurcation diagram should be constructed by carefully scanning the combustion instability properties for all n and β. In order to get the return path of the bifurcation diagram, the decreasing process of n should be conducted.

Rijke Tube
To get the eigenfrequency of Rijke tube, the Helmholtz equation was solved by the commercial software COMSOL Multiphysics. It should be noted that the heat release rate q in the FDF equation (Equation (4)) is the local value. The local mean heat release rate equals q =Q/V f , whereQ and V f are the entire mean heat release rate and the volume of the flame. For the Rijke tube,Q =ρ 1ū1 c p (t 2 −t 1 ) = γū 1 /(γ − 1)p 1 (t 2 /t 1 − 1), wheret,ρ 1 ,ū 1 and c p are the mean temperature, mean flow density, velocity and heat capacity at constant pressure, respectively. It is assumed a low Mach number flow in the Rijke tube and the inlet mean flow velocity equalsū 1 = 1 m/s. The mean pressure thus can be considered constant along the duct. The local heat release rate perturbation for this configuration thus becomes: The fundamental frequency of pressure perturbation in the Rijke tube is 74 Hz when there is no heat source (i.e., n = 0). This mode shape is presented in Figure 6a, with low pressure oscillation amplitudes near two ends and high amplitudes in the middle. The bifurcation diagram shown in Figure 6b is plotted under the condition that the time delay of the FDF τ is 7 ms, and a supercritical bifurcation presents. These results qualitatively match those results shown in [53], validating the methodology used in this paper. The value of n is increased by a step of 0.003. When n < 0.2933, the thermoacoustic system remains stable, and β * = 0. n H = 0.2933 thus is considered as the Hopf point. When n exceeds 0.2933, the envisaged acoustic mode becomes unstable, and a limit cycle is established when β = β * . It is noted that β * increases monotonously. Furthermore, the limit cycle value β * increases rapidly when n is close to the Hopf point, and then the increased speed decreases with increasing n.

BKD Combustor with Uniform Flame 0 Distribution Along the Radius
Due to the more complicated geometries and flow field in the BKD combustor, the FEM calculation becomes more complex. The grid number used for the FEM calculation is 112,008. The grid in the injection system is dense due to the small radius, while that in the combustion chamber remains coarse due to the sufficiently large acoustic wavelength. For this configuration, the global mean heat release rate is calculated by the RANS simulation; the local mean heat release rate is assumed to equalQ/V f , where V f again represents the volume of the flame. Thus the local heat release rate perturbationq equalŝ p/p ccQ /V f F(ω, β). The Helmholtz equation was solved by the COMSOL Multiphysics to get the modal frequencies and growth rates, and further the mode shapes.
In this work, one firstly considers the situation that the flame length distribution along the radius is uniform. The frequencies of the first longitudinal and first transverse thermoacoustic modes equal 2919.3 Hz and 9096.2 Hz respectively, and their mode shapes can be viewed in Figure 7a,b. The frequency of the first transverse mode 9096.2 Hz is similar to that of experiment 10,260 Hz [27], the offset may be because of physical parameter error from RANS. It is noted that the pressure perturbations in the oxygen injection unit are both smaller than those in the combustion chamber for the two modes. For the first transverse mode, a maximum pressure perturbation is attained near the injection plate. Under these conditions, the combustor is more prone to combustion instabilities. For the first longitudinal mode, the combustion instability in the combustor behaves more like that in the Rijke tube. Furthermore, the transverse combustion instabilities are more prone to occur in the LRE combustor and are more harmful. One thus pays more attention to the first transverse mode. One firstly considers the linear stabilities of the first transverse mode, i.e., only the linear part of the FDF is considered and µ = 0 in Equation (4). For the sake of simplicity, one defines a time period T 1T = 1/( f 1T ), where f 1T = 9096.2 Hz. It is then possible to normalise the time delay of the FDF by this time period, as it is typically found that the modal frequency and growth rate change periodically and the period typically approximately equals the inverse of the modal frequency [12,53,66]. Figure 8 shows the linear results of the first transverse mode for different time delays τ. For a small time delay τ/T 1T = 0.1, the first transverse mode remains stable and the growth rate α < 0 when n < 0.08 . With increasing n, the growth rate α increases and equals 0 when n = n H . One considers n H as the index value for the marginal stable situation. When n is further increased, this mode becomes unstable, and the nonlinear term should be considered. It is then possible to increase the time delay τ, and results show that the trajectory of n H also features a periodic motion. n H reaches the maximum value of 6.787 when τ/T 1T = 1 and has the minimum value of 0.005 when τ/T 1T = 0.5. It should be noted that the mode is more stable for a larger value of n H , i.e., the first transverse mode is more stable when τ/T 1T approaches 1 and more unstable when τ/T 1T approaches 0.5. It is also interesting to note that, when τ/T 1T ≈ 1, a small change of the time delay may lead to a large change of n H , indicating that the stability of the first transverse mode is more sensitive to τ when τ/T 1T ≈ 1 and is more inert when τ when τ/T 1T ≈ 0.5. It is also possible to examine the change rate of the growth rate α with n when n = n H . The trajectory of this value dα/dn with τ is also shown in Figure 8. When τ/T 1T ≈ 0.5, the value of dα/dn reaches the maximum value of 85 s −1 , indicting that the stability is more sensitive to n. When τ/T 1T ≈ 1, dα/dn has the minimum value, and the modal stability is less sensitive to n. One now considers the nonlinear behaviours of the first transverse mode by accounting for the full expression of the FDF (Equation (4)). Figure 9 shows the bifurcation diagrams of the first transverse mode of the BKD combustor for different time delays τ. One firstly takes the diagram when τ/T 1T = 0.2 to illustrate the bifurcation. When n is smaller than the value of Hopf point n H , the mode is stable. With increasing the value of n, the mode becomes more unstable. When n exceeds n H , the mode becomes unstable, and a limit cycle is established when β = β * , where the growth rate again equals 0 and a limit cycle is established. It is noted that the value β * arises with increasing n. Similar to that for the Rijke tube, a supercritical bifurcation also presents. When τ changes, the change in the shape of the bifurcation diagram remains small while the location changes. Since the calculation is quasi-symmetric with the symmetric axis τ/T 1T = 0.5, a period motion of the entire diagram with τ is also observed, which agrees with the linear analysis. It is further noted that, when τ/T 1T is close to 0 or 1, the Hopf point n H is larger and the value of β * is smaller. When τ/T 1T is close to 0.5, the Hopf point n H is smaller and the β * increases more rapidly when n is close to n H , which agrees with the results for dα/dn in Figure 8.
It is also interesting to note that the 3D diagram shown in Figure 9 degenerates to the linear stability results shown in Figure 8 when β = 0. It is obvious that linear analysis is helpful even in the nonlinear studies. When using the nonlinear flame model, all disturbances finally grow or decay to zero or the values corresponding to the limit cycle state. When the nonlinear flame model is close to real dynamics flames, if possible, and the resolution of the data is high enough, the state that an arbitrary disturbance will finally reach could be predicted precisely. It is also possible to examine the relation between τ and β * or the τ − β * bifurcation diagram based on the τ − β * plane for a fixed n from the 3D diagram, as also shown in [67].

Effects of Flame 0 Distribution on the Bifurcation Diagram of the BKD Combustor
In a supercritical bifurcation system like that shown in Figure 9, there is a consistent one-to-one match between a pair of τ − n and one value of β * ; the three dimensional bifurcation diagrams of the first transverse mode for the BKD combustor thus could be shown as contour maps as shown in Figure 10. The differences among the contours for the three flame length distributions are obvious. When τ/T 1T equals 0.5, the limit cycle for first transverse mode is established earlier for the non-uniform distributions than that for the uniform distribution. For the concave distribution, the marginal stable border is flatter, indicating that the marginal stability is more insensitive to the change of τ. At the same time, the more flat concave distribution illustrates that the system would become unstable when the gain index n is low. Maybe under this situation, the heat release rate and large oscillation region of pressure disturbance are in phase well, having more volume superimposed region. Furthermore, it is noted that the first transverse mode for the BKD combustor with a non-uniform flame length distribution is more likely to attain a limit cycle with a larger perturbation amplitude.
In order to quantitatively compare the differences among results for the three flame length distributions, β * − n plots under different τ are extracted and shown in Figure 11. Due to the quasi-symmetric shape of the contour as a function of n (Figure 10), only results for half part are presented. Results show that the first transverse mode is more unstable for non-uniform flame length distributions. For example, when τ/T equals 0.1, the first transverse mode of the BKD combustor with concave and convex shapes enter the limit cycle when β equals 0.018 and 0.028, and the value of β * is 0.0455 for the uniform distribution. Furthermore, the limit cycle value β * increases with increasing τ and attains the maximum value when τ/T 1T is close to 0.5. This also agrees with the periodic motion as illustrated in Figure 8. Further, the concave, convex and uniform flame length distributions represent different flow rate distributions. Once the law was commanded, combustion instability could be attempted to eliminate by designing flow rate distribution. Therefore, the simulation of this paper suggested that a uniform flame shape may be a good way to control combustion instability.

Bistable System
To investigate more complicated nonlinear behaviours of thermoacoustic instabilities, the fifth order flame describing function was accounted for in the BKD combustor with the uniform flame distribution. The phenomena of subcritical bifurcation were captured and plotted in Figure 12. Solid lines represent the stable states or limit cycles, while dashed lines indicate unstable limit cycles. Four time delays were accounted for. The green and black arrows indicate the route when the flame describing function gain index n increases and decreases, respectively. The Hopf points and the fold points were indicated by red stars and light blue triangles. Under the situation of subcritical bifurcation, these points are called subcritical Hopf bifurcation points and fold bifurcation points. When the flame describing function gain index n is between the values for the Hopf point n H and fold point n f , the system is linear stable but nonlinear unstable, indicating that the system remains stable if the pressure perturbation is smaller than that for the unstable limit cycle, while the system will increase to a stable limit cycle if the pressure perturbation is larger than the value for the unstable limit cycle. Figure 12. The bifurcation diagrams for the fifth order FDF. The green arrows represent the system developing direction when the gain n increases, the black arrows represent the system developing direction when the gain n decreases.
The region of stability analysis for the variation of flame describing function gain index n and the time delay τ is shown in Figure 13. The black solid line and the red solid line indicate the fold points and the Hopf points, respectively. Between them, different regions represent the different stable states. Similar to the contour of the pressure perturbation amplitude in the previous section, the lower part of Figure 13 is a globally stable region independent of the pressure perturbation amplitude. Similarly, the upper part is a globally unstable region independent of the pressure perturbation amplitude. Between them, a region that is neither globally stable nor globally unstable emerged, which considered as the bistable region. A stable steady state and a pair of stable and unstable limit cycles are included in this region. The converged state is determined by the initial pressure disturbance.

Conclusions
This paper analysed the nonlinear combustion instability, especially the Hopf bifurcation, of two combustors with a pressure dependent flame describing function based on the Helmholtz equation numerically resolved by the commercial software COMSOL multiphysics. The numerical method was firstly implemented on a Rijke tube combustor for validation, and a supercritical bifurcation for the first longitudinal mode was observed. The nonlinear responses and further the bifurcation diagrams of the first transverse mode for different time delay and gain index of the FDF were analysed. The change law between n H and τ was obtained. The trajectory of Hopf points and even the bifurcation diagram feature period motions with increasing the time delay τ of the FDF. For the envisaged flame model, the supercritical bifurcation and the subcritical bifurcation can be established. The effect of flame length distributions on the bifurcation diagrams was analysed by considering a non-uniform flame length distribution model. Results showed that the first mode is more easily to attain the limit cycle and is more insensitive to the time delay for the non-uniform distribution. Furthermore, the first transverse mode for the BKD combustor with a non-uniform flame length distribution is more likely to attain a limit cycle with a larger perturbation amplitude. Finally, a subcritical bifurcation is investigated and the bistable region, linear stable but nonlinear unstable, is analysed. Due to the large scales and multiple injectors in the liquid rocket engines, the direct simulation of the combustion instabilities within these devices is much more time consuming and impossible for the engineering design. The present work decouples the flame responses simulation from the entire simulation and can be used to reduce the computing cost. Although calculation of the flame responses have not been conducted, which is replaced by a nonlinear flame response model, these results in this work is worthy for the analysis of the nonlinear behaviors of the combustion instabilities, especially for the transverse modes. The discussion on the flame response distribution may be helpful for the injector's optimization to eliminate the combustion instability.

Acknowledgments:
We would like to express our thanks to the editors of Aerospace and the anonymous reviewers for their work in processing this article.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

FDF
Flame describing function CI Combustion instability RANS Reynolds average Navier-Stokes LRE Liquid Rocket engine