Open access peer-reviewed chapter

Combustion Processes with External Harmonic Excitation using Extended Lindstedt-Poincare Method with Multiple Time Scales

Written By

Rudolf R. Pušenjak and Igor Tičar

Submitted: 21 October 2015 Reviewed: 21 March 2016 Published: 05 October 2016

DOI: 10.5772/63220

From the Edited Volume

Developments in Combustion Technology

Edited by Konstantinos G. Kyprianidis and Jan Skvaril

Chapter metrics overview

1,336 Chapter Downloads

View Full Metrics


This work treats the control of nonstationary pressure oscillations of combustion processes subject to external harmonic excitation. The low-order van der Pol model of nonstationary combustion processes with two degrees of freedom and a multiplicative feedback control function is considered. Nonstationary oscillations of the combustion process with asynchronous and synchronous excitation, respectively, are analyzed by means of the extended Lindstedt-Poincare (EL-P) method with multiple time scales. The steady-state amplitudes of competitive quenching with asynchronous excitation are not so much affected in comparison with steady-state amplitudes of an autonomous system, opposite to the synchronous excited system. The influence of control parameters on the quenching of the externally excited oscillations is investigated in details.


  • combustion processes
  • pressure oscillations
  • external harmonic excitation
  • van der Pol model
  • EL-P method with multiple time scales
  • control

1. Introduction

The phenomenon of the instability of the combustion process is commonly explained by the nonstationary flame causing an acoustic wave, which reflects from the walls of the combustion chamber back to the combustion process. The thermoacoustic coupling between acoustic waves and unsteady heat release using a generalized energy equation, followed by feedback control concepts is presented by Dowling and Morgans [1]. Using theory of systems, this coupling can be explained alternatively by means of a feedback with a negative damping between the heat release process and acoustic wave propagation in the combustion chamber. The recent researches in the area of combustion process attempt to eliminate the unstable combustion and reduce the emission pollution, especially the emission of nitrogen-oxygen compounds at the same time. The reduction of pollution can be achieved by decreasing the fuel-to-air ratio; however, this causes the instability in the form of self-excited oscillations, which in general contain incommensurate frequencies. The problem of instability can be solved by the control of the fuel influx into combustion chamber, which can be realized by introduction of multiplicative control function performing the desired modulation of heat release rate. In order to obtain an efficient model of active control in this chapter, the reduced-order model of combustion instability process is proposed in the form of generalized, coupled van der Pol differential equations [2]. This model is built using a grey-box system approach [3], which is based on combination of physical knowledge determining the structure and identification of unknown parameter values, which are fitted to the measured input/output data of the system. In the proposed low-order model, the physics is present only in a feedback coupling between acoustic wave propagation and heat release giving rise to a negative damping as the cause of the combustion instability. The presented model of combustion instability is based on the linear acoustics and propagation of standing acoustic waves in the combustion chamber. Due to the presented nonlinear coupling between linear acoustics and nonlinear heat release rate, theoretically an infinite number of acoustic modes exist. Based on experimental results of Sterling [4] and Dunstan [3], the following analysis is considering the temporal variations of first two acoustic modes with corresponding natural frequencies. Unlike the models of Peracchio and Proscia [5], Dunstan and Murray et al. [6], respectively, natural frequencies in the present model can be incommensurate in general. From the system point of view, acoustic waves are modeled as pure linear oscillators and the heat release rate is modeled as the static nonlinearity, which can be accurately adapted to the real situation by identification of system parameters. On this basis, the model can be stated as a set of partial differential equations (PDEs) in terms of time scales, which are ideally suited to be attacked by the extended Lindstedt-Poincare (EL-P) method. The spatial scales are not present in the model, because it is assumed that all heat is released into the fluid at a point [5].

This chapter is organized as follows. In Section 2, the control structure and the governing coupled van der Pol equations of controlled combustion process with external harmonic excitation are presented. Section 3 applies the EL-P method with multiple time scales [7] to control problems of combustion process with asynchronous external harmonic excitation. Section 4 deals with control problems of combustion process with synchronous external harmonic force. Section 5 focuses on results and discussions of competitive quenching phenomenon. Section 6 discusses the validation of the proposed van der Pol model of combustion processes, and Section 7 outlines practical aspects of the model implementation. Section 8 contains conclusions.


2. Van der Pol model of combustion process with external harmonic excitation

In this section, the phenomenon of instability of the combustion process due to the feedback with negative damping is described by means of the low-order model, which consists of generalized two degrees of freedom (2-DOF) van der Pol oscillator with external harmonic excitation. The coupled van der Pol oscillator model is experimentally verified and it is possible to analyze the appearance of excited oscillations, which can evolve to limit cycles in accordance with experimental observations [8]. The dynamics of the generalized 2-DOF coupled nonautonomous van der Pol oscillator is governed by the following equation:


