Discussing the limitations of unconditional stability indicators in evaluating thermoacoustic quality of burners with flames

Network modeling and dispersion relation

To apply one-dimensional linear acoustic network modeling in the frequency domain, the frequency range of interest should be smaller than the cut-off frequency of the first transverse mode ¹⁶ that is usually valid in thermoacoustics. There are several possibilities for the variables used to express the transfer matrices of the elements. Here let us focus on the Riemann invariants f j and g j where are written in terms of the complex acoustic pressure p j ′ and the acoustic particle velocity u j ′ . Where j = 1 , 2 , and the numbers 1 and 2 indicate upstream and downstream, respectively. At the limit of zero Mach number, we have

f j = 0.5 ( p j ′ z j + u j ′ ) , g j = 0.5 ( p j ′ z j − u j ′ )

(1)

[ f 2 g 2 ] = [ T 11 T 12 T 21 T 22 ] ⏞ [ T ] [ f 1 g 1 ]

(2)

[ g 1 f 2 ] = [ S 11 S 12 S 21 S 22 ] ⏞ [ S ] [ f 1 g 2 ] = [ − T 21 T 22 1 T 22 T 11 − T 21 T 12 T 22 T 12 T 22 ] [ f 1 g 2 ]

(3)

OPEN IN VIEWER

Figure 1.

Generic 1D acoustic network model of a combustion system; (a) MIMO layout; SISO representation, (b) in the input plane; and (c) in the output plane. MIMO: multiple-input-multiple-output; SISO: single-input-single-output; 1D: one-dimensional.

In thermoacoustics, the common situation is encountered when the heat release is concentrated in a region much smaller than the acoustic wavelength. Therefore, the conservation equations for the flame can be reduced to jump conditions known as Rankine-Hugoniot (R-H) relations. To obtain the relationship between the acoustic quantities up/down-stream of the flame, one needs a closure relationship. For laminar premixed flames, it is generally assumed that the heat release fluctuations are mainly due to the velocity fluctuations upstream of the flame. In this case, one can define a flame TF ( F ( ω ) ) as the normalized ratio of the heat release rate fluctuation Q ˙ ′ / Q ˙ ¯ to the normalized ration of acoustic velocity excitation u ′ / u ¯ at some reference location upstream of the burner.

F ( ω ) = Q ˙ ′ / Q ˙ ¯ u ′ / u ¯

(4)

T = [ T 11 T 12 T 21 T 22 ]

(5a)

T 11 = T 22 = 0.5 ( ϵ + 1 + θ F )

(5b)

T 12 = T 21 = 0.5 ( ϵ − 1 − θ F )

(5c)

[ − 1 R u p 0 0 T 11 T 12 − 1 0 T 21 T 22 0 − 1 0 0 R d n − 1 ] [ f 1 g 1 f 2 g 2 ] = 0

(6)

Π T = T 22 − T 12 R d n + T 21 R u p − T 11 R u p R d n = 0

(7)

Stability analysis in frequency domain

As the localization of the roots (zeros) of dispersion relation in s-plane is crucial to resolve the dilemma of the system stability-instability, one can use the corresponding analysis methods developed in the complex function analysis theory to approach this issue. In this section, we will begin by offering an extensive introduction to fundamental principles originating from the complex function theory. Afterward, we will explore how these principles can be practically applied to analyze the challenge of stability and simplify the task of stabilizing flames.

Cauchy’s argument principle ¹⁸

If Π T ( s ) is a generic holomorphic function and nonzero at each point of a simple closed contour C and is meromorphic inside C , then

Z − P = N = 1 2 π i ∮ C 1 Π T ( s ) d Π T ( s ) d s d s = − 1 2 π Ξ ( Π T )

(8)

where Z and P are accordingly the number of zeros and poles of Π T ( s ) evaluated in the domain encompassed by the contour C . N is known as the winding number of Π T ( C ) , and counts the number of turns that the curve Π T ( C ) winds clockwise around the origin. Ξ ( Π T ) indicates the unwrapped argument of Π T . Figure 2 shows a clarification example of Cauchy’s argument principle when Π T is evaluated along contour C , in which blue and red circles represent zeros and poles of Π T ( s ) , respectively. The number of encirclement can be easily calculated as N = 2 − 3 = − 1 .

OPEN IN VIEWER

Figure 2.

