34 | Izvirni znanstveni članek TEHNIKA - nestabilni termoakustični procesi zgorevanja Datum prejema: 30. julij 2014 ANALI PAZU 4/ 2014/ 1: 34-40 www.anali-pazu.si Modeling and Control of Instabilities in Combustion Processes Modeliranje in upravljanje nestabilnosti v procesih zgorevanja Rudolf Pušenjak 1, *, Maks Oblak 1 Faculty of Industrial Engineering Novo mesto / Šegova ulica 112, SI-8000 Novo mesto E-Mail: rudi.pusenjak@fini.unm.si * Corresponding author: Full Prof. Dr. Rudolf Pušenjak, Kidričeva 5, 9240 Ljutomer, Slovenia; Tel.: +386-2-582-14-06 Abstract: This paper treats the control of nonstationary oscillations of acoustic pressure in the combustion instability process, which appears in the combustion chamber. Two models of nonstationary combustion process control are presented, where first of them is a classic model, which applies the Rayleigh criterium of phase matching between heat-release rate and acoustic pressure. The second model is an alternative van der Pol model exhibiting selfexciting oscillations due to the feedback with a negative damping. The paper outlines the model of fuel control, which is successfully applied for quenching of selfexcited oscillations. Efficiency of fuel control to quench selfexcited oscillations is shown in the phenomenon of competitive quenching, where the influence of various control parameters is explained. Key words: instability of combustion process, selfexcited pressure oscillations, multiplicative feedback control function, EL-P method. Povzetek: Članek obravnava upravljanje nestacionarnih nihanj zvočnega tlaka pri nestabilnem procesu zgorevanja v zgorevalni komori. Predstavljena sta dva modela upravljanja nestacionarnega procesa zgorevanja, pri čemer je najprej predstavljen klasičen model, ki temelji na uporabi Rayleighovega kriterija faznega ujemanja sproščene toplote z zvočnim tlakom. Kot alternativa je predstavljen van der Polov model, v katerem samovzbujena nihanja zvočnega tlaka povzroča povratna zveza z negativnim dušenjem. Članek prikazuje model upravljanja dotoka goriva, s katerim lahko samovzbujena nihanja popolnoma ali vsaj v dobršni meri uspešno zadušimo. Učinkovitost upravljanja dotoka goriva pri dušenju samovzbujenih nihanj je prikazana v pojavu tekmovalnega dušenja, v katerih je pojasnjen vpliv posameznih parametrov upravljanja. Ključne besede: nestabilnost procesa zgorevanja, samovzbujena nihanja tlaka, multiplikativna funkcija povratnozančnega upravljanja, RL-P metoda. | 35 MODELING AND CONTROL OF INSTABILITIES IN COMBUSTION PROCESSES 1. Introduction The complete description of combustion process consists from the set of governing partial differential equations (PDE's), which describe the heat release dynamics and acoustic wave propagation in feedback, hydrodynamics, kinetics of chemical processes, the heat transfer and entropy in a large scale model. The novel research is focussed on elimination of the unstable combustion process with reduction of environmental pollution, in particular, with a reduction of nitrogen oxide emissions. Reduction of environmental pollution can be achieved by a reduction of fuel to air ratio inside the combustion chamber, however, the combustion process becomes unstable due to the selfexcited oscillations of acoustic pressure. The problem of instability can be eliminated by the control of the fuel influx, which is realized by introduction of the multiplicative control function to implement an appropriate modulation of the heat release rate. An efficient control of unstable combustion process can be achieved by using a reduced model instead of a large scale model. By using a reduced model, the feedback effect of acoustic wave in the burning plane on the heat release dynamics is described in the paper. Interaction between acoustic wave propagation and heat release is usually modeled through mass, momentum and energy conservation laws respecting the Rayleigh criterion, which provides an increase of energy of the acoustic wave, when the heat is released in the phase with acoustic pressure and a decrease of the acoustic energy, when the heat is released in the antiphase with acoustic pressure. An alternative reduced model of combustion process is presented in the paper by using two coupled van der Pol oscillators, where selfexcited oscillations of acoustic pressure are produced as a result of negative damping rather than as a result of Rayleigh criterion. The efficiency of the control to assure the quenching of selfexcited oscillations of acoustic pressure is analyzed in the van der Pol model by a study of adjustment of various control parameters. 2. Nomenclature ρ = density of mass u = vector of velocity p = acoustic pressure q = heat release rate γ = coefficient of the specific heat t = time θ = time delay ε = gain of the heat release process ω n , ω ni , (i=1,2) = natural angular frequencies ωi 0 , (i=1,2) = linear angular frequencies ω i , (i=1,2) = nonlinear angular frequencies τ 1, τ 2 = fast time scales τ 3 = slow time scale A i (τ 3 ), (i=1,2) = amplitude of oscillations as function of the slow time scale Ф i (τ 3 ), (i=1,2) = phase angle of oscillations as function of the slow time scale 3. Governing equations of the classic model of instabilities in a combustion process Governing equations of the classic model of instabilities in a combustion process and derivation of Rayleigh criterion is attributed to Culick [1]. In this paper, the detailed derivation of governing equations is omitted, however, the frame in which this derivation is made, is presented for the sake of convenience. Under assumption, that selfexcited oscillations of the acoustic pressure arise mainly due to the nonlinear interaction between the heat source and acoustic wave, while body forces are neglected, the governing equations of mass, momentum and energy conservation for the inviscid and nonconductig fluid are written in the form where ρ is the mass density, p is the acoustic pressure, u stands for the velocity vector, q is heat release rate and γ denotes the coefficient of the specific heat. Separating these quantities into two components, where the first component corresponds to the mean value, denoted by overbar and the second component, denoted by the prime means fluctuation about the mean value, the wave equation of fluctuating component of acoustic pressure p' can be derived: together with the corresponding boundary conditions: where vector n denotes the normal on the boundary surface, directed in the outside of the combustion chamber. Solution of the wave equation is expressed in the form of series, which is composed by products of orthogonal acoustic mode shapes ψ i (r) with time variable amplitudes η i (t): where denotes the mean value of the acoustic pressure. Mode shapes ψ i (r) are solutions of the eigenvalue problem: , (1) , (2) , (3) , (4) , (5) , (6) 36 | ANALI PAZU, 4/ 2014/ 1, str. 34-40 Rudolf PUŠENJAK, Maks OBLAK where k i is the wavenumber of the standing wave, which corresponds to the natural angular frequency ω ni of i th mode shape. To determine time variable amplitudes η i (t), we proceed as follows. First, we insert ansatz (6) into Eq. (4) and multiply by ψ i . From resulting equation then we subtract Eq. (7), multiplied by ansatz (6) and integrate the obtained difference over the volume V of the combustion chamber. By considering orthogonality of acoustic mode shapes, by using Green theorem and by fulfilment of ideal boundary condition f = 0 we finally derive the system of ordinary differential equations for modal coefficients: The obtained system of ODE's (9) describes temporal variations od acoustic pressure inside the combustion chamber in the form of coupled oscillators, which are excited by temporal variations of heat release under integral operator on the right hand side of Eq. (9). When the nonlinear model of temporal variations of heat release is known, then system of ODE's (9) becomes autonomous. For solving autonomous system of oscillators, Culick [1] has proposed the method of Krilov-Bogoliubov (KB method). Instead of KB method, an alternative Linstedt-Poincare method with multiple time scales (EL-P method) is proposed in this paper, which has proven very efficient in solving problems of the control of unstable selfexcited oscillations of acoustic pressure in combustion chamber [3]. To demonstrate the usefullness of EL-P method with multiple time scales in classic model of instabilities in a combustion process, that is, to apply EL -P method on Eqs. (9),(10), we first derive a simplified model of heat release. 3.1. A simplified model of heat release The simplest model of heat release assumes that all heat is released in the fluid in a single point that corresponds to the flame front x f . For simplicity we assume also that only one mode shape dominates in propagation of the acoustic wave. This mode shape is either Helmholtz wave, which spreads uniformly in all directions or longitudinal wave. In both cases we can neglect the spatial variability of acoustic pressure in the area of the flame. Because of this, the volume integral on right hand side of Eq. (9) is simplified and results in expression: where q(t) represents the total heat release rate, that is the heat release rate, which is integrated over the total volume of the combustion chamber. This heat release rate can be expressed in the form: where denotes the instantaneous mass flow rate of unburned mixture of fuel with air: Instantaneous mass flow rate of unburned mixture of fuel with air is determined by motion of the flame front, where ρ m denotes the density of the mixture, S u is empirical turbulent flame speed and A f (t) is the flame area. The symbol ΔH m (x f ,t) represents the heat release per unit of mass of the mixture, which reaches the flame front in the time t. Convection of the mixture from the nozzle to the flame front takes the time of θ seconds and results in delay ΔH m (x f ,t)= ΔH m (x o ,t-θ), where x o indicates the location of the exit of the fuel from the nozzle. Similarly as the total heat release rate q(t) in Eq. (12) is expressed by the sum of the mean value and fluctuation around of mean value, so is the instantaneous mass flow rate in Eq. (13) given by the sum of the mean value and fluctuation around the mean value. By introducing dimensionless quantities , , , where R denotes the combustor radius and where , [2], the dimensionless acoustic velocity can be expressed by means of product of acoustic admittance Y and time derivative as: which is used to describe the fluctuating component of the flame area as [2]: In presented model of heat release, the heat release per unit of mass of the mixture ΔH m is usually expressed through dimensionless quantity , where and HHV is the higher heating value, while denotes the fuel to air ratio at stoichiometric conditions. Using dimensionless quantities, Eq. (11) now can be rewritten in the form: Simplified case of Eq. (16) and thus the simplified model of heat release in a whole is obtained, when fluctuations of mas flow rate of mixture subtended by flame front are neglected, and a linear relationship between the dimensionless heat release per unit of mass , (7) , (8) , (9) , (10) , (11) , (12) , (13) , (14) , (15) , (16) | 37 MODELING AND CONTROL OF INSTABILITIES IN COMBUSTION PROCESSES of the mixture and dimensionless acoustic velocity is assumed: where means the interception point and the slope of the straight line. By introducing Eq. (17) into Eq. (16) and considering Eq. (14), one obtains: Equation (18) is a delayed ordinary differential equation (DODE), where a delay θ appears in the second derivative of pressure amplitude . In Ref. [4] it was shown, that DODE's of Duffing type can be successfuly solved by applying EL-P method. 4. Van der Pol model of instabilities in combustion process An alternate approach to the solving of instability problem in combustion process represents a van der Pol model that is not based on the Rayleigh criterion, in which the instability of the process occurs due to the increase of the energy of the acoustic wave, when the heat is released in phase with the pressure, but on the feedback with a negative damping. Van der Pol model is two-dimensional and is in comparison with the classical model improved with modulation of the fuel influx. Modulation of the fuel influx represents a control action, which allows an efficient quenching of selfexcited pressure oscillations. Model contains two dominant natural frequencies, which may be incommensurate in general. These natural frequencies correspond to two mode shapes in accordance with the experimental findings. Governing equations of the model with control function are generalized van der Pol equations of two modal amplitudes forming the amplitude of pressure oscillations [3]: where the control function is chosen in a following polinomial form: In Eqs. (19) and (20), symbols ω n1 , ω n2 mean incommensurate natural frequencies (in general) and function represents a nonlinear function of heat release. Parameter φ 0 in this function denotes a drift and can be an arbitrary positive or negative constant. Parameter φ 1 is the slope of function q about origin and can be an arbitrary positive constant. A very important term of heat release function is the term , which causes the saturation of heat release and where the parameter φ 3 determines the level of saturation. Parameter ε is a small positive number, which determines the gain of the heat release process. Values of parameters ε, φ 0 , φ 1 and φ 3 are determined by using identification methods to reach the best fit to real experimental data of combustion process in combustion chamber. Parameter K 1 of control function is the amplifying factor, which must have the same sign as drift φ 0 , while the parameter K 2 is a constant, which must be greater than unity, K 2 >1. Block scheme of control of instabilities in combustion process is shown in Fig. 1. Inclusion of low pass filter in block scheme is optional, however, the use of low pass filter dramatically improves properties of the control [3]. 5. Extended Lindstedt-Poincare method with multiple time scales for control of instabilities in combustion process with van der Pol model Extended Lindstedt-Poincare (EL-P) method with multiple time scales is a perturbation method that allows computation of an approximate analytical solution of the given problem of control of instabilities in combustion process for small values of gain of heat release process ε. In EL-P method of twodimensional problem we introduce two fast time scales τ 1 =ω 1 t, τ 2 =ω 2 t, which correspond to the nonlinear frequencies ω 1 and ω 2 , respectively, which may be incommensurate in general and an additional slow time scale τ 3 =εt . By means of introduced time scales, we replace differential operators d/dt and d 2 /dt 2 by operators: , (17) , (18) , (19) (20) Figure 1. Van der Pol model of control of instabilities in combustion process. 38 | ANALI PAZU, 4/ 2014/ 1, str. 34-40 Rudolf PUŠENJAK, Maks OBLAK and convert Eq. (19) into system of nonlinear partial differential equations (PDE's): An approximate analytical solution od the system of PDE's (22) is sought in the form of power series: Nonlinear frequencies ω 1 and ω 2 are also expressed with power series around linear frequencies ω 10 and ω 20 : By introduction of ansatz (23) and power series (24) into Eq. (22) and by equating coefficients of equal powers of parameter ε on both sides, one obtains a system of PDE's, which are all linear when are solved successively. For example, equation of zeroth order, which is obtained by assembling terms at power ε 0 , is as follows: Solution of equation of zeroth order is sought by using ansatz for almost periodic solution: where A i (τ 3 ) are slowly varying amplitudes, which are changed in dependence on the slow time scale and Φ i (τ 3 ) are corresponding slowly varying phase angles. Introduction of the ansatz (26) into Eq. (25) gives an important relationship between linear and nonlinear frequencies: From the relation (27) it follows that linear frequencies are equal to the natural frequencies of coupled van der Pol oscillators. Continuation of the perturbation procedure to collect terms, which appear at power ε 1 , requires the assertion of the solution (26) in the right hand side of PDE of the first order. In this procedure, secular terms appear, which must be eliminated. Elimination of secular terms is achieved by fulfilment of so called solvability conditions, which have in the case of PDE of the first order the following form: Solution of equations (28.a,b) reveals the basic properties, which occur in control of instabilities in combustion process. Equation (28.a) shows that trivial solutions A i (τ 3 )=0 always exist. In the case of nontrivial solutions A i (τ 3 )≠ 0, Eq. (28.a) reduces on equation Ref. [3] it is shown, that phase angles change very little and therefore nonlinear frequencies ω i1 are approximately zero, . It is also shown, that an analytical solution of Eq. (28.b) is obtainable in some special cases [3], when the combustion process is not controlled and thus the control parameters are both equal zero, K 1 =K 2 =0. However, to solve equation in a general case, a numerical integration must be applied. When the slow time scale τ 3 increases towards infinity, τ 3 →∞, asymptotic solutions are obtained, which correspond to the limit cycles of selfexcited oscillations. Asymptotic solutions are solutions of Eq. (28.b), in which we request , or equivalently, solutions of the algebraic equation where mean steady state solutions as τ 3 tends to τ 3 →∞. These solutions are: a) trivial solutions , b) nontrivial solutions and c) combinations of a nontrivial solution with a trivial solution . In continuation of the perturbation procedure we seek an almost periodic solution of the first order: (21.a,b) (22) (23) (24) (25) (26) (27) (28.a,b) (29) . In , | 39 MODELING AND CONTROL OF INSTABILITIES IN COMBUSTION PROCESSES Coefficients and respectively are determined as follows. The ansatz (29) is inserted into PDE of the first order (that is into PDE, which is obtained by assembling terms at ε 1 ). On the right hand side of obtained equation, the secular terms appear, which must also be eliminated. After elimination of secular terms, unknown coefficients are determined by equating like terms on both sides of equation. When all unknown coefficients are determined, the solution (29) is obtained. Both solutions (26) and (29) are then inserted into power series (23). Usually first two terms of the series are sufficient and thus close the procedure. Otherwise, the procedure must be continued with approximation of the second order, etc. 6. Laplace transform solution of combustion process equation with simplified model of heat release Although the application of the EL-P method with multiple time scales to solve delayed ordinary differential equation (DODE) is quit possible, there is no need to do this. Instead of the EL-P method, we can use the Laplace transformation. By introducing the Laplace variable s=σ+iω, (i=√-1), we have the following Laplace transformation of Eq. (18) [5]: where denotes the Laplace transform of η(t). By rearrangement we have: Solution η(t) in the time domain is obtained by using inverse Laplace transformation. Because Laplace transform (31) contains a transcendental function, computation of the inverse is not a simple task and this computation will not be presented here. Nevertheless, the inverse Laplace transform can be conveniently computed in a numerical way by using the inverse Fast Fourier Transformation (FFT). What is more important, the stability analysis can be performed by investigation of the denominator of the Laplace transform (31) or by performing the Bode analysis on the open loop transfer function [2]. 7. Results A comprehensive analysis of the control of instabilities in combustion process using van der Pol model was presented in Ref. [3]. Because of this only the phenomenon of competitive quenching is treated in this paper. The phenomenon of the competitive quenching is important for two reasons. Firstly, the phenomenon of the competitive quenching is characterized by the excitation of one of the resonators, while oscillations of the other resonator are quenched in dependence on initial conditions. Secondly, the competitive quenching occurs in the presence of two incommensurate natural frequencies and thus the general theory of the Section 5 must be applied. The analysis is performed for values of parameters φ 0 =1, φ 1 =178, φ 3 =1.24×10 7 ,ε=1 and for natural angular frequencies of selfexcited oscillations ω n1 =2π×210 s -1 , ω n2 =2π×740 s -1 . The choice of the value of the parameter ε=1 allows that the real time t corresponds to the time scale τ 3 , τ 3 =t. (30) (31) Figure 2. The time response of selfexcited oscillations of the acoustic pressure. a) Without the control, b) with applied control , 40 | ANALI PAZU, 4/ 2014/ 1, str. 34-40 Rudolf PUŠENJAK, Maks OBLAK Oscillations of the acoustic pressure of both resonators wihout the applied control are depicted in Fig. 2,a, while the perfect asphyxiation of selfexcited oscillations due to the activation of the control in the switching time t on =0.105 s is shown in the Fig. 2.b. Slowly varying courses of amplitudes A i (τ 3 )=A i (t), which corresponds to the solution (26) are shown in Fig.2.a,b by red dashed line, while high-frequency oscillations according to the two-term solution (23), that is to the sum are plotted by continuous line. 8. Conclusions Two models of instabilities in the combustion process are presented in the paper. The first model is a classical model of the reduced order, which is based on the mass, momentum and energy conservation laws and utilizes the Rayleigh criterion. Solving of the instability problem of the combustion process by applying the classic model is performed in the paper with using one acoustic oscillator and a simplified heat release function. A much more complex analysis of the combustion process, which includes the control of selfexcited oscillations of acoustic pressure is done in the second model, which utilizes two generalized coupled van der Pol oscillators. The EL-P method with multiple time scales is successfuly extended to the analysis of instability problem in the combustion process in the model of autonomous coupled van der Pol oscillators. In the simplified model of heat release, a straightforward Laplace transformation of the combustion process equation is applied instead. The first order approximation is performed in the EL -P analysis in order to compute the self-excited oscillations in a phenomenon of the competitive quenching. The presented paper shows, that EL-P method can be advantageous extended to analyze the problems of the automatic control. References [1] F. E. C. Culick, Combustion Instabilities in Liquid- Fueled propulsion Systems-An Overview, AGARD CP-450, 1989. [2] A. A. Peracchio, W. M. Proscia, Nonlinear Heat- Release/Acoustic Model for Thermoacoustic Instability in Lean Premixed Combustors, ASME Journal of Engineering for Gas Turbines and Power, Vol. 121, 415--421, 1999. [3] R. R. Pušenjak, I. Tičar, M. M. Oblak, Self-excited oscillations and Fuel Control of a Combustion Process in a Rijke Tube, Int. J. Nonl. Sci. Num. Sim., 15(2), 87--106, 2014. [4] R. Pušenjak, M. Oblak, The Control of Nonlinear Oscillatory Systems with Delay – Upravljanje nelinearnih nihajočih sistemov z zakasnitvami, Anali PAZU, 3(1), 15-24, 2013. [5] R. Pušenjak, M. Oblak, M. Fošner. Osnove modeliranja dinamičnih procesov v logistiki II. Celje: Fakulteta za logistiko, 2013.http:// blend.fl.uni-mb.si.