where the left-hand side of Eq. (1) models the linear acoustics of the combustion chamber represented by two harmonic oscillators with ωn1 and ωn2 being the natural frequencies of first two temporal pressure oscillation modes y1 and y2, respectively. The first term on the right-hand side of Eq. (1) models the nonlinear heat release process, where the sum y1+y2 denotes the temporal variations of the downstream pressure p at the burning plane and function q=φ0+φ1(y1+y2)13φ3(y1+y2)3 represents the flame heat release rate in dependence on pressure p=y1+y2 . Parameters φ0, φ1 and φ3 are determined experimentally through system identification. The parameter φ0 stands for a drift, which is different from zero, so that φ0 is an arbitrary (positive or negative) constant. The parameter φ1 is a slope of function q about origin, which corresponds to a positive constant. The term 13φ3(y1+y2)3of the heat release rate static characteristic describes saturation phenomenon, where parameter φ3determines the level of saturation. This level is also defined by a positive constant. The parameter ε is the gain of the heat release process. To justify the use of EL-P method in the analysis of combustion process, the feedback nonlinearity is assumed to be small. This requirement is achieved by multiplication of the time derivative of heat release rate by a small positive number ε>0. Note, that derivative term dqdt=d[φ0+φ1(y1+y2)13φ3(y1+y2)3]/dt on right-hand side of Eq. (1) represents a negative damping, which in the absence of external excitation can drive self-excited oscillations. Although temporal acoustic oscillations remain linear, nonlinearity occurs in the heat release rate, which affects the initial growth of excited oscillations until a limit cycle is reached due to the heat release saturation. The second term F cos ωt on the right-hand side of Eq. (1) models the external harmonic excitation, where F is the excitation amplitude and ω is the excitation angular frequency. This term can be seen as a disturbance, which allows analyzing the system subjected to the noise. The proximity of excitation frequency ω to natural frequencies ωn1 and ωn2decides whether an interconnection between these frequencies exists. If the excitation frequency is far away of natural frequencies, all three frequencies are not interconnected and the system is driven by asynchronous force. On the other hand, when ω is close to one of two natural frequencies, there exists an interconnection between both frequencies and the system is driven by synchronous force. The thermoacoustic model of the combustion process is depicted in Figure 2a, where harmonic oscillators modeling acoustics are placed in two parallel forward branches and the heat release process is modeled in a feedback loop. On the basis of assumption about the spatial coherence, the proposed model captures the temporal variations of pressures and heat release rate, while the spatial distributions are omitted.

Figure 1.

Van der Pol model of the combustion process with external harmonic excitation. (a) Model without control, (b) Model with control.

The coupled van der Pol oscillator in Eq. (1) can be represented by the control structure, shown on Figure 1a, where the control plant consists of two harmonic oscillators in forward branches, which are harmonically excited and the feedback consists of derivative of the flame heat release rate multiplied by ε. By expanding the derivative of the flame heat release rate, a negative damping term appears, which is the cause of instability of the combustion process. The instability is reflected in self-excited/driven oscillations, which evolve to the limit cycles. The limit cycle amplitudes are determined by the nonlinearity of the heat release rate function. The growing self-excited/driven oscillations can be so intense as to lead to a structural damage. To prevent such consequences, the control of oscillations can be introduced by means of multiplicative control function Ψ(y1+y2), which is based on measurement of pressure p=y1+y2and placed in an additional feedback branch as shown in Figure 1b. Multiplicative control function Ψ(y1+y2) is chosen in order to achieve the modulation of a fraction of the fuel influx into the combustion chamber and consequently to affect the heat release rate. Thus, by applying the control function Ψ, the right-hand side of the governing Eq. (1) is modified in the form:


The suitable candidates of control function Ψ(y1+y2) are sought in the form of polynomial functions of pressure p=y1+y2. As suggested in reference [2], the control function Ψ(y1+y2) is proposed in the form:


where K1 denotes the amplifying factor, which is of the same sign as the constant φ0 and K2 is a positive constant. Besides, the constant K2 must be greater than unity, K2>1, as is shown later. When both K1 = K2 = 0, the system is without the control. The low pass filter, shown in Figure 1b, is optional. Its role is discussed in details in the paper by Pušenjak et al. [2] and will not be repeated here. Due to the low-order van der Pol model of combustion process, the implementation of the control algorithm is very fast and the active control is very effective.

The stability analysis in the scope of the active control can be applied to reduce NOx pollution by extending the stability domain of computed pressure oscillations [9]. By applying the presented model, a realistic modeling of the flame heat release rate is quite possible, when parameters φ0, φ1 and φ3, respectively, are properly adjusted by using identification methods. By such an approach, the present method can be applied to model combustion instabilities in gas turbines with lean premixed combustors and power plants.

In the next section, the EL-P method with multiple time scales [7] is applied for computation of pressure oscillations with external harmonic excitation arising in the coupled nonautonomous van der Pol oscillator. For the sake of compactness, we introduce new parameters:


and rewrite governing equation (2) with control function (3) in the form:


3. Application of the extended Lindstedt-Poincare method with multiple time scales in the combustion process with asynchronous external harmonic excitation