A clarification example of Cauchy’s argument principle. The blue and red circles represent zeros and poles of Π T ( s ) in the s -plane, respectively, and so the number of encirclements in Π T ( s ) -plane around the origin shown with a black star can be calculated as N = 2 − 3 = − 1 .

Thermoacoustic stability analysis using Cauchy’s argument principle

To assess the exponential stability of a linear thermoacoustic system, one must analyze the presence or absence of RHP zeros in the dispersion relation. The analysis based on the Cauchy’s argument principle enables the identification of RHP zeros and poles of the characteristic function Π T ( s ) . To perform the corresponding evaluation, as depicted in Figure 3, Cauchy’s argument principle is applied by selecting the contour C to enclose the entire RHP ( C = D + ). Accordingly, the function Π T ( s ) must to be known along the imaginary axis and at infinity. In many physical systems, the frequency response Π T ( i ω ) is sufficient for assessing stability, especially when evaluating a (semi-)proper target function with the property lim s → + ∞ Π T ( s ) = constant . In control theory, this approach is widely employed with the assistance of a graphical tool, namely the Nyquist plot. ¹⁹

OPEN IN VIEWER

Figure 3.

A clarification schematic of using Cauchy’s argument principle in right-half plane for stability analysis. The blue and red circles represent zeros and poles of Π T ( s ) , respectively.

It is important to highlight that stability can also be determined by directly analyzing the argument of the dispersion relation. Recently, multiple thermoacoustic stability criteria based on different forms of the dispersion relation, along with their corresponding conservative forms, are established using the principle of argument. ²⁰ Before applying the argument principle in RHP of the complex domain, it is good to make two remarks:

Remark 1. For a passive single-input-single-output (SISO) system, the TF has no RHP poles. Since we are usually dealing with acoustically passive terminations, we can conclude that R u p and R d n have no RHP poles.

Remark 2. Direct measurability of a SISO function is an indication of the absence of RHP poles. Therefore, the TF F and the coefficients of the transfer matrix T i j have no RHP poles.

Unconditional/inherent stability criteria for an active two-port device

In the design process of an active two-port device, achieving unconditional or inherent stability is a fundamental concept. It implies that the device will not oscillate under any combination of possible passive terminals. However, attaining unconditional stability in acoustics is a complex task. In this article, we aim to shed light on two widely used criteria for achieving unconditional stability: the Rollett conditions ( k ) and the Edwards-Sinsky ( μ ) factor. Additionally, we will discuss the limitations associated with these requirements. Let us begin with the mathematical background of these requirements. Equation (7) represents the dispersion relation expressed in terms of a multiple-input-multiple-output (MIMO) representation. It can also be formulated as a SISO representation of the system, as depicted in Figure 1(b). In order to adopt the current approach, we reorganize the terms of the equation as follows. By dividing the dispersion equation (equation (7)) by T 22 − R d n T 12 , and using the transformation relation between transfer matrix and scattering matrix we obtain the following expression:

Π T = Π I S ⋅ den [ R i n S ] ⋅ T 22 = 0

(9)

where

Π I S = 1 − R u p R i n S = 0 , R i n S = S 11 + S 21 S 12 R d n 1 − S 22 R d n

(10)

where the operator den [ ⋅ ] indicates denominator of function [ ⋅ ] , and variables Π I T and Π I S are named input dispersion relations based on the transfer matrix and scattering matrix, respectively. The terms R i n T and R i n S can be referred to as the input reflection coefficients based on the transfer matrix and scattering matrix, respectively, which are defined as R i n T = R i n S = g 1 f 1 (see Figure 1(b)). The SISO form of the dispersion equation is widely used in microwave theory,^13,21,22 and has already been applied to thermoacoustic problems.^23–26

Discussion I. From a broad perspective, the resemblance between microwave theory and thermoacoustics opens up possibilities for utilizing the extensive formalism of systems theory to establish a rigorous framework for studying thermoacoustic instabilities. However, there exists a fundamental difference in the approach employed by the two fields. In microwave theory, scattering matrices are directly measured, and as highlighted in Remark 2, it can be ensured that the scattering coefficients do not possess any RHP poles. In other words, according to the definition in equation (3), T 22 does not have any RHP zeros. Conversely, in thermoacoustics, the zeros of T 22 are referred to as intrinsic thermoacoustic (ITA) modes. Consequently, it can be concluded that the present notation is derived under the assumption of the absence of unstable intrinsic modes.