The instability problem of the combustion process, which is posed in the form of coupled van der Pol Eq. (5) with external harmonic excitation, can be analyzed by EL-P method with multiple time scales [2]. To provide a perturbation analysis of the generalized model of instability, it is assumed at first that interaction between both resonators and external harmonic excitation does not exist. In such a case, the natural frequencies ωn1 and ωn2 are arbitrary and do not meet any of the special conditions ωn1=ωn2, ωn1=3ωn2 or ωn1=⅓ωn2. In addition, it is assumed, that the excitation frequency ω is far away of both natural frequencies ωn1 and ωn2. When these conditions are met, we say that such a combustion process is driven by asynchronous external harmonic force. As we will see, the pressure oscillations developed in this case are characterized by the phenomenon of the competitive quenching. In this phenomenon, one of the resonators is excited and oscillations of the other resonator are quenched in dependence on initial conditions. On the other side, the excited oscillations under special conditions ωn1=ωn2, =3ωn2 or ωn1=ωn2correspond to the several kinds of interaction between both resonators, where mutual synchronization phenomenon with close frequencies or multiple harmonics, respectively, appears. Owing to the limited space, perturbation analysis of such oscillations will not be treated.

Following the EL-P method with multiple time scales, two fast time scales τ1=ω1t and τ2=ω2t are introduced, where ω1 and ω2 are termed the nonlinear frequencies of coupled van der Pol oscillators in accordance with Eq. (12). In addition, the time scale τ3=ωt, which corresponds to the excitation frequency ω and a slow time scale τ4=εt with a small positive expansion parameter ε are introduced. Using time scales, the time derivatives d/dt and d2/dt2 are replaced by differential operators:


Owing linear acoustics, following substitutions are proposed:


and by using differential operators (6a, 6b), governing equations of the coupled van der Pol oscillators [Eq. (5)] are rewritten into partial differential equations:


After performing indicated partial derivations, one obtains:


When the system is not subject of control, K1 = K2 = 0, then γ0 = φ0, γ1 = φ1, γ2 = 0, γ3 = −⅓φ3, γ4 = 0, γ5 = 0, γ6 = 0 and Eq. (9) is simplified in the form:


The approximate solution of Eq. (9) in dependence on three time scales τ1, τ2 and τ4 is sought in the form of power series:


Similarly, the nonlinear frequencies ωi, (i=1,2) are expressed in the form of power series about linear frequencies ω10 and ω20, respectively:


where frequencies ωik, (i=1,2) of power series (12) may be slowly varying parameters in dependence on the slow time scale τ4. Frequencies ωik, (i=1,2) are unknown in advance and must be determined successively during the perturbation procedure. As we will see later at perturbation step of zeroth order in Eq. (9), frequencies ωi0, (i=1,2) are linearly dependent on natural frequencies ωni and thus are termed linear frequencies. Frequencies ωik,(k ≥ 1) are determined at perturbation step of kth order from so-called solvability conditions and may be in general nonlinear functions of slowly varied amplitudes. Therefore, frequencies ωi, (i=1,2) are termed nonlinear on the basis of Eq. (12). By substituting the assumed solution [Eq. (11)] and the power series [Eq. (12)] into Eq. (9) and then equating terms of equal power of ε on both sides of equation, one obtains a set of PDEs, which must be successively solved as perturbation steps of ascending order:


for i=1,2.

The right-hand side of Eq. (13) is only a function of time scale τ3, while the left-hand side of Eq. (13) contains dependencies on time scales τ1, τ2 and τ4. Thus, both sides of Eq. (13) are independent on each other and the right-hand side must satisfy the condition F + (ω2ωni2)Ei = 0. According to the supposition, the excitation frequency ω is far away from each of natural frequencies ωni, (i=1,2) and amplitudes Ei, (i=1,2) of solution (7a) can be determined from this condition as:


As mentioned earlier, in such circumstances we have deal with the asynchronous case of external excitation. When the excitation frequency is close or equal to one of the natural frequencies, the corresponding amplitude Ei cannot be uniformly determined by means of Eq. (16). In this case, we have a synchronous excitation and the corresponding solution (7a) requires setting of the amplitude Ei equal to zero. The perturbation procedure for the synchronous excitation force will be treated in details in the next section. In the sequel, the general asynchronous case is analyzed.

The general solution of zeroth order, which allows two independent natural frequencies ωn1 and ωn2, is sought in the form:


from where it follows that the self-excited/driven oscillations of coupled autonomous van der Pol oscillator, which correspond to the zero-order solution, are foreseen as almost periodic solutions with slowly varying amplitude and phase. Inserting assumed solutions into zero-order PDE gives an important relation between linear and natural frequencies of the oscillator:


which reveals, that linear frequencies are equal to the natural frequencies of the coupled van der Pol oscillator under assumption that linear frequencies can be only positive. Consequently, it follows from Eq. (12) that nonlinear frequencies ω1 and ω2 are independent in this case, too. By inserting the solution (17) into first-order PDE (14), one obtains:


Secular terms on the right-hand side appear, which causes an unbounded increase of the solution. To obtain an uniform solution, secular terms must be removed. The elimination of secular terms is achieved by assembling terms, which appear at functions cos[τi+Φi(τ4)] and sin[τi+Φi(τ4)], (i=1,2), respectively, and equating obtained expressions by zero. This task is accomplished by symbolic computation in Mathematica® [2] and results in solvability conditions of first order:


for i=1,2. If the combustion process is not externally excited, the system of governing Eq. (9) is autonomous, where p = 0, Ei = 0, (i=1,2) and E = 0. By inserting expressions (4) for parameters γj, (j=1,…,6), solvability conditions (21) can be written in the form:


Solvability conditions, given by means of Eqs. (20) and (22) for an autonomous case, agree perfectly with solvability conditions of first order as were derived in reference [2].

From Eq. (20) it follows that trivial solution A1(τ4) = A2(τ4) = 0 always exists. In the case of nontrivial solution Ai(τ4) ≠ 0, (i=1,2) one has:


and frequencies ωi1, which appear in Eq. (23) are unknown at this stage of approximation and their computation is deferred to the second-order approximation. The nontrivial solution Ai(τ4) ≠ 0, (i=1,2) can be obtained by numerical integration of Eq. (21), although the closed form of solution is obtainable in exceptional cases [2].

The solution of Eq. (21) tends to reach the asymptotic values, when τ4→∞. In this case, limτ4Ai(τ4)=0 and steady-state solution Ais=limτ4Ai(τ4) can be obtained by solving equation:


where following possibilities exist: (i) a trivial solution Ais = 0, (i=1,2), (ii) a nontrivial solution Ais0(i=1,2),A1s=±A2s,Ais|1,2=±23γ1γ332E2 and (iii) a nontrivial solution Ais0,(i=1or2), coupled with trivial solution A(3-i)s = 0, where i=1or2. The solution in this case isAis|3=±4γ13γ32E2,A(3i)s|4=0. By substituting expressions for constants γ1 and γ3, one obtains asymptotic values in explicit form:


From Eq. (25) it follows that by applying the control, the parameter K2 may not take values in the range 0< K2≤1 in order to prevent complex roots or the divergent behavior of asymptotic values for K2 = 1, respectively. Therefore, the parameter K2 must be greater than unity, K2>1. If control is not applied, both constants K1 and K2 are equal to zero, K1 = K2 = 0, giving the steady-state (asymptotic) values:


From Eqs. (25) and (26), respectively, it follows that all nontrivial steady-state values of the combustion process with external harmonic excitation are affected to some extent, when compared by asymptotic values of self-excited oscillations [asymptotic values of self-excited oscillations can be obtained by inserting value E = 0 into Eq. (25) or Eq. (26), respectively].

To proceed with the perturbation procedure, generally the almost periodic solution of the first order is constructed, which contains all possible combination tones and has the form:


Coefficients Bim,(2m+1)k+j (τ4) and Cim,(2m-1)k+j (τ4), respectively, can be computed by putting the solution (27) into Eq. (19), removing secular terms from its right-hand side and equating like terms on both sides.

The second-order approximation is constructed in the same manner as the first-order approximation. To do this, one inserts solutions (17) and (27) into second-order PDE (15). The right-hand side of the obtained equation contains secular terms, which must be removed to prevent the unbounded increase of the solutions. The elimination of secular terms, which appear at functions cos[τi+Φi(τ4)], (i=1,2), is again achieved by assembling terms and equating obtained expressions by zero. Because this task is extensive, it is accomplished by symbolic computation in Mathematica®. The resulting solvability conditions offer a relationship between unknown nonlinear frequencies ωi2, (i=1,2) and coefficients Ai(τ4), Bim,(2m+1)k+j (τ4) and Cim,(2m-1)k+j (τ4), respectively. Due to the lengthy formulas, they will not be presented here. The elimination of secular terms, which appear at functions sin[τi+Φi(τ4)], (i=1,2), respectively, and equating obtained expressions by zero, leads to solvability conditions:


After determination of solvability conditions of the second order, the perturbation procedure may be continued with the computation of the second-order solution xi2(τ1,τ2,τ4). Under hypothesis of the small expansion parameter ε1, one observes, that this task is sometimes unnecessary on the basis of the following test. By computing frequencies ωi2, (i=1,2) from solvability conditions at a given instant of the slow time scale, one can check if the terms ε2 ωi2 of the power series (12) are sufficiently small under the prescribed tolerance. If they are so on the entire interval of the slow time scale, then the first-order perturbation analysis is adequate. Now, by considering the solvability condition of the first order (20) for the nontrivial case Ai(τ4)≠0, (i=1,2), that is the simplified condition ωi1+Фi'(τ4)=0 in Eq. (28), one obtains:


By considering the above result in solvability condition (20) again, it follows:


From Eq. (30), it follows that frequencies ωi1 are constant. Thus, Eq. (29) can be integrated as follows:


From Eq. (31), one concludes that strictly speaking, phases are linear, slowly varying functions of the time scale τ4. However, keeping in mind the definition of the slow time scale τ4 = εt, we can approximate Eq. (31) with Φi(τ4)≈Φi(0), which implies that derivatives of phase angles are equal to zero, Φi(τ4)=Φi(0)=0. By considering this in Eq. (30) again, we finally obtain:


The obtained result simply means, that nonlinear frequencies ωi of the first order analysis of van der Pol type nonlinear system can be approximated by the zero order term of Eq. (12). Thus by considering Eq. (18), nonlinear frequencies ωi can be approximated by the corresponding natural frequencies ωni:


Now, the following important observation can be made. Substituting Eq. (32) into Eq. (20) and assuming K1 = K2 = 0 in Eq. (22) to consider the uncontrolled case, the solvability conditions of the first order are simplified on the form:


for i=1,2. This result is exactly the same as derived by Krilov-Bogoliubov averaging method (K-B method) in Landau et al. [8] (with ε = 1, so that the slow time scale τ4 = εt in the present method corresponds to the time t in the K-B method). Consequently, there is no difference in computed zero order amplitudes and phases of pressure oscillations, regardless of whether presented EL-P method or Krilov-Bogoliubov averaging method, respectively is applied.


4. Application of the extended Lindstedt-Poincare method with multiple time scales in the combustion process with synchronous external harmonic excitation (ω close to ωn2)

In this section, the combustion process with synchronous external harmonic excitation is treated. Although we will assume in the sequel, that excitation frequency ω is close to the second natural frequency, ωωn2, the case with ωωn1 can be treated in an analogous way. To stress the nearness of the natural frequency ωn2 to the excitation frequency ω, we introduce the parameter Δ by means of equation:


where ε is the same small positive expansion parameter as appears in definition of the heat release process, Eq. (1). Opposite to the asynchronous case, the second mode must be driven by the soft harmonic excitation force, which is given in the form εf cos ωt, (F = εf). For convenience, the excitation frequency ω is renamed to ω2 in order to reduce the number of time scales. With this device, three time scales τ1 = ω1t, τ2 = ω2t = ωt and τ3 = εt are used in the perturbation analysis. Under these conditions, the synchronization phenomenon can be described by means of equations:


Equation (7.a) for the first mode, i=1, is now rewritten in the reduced form:


where solution of PDEs (37.a,b) is sought by means of power series:


Performing perturbation steps, which are presented in details in the asynchronous case and finally assembling secular terms, which appear at functions cos[τi+Φi(τ3)] and sin[τi+Φi(τ3)], (i=1,2), respectively, one obtains following solvability conditions of first order:


Solvability conditions (40a) and (40c) for the first mode pressure oscillation, i=1, are little changed, when compared by solvability conditions (20) and (21), respectively for asynchronous excitation. Formally, amplitude E = E1+E2 in Eq. (21) must be replaced with E = E1, (E2 = 0) in solvability conditions (40c, 40d). However, solvability conditions (40b) and (40d) for the second mode pressure oscillation, i=2, are much more changed in comparison with asynchronous case. There is a substantial difference between phase courses Φ13) and Φ23) provided by solvability conditions (40a) and (40b), respectively. While Eq. (40a) produces a linear dependence of phase angle Φ13) on the slow time scale τ3 [compare Eq. (31) in Section 3], Eq. (40.b) produces a variable phase angle Φ23) with complicated dependence on the slow time scale, which cannot be obtained analytically.


5. Results

In this section we investigate self-excited as well as externally driven oscillations of coupled van der Pol model of combustion processes. Besides this, the effect of the control in order to suppress the nonzero amplitude of the limit cycles is analyzed. The control function is activated, when the switching time ton is reached, which is defined as the time when the slowly varying amplitude approaches a steady-state value within the prescribed tolerance range. The parameter tol is chosen to define the allowed tolerance range and has the value tol = 10-6 in all examples, which are analyzed bellow. Inside the time interval 0 ≤ t < ton, the control is not active and both parameters K1 and K2 are zero valued (K1 = K2 = 0), while in the time interval ton ≤ t ≤ ∞ the control is applied with a positive valued parameter K2 and a gain factor K1 with the same sign as the parameter φ0. In what follows, results of the analysis of the competitive quenching phenomenon by using coupled van der Pol model of combustion process are presented. Although the presented general theory is valid for a small expansion parameter ε1, the reasonable results can be computed even for the value of the expansion parameter ε = 1. Nevertheless, Lindstedt-Poincare method with multiple time scales can be modified with strong nonlinearities and large expansion parameter ε [10]. The parameters of the combustion process analyzed in this section have values: φ0 = 0.5,φ1 = 1.78 × 102,φ3 = 1.24 × 107 and ε = 1. Thus, the slowly time scale τ4 = εt can be replaced by the time t for the sake of simplicity. By using this agreement, dependencies of all variables on the slow time variable τ4 are denoted in Figures 2 and 3 as dependencies on the time t.

The phenomenon of competitive quenching arises, when nonlinear frequencies ω1 and ω2 do not have interconnection between itself nor between first odd subharmonic or first odd superharmonic, respectively. In other words, the phenomenon of competitive quenching appears, when ω1 is not close to ω2 nor 3ω2 or ⅓ω2. The general theory presented in Section 3 is applicable to this phenomenon on the whole, where important conclusions about the phenomenon are analyzed in details.

Figure 2.

Comparison of pressure oscillation in the uncontrolled combustion process. (a) Pressures x1(t) and x2(t) without external excitation. (b) Pressures y1(t) and y2(t) with external harmonic excitation.— Oscillation computed by EL-P method. --- Envelope computed by EL-P method.

Characteristic for the competitive quenching phenomenon is excitation one of resonators, while oscillations of the other resonator are quenched in dependence on initial conditions. Such oscillations are analyzed in the case of natural frequencies ωn1 = 2π × 210 s-1, ωn2 = 2π × 740 s-1. The characteristic properties of the competitive quenching phenomenon at selected initial conditions A1(0) = 4 × 10-3, A2(0) = 2 × 10-3 are computed for the value of expansion parameter ε = 1 and are shown in Figures 2 and 3.

The combustion process is first analyzed without applied control by setting the control parameters K1 and K2 to be equal to zero, K1 = K2 = 0. Computed pressure oscillations of combustion process without control are plotted in Figure 2. In Figure 2a, first- and second-mode pressure oscillations x1(t) and x2(t), respectively are shown, which are computed for the combustion process without external excitation. When external harmonic excitation is applied with excitation amplitude f = 104 and with corresponding excitation frequency ω = 2ωn2, respectively, the corresponding pressure oscillations y1(t) and y2(t) look like on Figure 2b. Comparison of oscillations x1(t) and y1(t) shows, that external harmonic excitation has a very small impact on the first oscillation mode, so that differences are not visually observable. However, computation of steady-state amplitudes shows a small deviation. The steady-state amplitudes for the uncontrolled combustion process are computed by using Eq. (26), which determines the steady-state amplitude of harmonically excited response y1(t) by means of expressionA1s|3=±2φ1/φ3E2/2, while the steady-state amplitude of response x1(t) without the external excitation is obtained by the expression A1s|3=±2(φ1/φ3). When harmonic excitation is not applied, E = 0, the steady-state amplitude is A1s|3=±0.007578and with the applied harmonic excitation, the steady-state amplitude is A1s|3=±0.007568. We can see, that the difference between steady-state values is very small, but the computed values can be approximately checked on Figure 2.

On the contrary, comparison of second mode oscillations x2(t) and y2(t) reveals significant differences. Steady-state amplitudes of both responses are in accordance with Eq. (26) equal to zero, A2s|4=0. Zeroed steady-state amplitudes of second-mode oscillation can be seen in Figure 2 by dashed curve, which presents the envelope of zero-order solution (17).

Figure 3.

Comparison of pressure oscillation in the controlled combustion process. (a) Pressures x1(t) and x2(t) without external excitation. (b) Pressures y1(t) and y2(t) with external harmonic excitation. — Oscillation computed by EL-P method. --- Envelope computed by EL-P method.

Combining zero-order solutions xi0(τi, τ4) by using Eq. (17) and first-order solutions xi1(τ1,τ2,τ4) by Eq. (27), respectively, in power series (11), one obtains corresponding courses of self-excited oscillations xi(τ1,τ2,τ4) ≈ xi0(τi, τ4)+ε xi1(τ1,τ2,τ4), (i=1,2) for the combustion process without external excitation taking E = 0 into account. Self-excited time responses are plotted by continuous lines in Figure 2a.

On contrary, harmonically excited counterparts are computed by means of Eqs. (7a, 7b), where Eq. (11) is considered in approximation xi(τ1,τ2,τ4) ≈ xi0(τi, τ4)+ε xi1(τ1,τ2,τ4), (i=1,2). The corresponding time responses are plotted by continuous lines in Figure 2b. The most important difference occurs between second-mode time responses x2(t) and y2(t), respectively, after the passage of the transient phenomenon.

While transient self-excited oscillation x2(t) without external harmonic excitation expires in a sinusoidal oscillation with corresponding limit cycle, transient oscillation y2(t) with external harmonic excitation expires in a composite oscillation of aperiodic nature.

As is clearly seen from Figure 2, the amplitudes of oscillations xi and yi, respectively, slightly differ from the zero-order amplitudes A1(τ4), A2(τ4) due to the first-order contributions xi1(τ1,τ2,τ4), (i=1,2) added to the zero-order solutions xi0(τi, τ4). Thus, the first-order contributions are of the moderate size than the zero-order solutions, and it can be expected that the contributions of the second-order solutions xi2(τ1,τ2,τ4), (i=1,2) are still smaller and thus can be neglected. From this reason, it follows that the first-order EL-P analysis is sufficient even if the expansion parameter takes the value ε = 1.

The characteristic properties of the competitive quenching phenomenon depicted in Figure 3 indicate the effect of the applied fuel influx control on the combustion process. The control is activated in the switching time ton, which is computed by means of the numerical integration of solvability conditions (22) until amplitudes A1(τ4), A2(τ4) approach the steady-state values, A1s|3=±2φ1/φ3E2/2 A(3-i)s|4 = 0 within the tolerance range tol = 10-6. By activating the control function in the switching time ton = 0.1, which is determined in a way as is described, the values of control parameters are changed to K1 = 0.5, K2 = 2.