Input (plane) definitive stability (IDS) criterion based on Π T = Π I S ⋅ d e n [ R i n S ] ⋅ T 22 = 0

²⁰

Using Cauchy’s argument principle in the RHP, we can rewrite equation (8) for the given dispersion relation in equation (9) as follows:

N Π T = Z Π T − P Π T

(11)

Z Π T and P Π T can be written as follows:

Z Π T = Z Π I S + Z d e n [ R i n S ] + Z T 22

(12)

P Π T = P Π I S + P d e n [ R i n S ] + P T 22

(13)

This means that to evaluate Π T , one has to separately evaluate Π I S , d e n [ R i n S ] and T 22 . One can think of two scenarios here.

First scenario: Z T 22 ≠ 0 and/or Z d e n [ R i n S ] ≠ 0 : This condition can be recognized as an artificial pole-zero cancellation case. A thermoacoustic system calls exponentially stable when Z Π T = 0 . If one substitutes this requirement in equations (11) and (12):

N Π T = − P Π T

(14)

Z Π I S = − Z d e n [ R i n S ] − Z T 22

(15)

The RHP poles of the T 22 is equal to RHP poles of the flame TF, so P T 22 = P F = 0 ; also with respect to the definition of R i n S , one can write P d e n [ R i n S ] = Z T 22 and P R i n S = Z d e n [ R i n S ] . Keeping in mind these equalities, after substituting equation (13) in equation (11):

N Π T = − Z d e n [ R i n S ] − Z T 22

(16)

Equations (15) and (16) can be summarized as follows:

N Π T = Z Π I S

(17)

Equation (17) implies that even in cases of pole-zero cancellation with 100% fidelity, one still needs to evaluate the encirclement of the basic form of the dispersion relation given in equation (7) to provide information on the stability dilemma. Therefore, struggling with the dispersion relation in the form of the input plane as given in equations (9) and (10) is not useful. Another key conclusion is that the absence of RHP zeros of the input dispersion relation ( Z Π I S = 0 ) does not necessarily indicate anything regarding the stability dilemma. We will return to this point later in this article.

Second scenario: Z T 22 = 0 and Z d e n [ R i n S ] = 0 : Using Cauchy’s argument principle in RHP, recalling equations (15) and (16) and substituting Z T 22 = 0 and Z d e n [ R i n S ] = 0 , the necessary and sufficient conditions for MIMO stability is the absence of encirclement for functions Π I S and d e n [ R i n S ] and it can be summarizes as follows:

N Π I S = N d e n [ R i n S ] = N T 22 = 0

(18)

In microwave theory N T 22 equals to zero due to the direct measurement of the scattering matrix. In the thermoacousics, in the case of direct assessment of scattering matrix (not flame TF) can be the case.

Input (plane) conservative stability (ICS) utilizing the concept of infinite phase margin ¹³

In the previous subsection, we derived sufficient conditions (though not necessary) to ensure the stability of a generic active two-port device in the absence of an unstable intrinsic mode, as described in equation (18). Now, the objective is here to establish a conservative stability criterion based on the concept of infinite phase margin. A system is said to have infinite phase margin if the magnitude of its open-loop SISO TF is always less than or equal to unity. In accordance with equation (18), we need to apply this concept to two components: the denominator of the input reflection coefficient ( den [ R i n S ] = 1 − L 1 , where L 1 = S 22 R d n ) and the input dispersion relation ( Π I S = 1 − L 2 , where L 2 = R u p R i n S ):

| L 1 | < 1 ⇒ | S 22 | < 1 | R d n |

(19)

| L 2 | < 1 ⇒ | R i n S R u p | < 1

(20)

Starting from equation (20) and substituting equation (10), we can rewrite it as follows:

| R i n S R u p | = | S 11 R u p + S 21 S 12 R u p R d n 1 − S 22 R d n | = | S 11 ′ − Δ ′ R d n | | 1 − S 22 R d n | < 1

(21)

where

S 11 ′ = S 11 R u p , S 21 ′ = S 21 R u p

(22a)

Δ = S 11 S 22 − S 12 S 21 , Δ ′ = Δ R u p

(22b)

If we set | R i n R u p | = 1 , we establish a boundary of R d n in the input plane that separates stable from unstable regions:

| S 11 ′ − Δ ′ R d n | = | 1 − S 22 R d n |

(23)

We now proceed by squaring both sides of equation (23) and collecting terms, resulting in

| S 11 ′ | 2 + | R d n | 2 ( | Δ ′ | 2 − | S 22 | 2 ) + 2 ℜ ( R d n ) ℜ ( ( S 22 − Δ ′ S 11 ′ * ) * ) + 2 ℑ ( R d n ) ℑ ( ( S 22 − Δ ′ S 11 ′ * ) * ) = 1

(24)

Which simplifies to:

| S 11 ′ | 2 + | R d n | 2 ( | Δ ′ | 2 − | S 22 | 2 ) + R d n ( S 22 − Δ ′ S 11 ′ * ) * + R d n * ( S 22 − Δ ′ S 11 ′ * ) = 1

(25)

Adding | S 22 − Δ ′ S 11 ′ * | 2 / ( | S 22 | 2 − | Δ | 2 ) 2 to both sides of equation (25) results in

| R d n | 2 + | S 22 − Δ ′ S 11 ′ * | 2 ( | S 22 | 2 − | Δ ′ | 2 ) 2 − R d n ( S 22 − Δ ′ S 11 ′ * ) * ( | S 22 | 2 − | Δ ′ | 2 ) 2 − R d n * ( S 22 − Δ ′ S 11 ′ * ) ( | S 22 | 2 − | Δ ′ | 2 ) 2 = ( | S 11 ′ | 2 − 1 ) ( | S 22 | 2 − | Δ ′ | 2 ) 2 + | S 22 − Δ ′ S 11 ′ * | 2 ( | S 22 | 2 − | Δ ′ | 2 ) 2

(26)

Equation (25) is of the form of a circle in the input plane for R d n , that is

| R d n − C i | 2 = | R i | 2

(27a)

where

C i = ( S 22 − Δ ′ S 11 ′ * ) * | S 22 | 2 − | Δ ′ | 2 , R i = | S 21 ′ S 12 | | S 22 | 2 − | Δ ′ | 2

(27b)

are the center and radius of this circle, respectively. The stability circles represent the boundaries of permitted terminations in their respective reflection coefficient planes. To find whether the interior or exterior of the stability circles is representing the stable region, we set R d n = 0 , which results in | R u p R i n | = | S 11 ′ | . So, we can, therefore, say that, if | S 11 ′ | < 1 , the origin of the downstream plane in the Smith chart (unit circle) must lie in the stable region. Stability circles are a visualization tool for analyzing the stability of a system using the complex unit circle. In microwave theory, the goal is to find a suitable set of up/down-stream termination. The stability circles represent the boundaries of permitted terminations in their respective reflection coefficient planes. Since unconditional stability means that there are no possibilities that a passive termination trigger instabilities, a number of researchers worked to obtain a quality factor to check whether an amplifier is unconditionally stable.^27–31,21 It will be discussed in the following subsections.

Edwards-Sinsky ( μ ) criterion

Taking into consideration the prerequisites for unconditional stability, Edwards and Sinsky ²⁸ considered a highly unfavorable scenario where | R u p | = | R d n | = 1 . By deriving the expression for stability circles, they were able to establish the conditions necessary for an unstable region to exist outside the downstream unit circle.