Figure 3 shows the courses of slowly varying amplitudes A1(τ4), A2(τ4), drawn by dashed lines, while the corresponding pressure oscillations xi(τ1,τ2,τ4) ≈ xi0(τi,τ4)+ε xi1(τ1,τ2,τ4), (i=1,2) without external excitation and responses yi with external harmonic excitation, respectively, are drawn by continuous lines.

From Figure 3a, it follows that pressure oscillations xi(t) without external excitation are successfully quenched to zero after the activation of the control in the switching time ton = 0.1.

On contrary, Figure 3b reveals that quenching of pressure oscillations yi(t) after the activation of the control is not perfect, because a small composite oscillation of aperiodic nature remains due to the constantly present external harmonic excitation. However, the remaining oscillation with very small amplitude is acceptable.


6. Validation of the van der Pol model of combustion processes

In the combustion chamber, theoretically an infinite number of acoustic modes of pressure wave propagation can exist. However, experimental researches of combustion processes, done by Sterling [4] and Dunstan [3], indicate the existence of two dominant longitudinal acoustic modes with corresponding natural frequencies. This fact validates the choice of 2-DOF van der Pol model of combustion processes, proposed in this chapter. Another experimentally verified feature of acoustic modes is occurrence of incommensurate natural frequencies [8]. Considering this feature in the proposed van der Pol model leads to the multiple time scales concept. Due to the strong nonlinear coupling, various combination tones emanate in addition to dominant natural frequencies. As mentioned in Section 3, first-order solution in Eq. (27) contains all possible combination tones with 21 corresponding Fourier coefficients Bm,(2m+1)k+ji(τ4) and 11 Fourier coefficients Cm,(2m1)k+ji(τ4), respectively, for both modes (i=1,2), which must be determined. Computation of Fourier coefficients using EL-P method is validated by comparison with K-B method [2]. By considering all combination tones of Eq. (27) in K-B method, an exact agreement of computed coefficients in both methods is found.

Figure 4.

Van der Pol model of the combustion process with transport time delay θ.

The choice of the van der Pol model of combustion processes is justified especially due to the existence of limit cycles. Limit cycles in combustion process arise as a result of heat release saturation. This feature of the van der Pol model is particularly valuable because it allows to predict the highest level of the acoustic pressures by means of the calculation of the steady-state values, see Eq. (25). The static characteristic of the heat release process is described by means of the generalized nonlinearity of van der Pol type, which includes the saturation term. The generalized static characteristic of heat release rate contains three parameters φ0,φ1and φ3,, respectively, which must be determined by identification methods so as to provide the best matching with a real heat release process.

A further validation can be obtained by modification of van der Pol model to consider a transport time delay θ from nozzle to flame surface of combustion processes. For this purpose, the time delay, differentiator block and low pass filter are included in the feedback path of a modified van der Pol model as is depicted in Figure 4. To show the effect of the time delay on self-excited pressure oscillations and to determine the steady-state values, external excitation is omitted, because we know now, that its influence is small. On the other side, the computational effort in such an analysis is much smaller.

In accordance with the proposed structure, the van der Pol model in Figure 4 is described by governing equations:


where p˙θ=x˙θ,1+x˙θ,2 is the time derivative of pressure p=x1+x2 at the burning plane, which is delayed for a time delay θ and LPF denotes the action of the low pass operator on the heat release process.

By performing the EL-P procedure as described in Section 3, one obtains the following set of PDEs:


Low pass filter is linear and the filter dynamics is assumed to be much faster than the slowly varying amplitudes and phases in dependence on the slow time scale τ4 (the index 4 in τ4 is retained here for convenience). Because the time delay θ can be assumed small, amplitudes and phases are considered independent of delay (that is Aθ,i(τ4)=Ai(τ4)and Φθ,i(τ4)=Φi(τ4)). Under this assumption, the zero-order solution of Eq. (42) reads:


The low pass filter produces The gain G(ω) and the phase shift φ (ω) at the phase shift ϕ(ω) at the specified frequency ω and affects the input signal in accordance with the following relationship:


By means of zero-order solutions (44), one deduces the following solvability conditions for elimination of secular terms in Eq. (43):


Equations (46) and (47) determine the evolution of slowly varying amplitudes Ai(τ4) and phases Φi(τ4), respectively. By comparing solvability conditions in Eqs. (20) and (22) (with K1 = K2 = 0!), one can conclude that evolution of amplitudes and phases, respectively, is strongly influenced by the transport time delay θ as well as by the transfer function of the low pass filter. This influence is reflected in the change of nontrivial steady-state values. Besides trivial steady-state values A1S=A2S=0 the following nontrivial steady-state values are obtained from Eq. (47):


The obtained results show that steady-state values of first- and second-mode pressure oscillations depend inversely on the values of natural frequencies ωni=ωi0[see Eq. (33)]. This outcome is extremely valuable in prediction of the highest level of the acoustic pressures and is in accordance with experimental observations [8].


7. Implementation of the van der Pol model of combustion processes