{ C i − R i > 1 , if | S 22 | 2 − | Δ | 2 > 0 R i − C i > 1 , if | S 22 | 2 − | Δ | 2 < 0

(28)

Using algebraic work to quantify the minimum distance from the center of the unit circle to the unstable region, in both cases one arrives the similar relation as follows:

| S 22 − Δ S 11 * | + | S 12 S 21 | | S 22 | 2 − | Δ | 2 > 1

(29)

The denominator of equation (29) can be eliminated and the expression can be further simplified by noting that Edwards and Sinsky: ²⁸

| S 22 | 2 − | Δ | 2 = | S 22 − Δ S 11 * | 2 + | S 12 S 21 | 2 1 − | S 11 | 2

(30)

Finally, by introducing equation (30) in equation (29) a criterion known as the μ factor can be drived. This factor serves as a measure of the required margin of stability:

μ = 1 − | S 11 | 2 | S 22 − Δ S 11 * | + | S 12 S 21 |

(31)

In Figure 4, we have depicted several scenarios for better visual understanding. It should be noted that μ usually has different values in a frequency range and this should be independently interpreted. The red circles schematically illustrate the stability circle, while the white circle with a radius of unity represents the complex unit circle encompassing all possible values of reflection coefficients for a passive termination. Referring to Figure 4, if the μ factor is greater than or equal to 1 across the entire frequency range, a generic target device (in this case, a burner with a flame) is deemed unconditionally stable. This implies that the system will remain stable irrespective of the upstream and downstream terminations. When the μ factor is less than unity but still positive (i.e. 0 < μ < 1 ), even if it’s limited to a narrow frequency range, the center of the unit circle denotes the region of termination that ensures stability. In this case, the system is considered conditionally stable, as it necessitates specific termination conditions within that frequency range to maintain stability. In the scenario where − 1 < μ < 0 , the center of the complex unit circle represents the unsafe region corresponding to the termination at that specific frequency. Lastly, when μ = − 1 or μ < − 1 , the unsafe region encompasses the entire area of the complex unit circle even at single frequency, indicating that it is not possible to consider this system conditionally stable with infinite phase margin under.

OPEN IN VIEWER

Figure 4.

Schematic diagram indication of unsafe regions (red circle) and their relation with the μ factor in the complex unit circle (white circle).

Edwards and Sinsky proved that in the case of unconditional stability, that is, μ > 1 , | S 22 | should be smaller than unity. Therefore, for unconditional stability there is no need to evaluate directly | S 22 | < 1 according to equation (19).

Rollett ( k ) factor

To obtain the Rollett factor one can follow a similar procedure to form the μ factor, but this time starting from

R i − | C i | > 1

(32)

Substituting equation (27a) and (27b) into equation (32) results in

| S 22 2 − Δ S 11 * | − | S 21 S 12 | | S 22 | 2 − | Δ | 2

(33)

Rearranging equation (33) and squaring both sides

2 | S 12 S 21 | | S 22 2 − Δ S 11 * | < | S 22 2 − Δ S 11 * | 2 + | S 12 S 21 | 2 − ‖ S 22 | 2 − | Δ | 2 | 2

(34)

Note that the term | S 22 2 − S 11 * Δ | 2 can be expressed as follows:

| S 22 2 − S 11 * Δ | 2 = | S 12 S 21 | 2 + ( 1 − | S 11 | 2 ) ( | S 22 | 2 − | Δ | 2 )

(35)

Substituting equation (35) into equation (34) yields

2 | S 12 S 21 | < 1 − | S 11 | 2 − | S 22 | 2 + | Δ | 2

(36)

This leads to the definition of the Rollett stability factor as follows:

k = 1 − | S 11 | 2 − | S 22 | 2 + | Δ | 2 2 | S 12 S 21 |

(37)

and the Rollett criterion to ensure unconditional stability as follows:

k > 1.

(38)

Discussion II. The μ and k criteria are derived based on the second scenario and the conditions described in equation (18). Hence, these formulations cannot be directly applied in the presence of an unstable ITA mode. It is important to note that the μ criterion solely represents the minimum distance between the unsafe zones of the upstream or downstream unit circle of reflection coefficients and the center of the unit circle. When μ > 1 , it indicates unconditional stability that the entire unit circle is within the safe region. In the next section, we demonstrate that it is generally not feasible to meet this requirement across the entire frequency range in thermoacoustic applications. When − 1 < μ < 1 , this criterion does not provide any information about the size of the safe area inside the unit circle. Additionally, both μ and κ are derived under the assumption of fully reflective acoustic boundaries ( | R u p | = | R d n | = 1 ). In control theory, this assumption implies that the thermoacoustic system operates with infinite phase margin and under the most critical conditions for the gain margin.

Insight into the μ and k criteria

To obtain the scattering matrices of the burners with flames, one can substitute the flame TFs provided in Table 1 into equation (5a) to (5c) and then equation (3). By substituting the entries of the scattering matrix into equations (31) and (37), respectively, the variables μ and k can be easily calculated. Figures 6 and 7 depict the variations of k and μ with frequency, respectively. As expected, these variables exhibit similar qualitative and quantitative behavior. This alignment was anticipated since we demonstrated in the previous section that k can be derived from μ through algebraic manipulation. Another important observation is that both variables vary between − 1 and 1 , which may be attributed to the symmetry of the transfer matrix in the vicinity of zero Mach number. Furthermore, this figure clearly demonstrates that neither μ nor k can satisfy the requirement for unconditional stability, which necessitates values greater than unity for both variables in the entire frequency range.

OPEN IN VIEWER

Figure 6.

The Rollett factor k versus the frequency. The information on different flame transfer functions (TFs) are given in Table 1.

OPEN IN VIEWER

Figure 7.

The Edwards-Sinsky factor μ versus the frequency. The information on different flame transfer functions (TFs) are given in Table 1.

As mentioned earlier, μ represents the minimum distance from the origin to the unsafe region inside the unit circle. This feature of μ provides an advantage compared to κ . Although μ does not directly quantify the size of the unsafe region, a previous study by Kojourimanesh et al. ¹⁵ demonstrated that the size of the unsafe area inside the unit circle varies monotonically with μ . This relationship can be utilized in the analysis. Therefore, even though μ does not directly provide information about the size of the unsafe region, it can still effectively assess stability and enable relevant conclusions in the analysis, considering its monotonic relationship with the size of the unsafe area inside the unit circle. This distinction is important to make. Instead of using the term “unstable,” it is more accurate to refer to the regions as “unsafe” because the criteria are derived based on the concept of infinite phase margin and unconditional stability of subsystem.

An insightful visual representation of the situation in hand can be obtained if to notice that μ can be treated as a real-valued function of complex valued variable, namely of the flame TF. Accordingly, it can be directly visualized using, for example, a contour plots. Figure 8 provides a comprehensive representation of the μ variable by evaluating it over a wide range of generic flame TFs F = G e i φ . The contour plot illustrates the variations of μ , while the dotted circles within the contour represent iso-flame-TF-gain, and the dotted straight lines represent iso-flame-TF-phase. As an example, the specific flame TF F 1 is shown by a solid black curve within this plot. To highlight better the more critical region, we have plotted − μ , showing negative μ values in red, on this graph. It can be observed that when the phase of the flame TF crosses − π , it corresponds to the most critical region based on this criterion. This observation aligns well with the concept of an ITA mode. However, it appears to be more conservative than the requirements for the ITA mode, and can be demonstrated through algebraic analysis.

OPEN IN VIEWER

Figure 8.

The contour plot of the − μ factor in by considering arbitrary values for flame transfer function (TF).

From the concept of ITA in the vicinity of zero Mach number, ³² where T 22 = 0 , it is known that the pure ITA mode is unstable when the gain of the flame TF is larger than the critical value ( G ≥ G c = ϵ + 1 θ ) at frequencies where the phase crosses − π ± 2 m π , with m being an arbitrary positive integer. Based on equations (19) and (20), we can summarize the necessary requirement for μ > 0 and the necessary but not sufficient requirement for μ > 1 as | S 11 | = | S 22 | < 1 . By substituting F = G e i φ into T 21 = 0.5 ( ϵ − 1 − θ F ) and T 22 = 0.5 ( ϵ + 1 + θ F ) , and considering the phase of the flame TF F crossing − π ± 2 m π with m being any positive integer, we can express the conditions as follows:

| S 11 | = | S 22 | = | ϵ − 1 θ + G ′ | | ϵ + 1 θ − G ′ | < 1 ⇒ | S 22 | = n ′ + G ′ G c − G ′ < 1

(41)

G ′ < 0.5 ( G c − n ′ ) = 1 θ

(42)

For further clarification, a contour plot of the argument of T 22 and a contour plot of the module of S 22 are presented in Figures 9 and 10, respectively. In Figure 9, the concept of ITA mode can be interpreted. The π -wrapped argument of T 22 is shown, and a phase jump transition is observed when the phase of the TF crosses − π and the gain of the flame TF exceeds a critical value ( G c r ) slightly < 0.7 . Through calculations, it has been determined that this critical gain represents the threshold for the stable ITA mode. According to the requirements for the unstable ITA mode, a phase crossing of − π should occur along with a gain greater than G c r . This graph clearly shows that the blue-to-red transition (phase jump) only happens for gains greater than G c r . The white region within the contour plot corresponds to the argument of T 22 in the vicinity of 0 , and it occurs when the phase of the flame TF crosses 2 m π , where m can be any positive integer, with no restrictions on the gain.

OPEN IN VIEWER

Figure 9.

Contour plot of the argument of transfer matrix coefficient T 22 . The critical gain G c r is shown with a black diamond.

OPEN IN VIEWER

Figure 10.

Evaluation of the scattering matrix coefficient | R d n ‖ S 22 | with | R d n | = 1 .

This visual representation highlights that for the ITA mode, if the flame TF crosses the negative real axis in the polar plot with a sufficiently large gain ( G > G c r , which is necessary and sufficient for the occurrence of encirclement according to winding number), the ITA mode is present. Hence, we have visually demonstrated the ITA mode through this representation.

In Figure 10, the region where | S 22 | > 1 is highlighted in red, as evident, the critical region for the ITA mode ( G > G c r ≈ 0.7 ) with a phase crossing of − π is entirely situated inside the red zone, as expected from equation (42). In addition, the sufficient condition to ensure μ > 0 is G < G ′ when − 0.5 π < φ < − 1.5 π . Returning to the details in the derivation of the μ criterion, it is evident that this formulation has been derived under the assumption of fully reflective terminations (i.e. | R u p | = | R d n | = 1 ), which represents the most critical condition in the absence of the unstable ITA mode. To demonstrate this, we visualize S 22 > 1 / | R d n | in Figures 11 and 12 when R d n ≠ 1 . It can be observed that decreasing | R d n | leads to the critical (red) region, where the requirement for an infinite phase margin of S 11 and S 22 is violated, converging to the point G = G c r . Since we only consider flame TFs in the absence of the ITA mode, the gain at the frequency corresponding to the phase crossing of − π is less than G c r , indicating that the case of fully reflective terminations (i.e. | R u p | = | R d n | = 1 ) is the most critical scenario for the risk of thermoacoustic instability.

OPEN IN VIEWER

Figure 11.

Evaluation of the scattering matrix coefficient | R d n ‖ S 22 | with | R d n | = 0.5 .

OPEN IN VIEWER

Figure 12.

Evaluation of the scattering matrix coefficient | R d n ‖ S 22 | with | R d n | = 0.2 .

Discussion on the limitation of μ and k to be used as figure-of-merit

Having discussed the features and limitations of the μ and k factors, we now explore their potential as quality factors in thermoacoustics (in the absence of pure unstable ITA mode) for comparing and ranking different burners. To address this question, a reference dataset is required for comparison. One possibility is to employ a MC simulation on solving the dispersion relation (equation (7)). In this simulation, a library of 5000 sets of frequency-independent passive upstream and downstream reflection coefficients is created uniformly and randomly distributed inside the unit complex circle (see Figure 13). Through an iterative MC simulation, the stability of the system equipped with the considered burners with flames in combinations with sets of upstream and downstream reflection coefficients can be calculated in numerous cases, enabling the determination of the probability of stability as P S = N u / N t , where N u and N t are the number of unstable cases and the total number of cases, respectively. For more detailed description, see Kornilov and de Goey ¹⁰ and Saxena et al. ¹² . Figure 14 presents the histograms depicting the probability density of instability ( P D I ) versus unstable frequencies for different burners with flames. The distributions are compared with the − μ and − k criteria, revealing a strong correspondence between the negative values of μ and k and the concentration of unstable frequencies within a specific frequency range. It is worth noting that, as defined, the center of the unit circle represents the unsafe region, and the results provide indirect evidence that a more negative value of μ or k corresponds to an enlarged unsafe area within the unit circle. This finding aligns with previous investigations conducted by Kornilov et al. ¹⁴ and Kojourimanesh et al. ¹⁵ . Furthermore, Figure 14(d) offers a clear visualization of how the − μ and − k criteria effectively capture the distribution of unstable frequencies, even in scenarios where multiple unstable modes are present. Figure 15 shows the probability of instability concerning the number of unstable modes, which is defined as the ratio of the number of unstable cases to the total number of examined cases in the MC simulation. The results illustrate that B 1 ( F 1 ) , B 2 ( F 2 ) , and B 3 ( F 3 ) have only one unstable eigenmode each. In contrast, B 4 ( F 4 ) , which has been modeled with the n − τ model, exhibits multiple modes with equal probabilities. We here only evaluate first four unstable modes.

OPEN IN VIEWER

Figure 13.

Randomly generated upstream and downstream reflection coefficients as an input library for the Monte Carlo (MC) simulation.

OPEN IN VIEWER

Figure 14.

The histogram of distribution of unstable frequencies obtained from MC simulation, and a comparison with − μ and − k criteria. The information on different flame TFs are given in Table 1. MC: Monte Carlo; TFs: transfer function.

OPEN IN VIEWER

Figure 15.

Probability of instability of different modes of the system for the four flame transfer functions (TFs) according to Table 1.

Figure 16 displays the overall probability of stability for different burners based on the MC simulations, allowing for easy comparison and ranking. Probability of stability can be obtained easily by calculating P S = 1 − P I . Upon closer examination, it is evident that the B 1 ( F 1 ) , B 2 ( F 2 ) , and B 3 ( F 3 ) show almost similar thermoacoustic quality. However, B 1 ( F 1 ) slightly outperforms B 2 ( F 2 ) and B 3 ( F 3 ) , and B 3 ( F 3 ) outperforms B 4 ( F 4 ) . In this context, “better” implies a higher likelihood of thermoacoustic stability when combined with arbitrary sets of upstream/downstream embeddings.

OPEN IN VIEWER

Figure 16.

Total probability of stability of the system for burners with flames. The information on different flame transfer functions (TFs) are given in Table 1.

It is important to note that the performed MC simulation generated a library of constant (frequency-independent) uniformly random reflection coefficients, which represents a specific class of reflection coefficients. Therefore, it is not appropriate to generalize the obtained probability as a definitive measure. The MC simulation based on the constant reflection coefficients may indicate a good burner with a low probability of stability, but it is not possible to identify inherently bad burners with a high probability of stability. Keeping in mind, this limitation of the MC simulation, we can compare the burner ranking obtained from the MC simulation with that derived from the μ (or k ) criterion.

Figures 6 and 7 demonstrated that the variables μ and k exhibit a strong dependence on the frequency and can qualitatively indicate the concentration of unstable frequencies. To compare the burners quantitatively, a single performance metric or quality factor needs to be calculated. The μ criterion represents the minimum distance to the unsafe region, where a higher value signifies a greater distance between the unsafe region and the origin, indicating a potentially safer operating condition. In order to obtain a single numerical value for each burner with a flame, the average μ value was computed over the entire frequency range, as depicted in Figure 17.

OPEN IN VIEWER

Figure 17.

The averaged value of μ to be used as figure-of-merit to rank different burners. The information on different flame transfer functions (TFs) are given in Table 1.

This averaging process facilitates a straightforward comparison among the burners, with higher average μ values may indicate a higher likelihood of thermoacoustic stability, and lower average μ values suggesting an increased risk of instability. Based on the rankings derived from the average μ values, the burners can be ranked in terms of their stability performance. B 1 ( F 1 ) , possessing the highest average μ value, is considered superior to B 3 ( F 3 ) . Similarly, B 3 ( F 3 ) surpasses B 2 ( F 2 ) due to its higher average μ value, indicating a more favorable stability profile. Finally, B 2 ( F 2 ) exhibits a higher average μ value compared to B 4 ( F 4 ) , thereby suggesting a better likelihood of stability. Upon comparing the results of the μ criterion with the MC simulation, it is worth noting that the ranking of burners based on the μ criterion does not align with the ranking based on the MC simulation. Specifically, burner 2, which is identified as the second-best burner according to the MC simulation with slight thermoacoustic difference with burner 1 and burner 3, does not exhibit the second-highest μ value. This discrepancy can be attributed to the fact that the μ criterion heavily relies on the concept of unconditional stability of subsystem (e.g. input reflection coefficient).

Upon closer examination of the flame TF of B 2 ( F 2 ) in Figure 5, it is apparent that its phase gradually decreases over a wide frequency range, remaining around − π . Consequently, the μ value for B 2 ( F 2 ) in Figure 7 is significantly negative. This behavior can lead to a misleading ranking when using the μ criterion, particularly for burners whose flame TF saturates around − π , almost irrespective of their gain values. A similar conclusion can be drawn for the k criterion as well.

Moreover, when comparing B 2 ( F 2 ) and B 3 ( F 3 ) based on the MC simulation, it is observed that they exhibit quite similar probabilities of stability. However, the μ criterion shows a significant difference between them, potentially leading to a misleading ranking. Additionally, the μ criterion suggests that B 4 ( F 4 ) has a high risk of instability, while the MC simulation demonstrates that it may actually operate acoustically stable in ∼ 60% of the situations. This discrepancy highlights the limitations of the μ criterion in accurately predicting the thermoacoustic stability of burners. In summary, the concept of infinite phase margin and unconditional stability of subsystem, which underlie the μ and k criteria, yield very conservative quality factors that can mislead to evaluate properly thermoacoustic quality of burners with their associated flames, even in the absence of pure unstable ITA mode. Nevertheless, they can still qualitatively identify the frequency range of potential instability correctly.