To eliminate growing combustion oscillations, the coupling between acoustic waves and the unsteady heat release must be interrupted by applying the control. To reach this goal, passive and active control methods are enforced. The passive control is effective in a limited range of operating conditions and design changes are usually expensive and time consuming. Earlier concepts of active control of combustion instabilities used loudspeakers to increase surface losses in a combustor [1]. The lack of loudspeaker actuators is insufficient robustness for use in industrial applications and prohibitive power requirements at large-scale combustion systems. The recent development of active control seeks to achieve a stable operation regime of a full-scale combustion system. The proposed van der Pol model of combustion processes applies a recent concept of an active feedback control of the heat release by controlling the fuel flow rate. The control function is multiplicative and is realized with the modulation of the injection of the fuel flow. The properties of a good actuation can be reached by the use of fuel valve with a linear response characteristic, large bandwidth and operating range to achieve the effective control of two oscillation modes and to affect limit cycle oscillation. In addition, the fuel valve should have a fast response time, good robustness and durability; however, it should not have hysteresis behavior.

Sensors are generally provided to measure the dynamic quantities, which are related to the combustion oscillations. In the present model, the pressure of the acoustic wave and the fuel flow are measured quantities, which are acquired by means of the pressure transducers and sensors, which measure the light emission of some chemicals in the flame or by chemiluminescence sensors, respectively. The use of pressure transducers is advantageous, because they can be installed outside the area of high temperatures at the burning plane on account of acoustic wave propagation across the entire combustion chamber. Pressure transducers have large bandwidth and a good robustness. The placement of pressure transducers is important, because the location of the possible pressure node must be avoided and on other side, the signal-to-noise ratio must be maximized.

Recent development in the heat release sensing offers the use of diode lasers. Earlier, they had been used for temperature measurements, but recently they can directly measure fuel concentrations. A comprehensive review of this topic is given in [9].

The proposed van der Pol model of combustion processes is based on the grey-box principle. This means that the controller must be designed on system measurements. The instability of the combustion process may have severe consequences, which makes measurement a difficult task. Due to existence of two oscillation modes, which span over a wide frequency range, the controller must be able to control both modes. The combustion process almost always involves time delays caused by fuel convection and acoustic wave propagation. When such a delay becomes larger than an oscillation period, assumptions made in Section 6 are violated and provide a great challenge to the controller design [11].


8. Conclusions

In this chapter, the EL-P method with multiple time scales is successfully applied to the analysis of instability problem in the combustion process with external harmonic excitation, which is modeled by means of coupled van der Pol oscillator. By applying the control of the fuel influx, the asphyxiation of pressure oscillations is reached and studied in the phenomenon of the competitive quenching.

The solvability conditions of the first order in EL-P method are derived, which determines slowly varying amplitudes and phases as well as steady-state values. With external harmonic excitation, both limit cycle oscillations and compound oscillations can be successfully suppressed to the acceptable level by using the multiplicative feedback control.


  1. 1. Dowling A. P., Morgans A. S. “Feedback control of combustion oscillations.” Annu. Rev. Fluid Mech. 37: 151-182, 2005.
  2. 2. Pušenjak R. R., Tičar I., Oblak M. M. “Self-excited oscillations and fuel control of a combustion process in a Rijke tube.” Int. J. Nonlin. Sci. Num. 15(2): 87-106, 2014.
  3. 3. Dunstan W. J. System Identification of Nonlinear Resonant Systems. PhD thesis, University of California, San Diego, 2003.
  4. 4. Sterling J. D. “Nonlinear analysis and modelling of combustion instabilities in a laboratory combustor.” Combust. Sci. Tech. 89: 167-179, 1993.
  5. 5. Peracchio A. A., Proscia W. M. “Nonlinear heat-release/acoustic model for thermoacoustic instability in lean premixed combustors.” ASME J. Eng. Gas Turbines Power. 121(3): 415-421, 1999.
  6. 6. Murray R. M., Jacobson C. A., Casas R., Khibnik A. I., Johnson C. R. Jr., Bitmead R., Perachio A. A., Proscia W. M. “System identification for limit cycling systems: a case study for combustion instabilities.” American Control Conference Proceedings, 1998.
  7. 7. Pušenjak R. R. “Extended Lindstedt-Poincare method for non-stationary resonances of dynamical systems with cubic non-linearities.” J. Sound Vib. 314(1-2): 194-216, 2008.
  8. 8. Landau I. D., Bouziani F., Bitmead R. R., Voda-Besançon A. “Analysis of control relevant coupled nonlinear oscillatory systems.” Eur J. Control. 14(4): 263-282, 2008.
  9. 9. Docquier N., Candel N. “Combustion control and sensors: a review.” Prog. Energy. Combust. Sci . 28: 107-150, 2002.
  10. 10. Pušenjak R. R., Oblak M. M., Tičar I. “Modified Lindstedt-Poincare method with multiple time scales for combination resonance of damped dynamical systems with strong non-linearities.” Int. J. Nonlin. Sci. Num. 11(3): 173-201, 2010.
  11. 11. Cohen J. M., Banaszuk A. “Factors effecting the control of unstable combustors.” J. Prop. Power. 19(5): 811-821, 2003.

Written By

Rudolf R. Pušenjak and Igor Tičar

Submitted: 21 October 2015 Reviewed: 21 March 2016 Published: 05 October 2016