Methane combustion kinetic rate constants determination: an ill-posed inverse problem analysis

ARTIGO

Methane combustion kinetic rate constants determination: an ill-posed inverse problem analysis

Bárbara D. L. Ferreira^I; Jaqueline M. Paulo^I; João P. Braga^I; Rita C. O. Sebastião^I,*; Fabrício J. P. Pujatti^II

^IDepartamento de Química, Instituto de Ciências Exatas, Universidade Federal de Minas Gerais, Av. Antônio Carlos, 6627, 31270-901 Belo Horizonte - MG, Brasil

^IIDepartamento de Engenharia Mecânica, Universidade Federal de Minas Gerais, Av. Antônio Carlos, 6627, 31270-901 Belo Horizonte - MG, Brasil

ABSTRACT

Methane combustion was studied by the Westbrook and Dryer model. This well-established simplified mechanism is very useful in combustion science, for computational effort can be notably reduced. In the inversion procedure to be studied, rate constants are obtained from [CO] concentration data. However, when inherent experimental errors in chemical concentrations are considered, an ill-conditioned inverse problem must be solved for which appropriate mathematical algorithms are needed. A recurrent neural network was chosen due to its numerical stability and robustness. The proposed methodology was compared against Simplex and Levenberg-Marquardt, the most used methods for optimization problems.

Keywords: methane combustion; kinetic rate constants; inverse problem.

INTRODUCTION

The utilization of renewable fuels derived from non-fossil carbon sources, such as Biogas, is rising in Brazil. There are different applications of this fuel for power generation such as gas turbines, boilers, steam generators and internal combustion engines.^1-4 In the last case, the use of Biogas as fuel is probably the most advantageous manner of transforming alternative energy into net power. Therefore, it is important to investigate methane combustion since it represents about 65% of this fuel.¹

Generally, models representing methane combustion use more than 300 reactions to define combustion products, requiring complex computational processing.^5-9 To simplify this problem, the well-known and tested Westbrook and Dryer model (WD-modified) can be used for chemical kinetic calculations.⁸ Simulated and real conditions such as biogas concentration, excess air ratio (λ) and temperature are considered in this simulation.^4,9-11

The use of simplified mechanisms is also very helpful in Computational Fluid Dynamics analysis (CFD) since, for the same reason, the modeling of industrial combustion requires great computational effort.^8,11 CFD analysis becomes an important industrial tool in which chemical reactions are often represented by a mixed-is-burned assumption or by a chemical equilibrium.

Kinetic rate constants of irreversible consecutive methane combustion are calculated from the product chemical concentration. For this problem, the differential equations set characteristics and the usage of experimentally noisy data implies an ill-conditioned inverse problem, which has to be solved by robust mathematical algorithms.^12-14 For the present study, a recurrent neural network was chosen.^15,16 The approach is numerically stable and robust with respect to deviations in the initial conditions or experimental noises.

THEORETICAL METHODOLOGY

Methane combustion kinetic modeling

The WD model consists of three reactions in which oxidation of [CO] to [CO₂] is treated as an irreversible reaction. The mechanism is represented by Equation (1):⁸

Chemical kinetic calculations, using this model to describe hydrocarbon oxidation, proved to be useful in CFD models for conventional combustion in air. This simplified mechanism can also be compared to the Detailed Chemical Kinetic Model (DCKM) by Glarborg and Bentzen, as presented in reference 8. The original WD model adequately predicts [CO₂] and [O₂] concentration from [CH₄] oxidation in flow reactors. Nevertheless, the [CO] concentration is not accurately predicted by this model in its original form. For this reason, the global model (WD) was modified by Andersen to improve this property simulation.⁸ The WD-modified model consists of a system of chemical rates:

r₁ = [CH₄]^0.7[O₂]^0.8

r₂ = [CO][O₂]^0.25[H₂O]^0.5

r₃ = [CO₂][H₂O]^0.5[O₂]^-0.25

and the kinetic model for methane combustion can be represented as Equation 2,

These coupled differential equations have to be solved from the chemical species concentrations to retrieve the kinetic rate constants. Although this is a simplified kinetic model, it can yield useful insights into the combustion process.

Neural network theoretical background

In the neural network approach, an analogy between artificial and biological neurons is established by a mathematical model. The neuron states, u_i, are evaluated considering the weighted contributions, T_ij, of all the neurons, u_j connected to them. To propagate the information, the neurons have to be activated by a monotonically increasing function, f(u_j)^14,15

with I_i(t) representing some external impulse and n the number of neurons in the network. During a learning time, t, information is changed in the network by the activated neurons until a convergence criterion is reached. For this, an energy function is defined as,^16-18

with C_cal and C_exp constituting the calculated and experimental properties, respectively.

The property, C_cal, is calculated based on the activated neurons, f(u_i(t)) and the derivative of the energy function can be written as

with e_j = (C_cal,j - C_exp,j) and equal to the number of experimental data.¹²

From Equation 4, a neural network convergence criterion can be established by imposing the condition,

Equation 5 is thus transformed to

From this, one can note it is necessary to adopt the transfer function as an increasing function of the neuron states, i.e. ∂f/∂u > 0, which will imply .^16,19 In this study, the hyperbolic tangent function was used to activate the neuron states.

In this algorithm, the Hopfield neural network (HNN) provides the solution that minimizes the error function, Equation 4, and consequently best reproduces the experimental property. In the present study, the calculated property is the chemical concentration while the neurons in the network represent the rate constants of the process, which were obtained by integrating Equation 6. A fourth order Runge-Kutta method was used in the neuron state evaluation.^20,21 Starting from an initial guess, the stopping point of the integration process is defined as when the error function reaches the minimum value, i.e. du_i/dt = 0.

RESULTS AND DISCUSSION

In the first part of this work, simulated data and the modified WD mechanism for methane combustion were adopted to test the efficiency of the proposed method. Simulated conditions of oxy-firing combustion, i.e. 28% O₂ and 72% N₂, and stoichiometric combustion at 1600 K, were selected. For this assumption, the methane and oxygen initial molar fraction percent concentrations were 11.18 and 23.02, respectively. Figure 1 shows the highest species concentration calculated in the direct problem using the modified WD model with k₁ >= 881 s^-1, k₂ = 2464 s^-1 and k₃ = 0.5 s^-1.⁸ The CO concentration data proved to be in agreement with the validated model, the Detailed Chemical Kinetic Model - DCKM, as proposed by Andersen.⁸

This model adequately predicts the [O₂], [CO] and [CO₂] levels under fuel-lean and oxy-firing conditions and can be used to describe post ignition fuel-consumption rate in isothermal plug-flow reactors.^4,22 The model also predicts the reactor residence time as the time scale for complete conversion of methane to [CO₂ ]. As can be noted in Figure 1, the time of 10^-3 s is in agreement with the literature.²²

The neural network states are established by the integration of Equation 6, with rate constants used as initial guesses. The activated neuron states correspond to inverted rate constants calculated from the error optimization process. The [CO] molar fraction percent concentration data, from Figure 1, was input to the network together with the initial condition k₁ >= 100 s^-1, k₂ = 1000 s^-1 and k₃ = 1 s^-1. These initial estimates were chosen considering the same magnitude as the rate constants used in the simulation step, where the initial residual error for these constants is 13.001. The inverted rate constants retrieved from this data are given in Table 1. Also, this table also shows the residual error for the concentration fit. The good agreement among the data can be verified in Figure 2, in which all concentrations were calculated using the kinetic constants retrieved by the network, as shown in Table 1.

The HNN inputs are established by initial guesses for the rate constants and the converged neuron states correspond to inverted rate constants. The network time convergence, the most used optimization technique, is related to the quality of this initial condition, and this property can be used to test the reliability of the methodology, since for an exact initial guess and theoretical data, the temporal variation of the activated neuron state has to be null. Using the rate constants provided in Table 1 and the simulated [CO] or [CO₂] concentration data, the retrieved constants do not change, thereby proving the reliability of the computational method. The [CO] (or [CO₂]) concentration should be used as the simulated data in the neural network error function, for its greater sensitivity with respect to the rate constants.

Simulated data with random noises was also used to analyze Hopfield neural network efficiency against Simplex and Levenberg-Marquardt algorithms. The Levenberg-Marquardt method requires matrix inversion and applying this inverted matrix to data presenting errors results in large fluctuation in the results due to the matrix condition. Matrix inversion is not performed for the Simplex algorithm, but this method is very sensitive with respect to experimental error. This is not the case for the HNN approach, which solves the problem using the dynamical character of the neurons. In this work, a study of the influence of error on the inversion procedure was performed considering these three methodologies. A comparison among the methods is shown in Table 2. Random noises of 3, 5, 7 and 10% were added to the simulated [CO] concentrations. Even at a noise level of 10%, the algorithm was stable and, in comparison with the Simplex and Levenberg-Marquardt algorithms, the HNN method had smaller errors in the rate constant calculations and in the residual error of the concentration reproduction, as shown in Table 2. The maximum residual error of the HNN was attained when 10% random noise was added to the simulated data. The Simplex and Levenberg-Marquardt methods cannot handle errors greater than 7%, giving unphysical results for the concentration.

Other tests of the algorithms entail providing different initial conditions to propagate the coupled differential equations. In this problem, initial guesses deviating by up to 30% for k₁ and k_2, and up to 300% for k_3, from the correct values were used and in all neural network inversions, the inverted rate constants produce an error of about 10^-6 in relation to the simulated [CO] concentrations.

For the WD kinetic differential equations, the sensitivity analysis involves an investigation on how changes in the concentration data affect the solution of the problem, the kinetic rate constants. Important information on the inversion process can be obtained by analyzing some aspects of the sensitivity curves. The region in which curves assume larger values will be more stable for obtaining the inverted rate constants and conversely, low sensitivity values correspond to a limited potentiality of the inversion procedure.

Figure 3 depicts the sensitivity curves, S_ij = ∂C_j/∂k, of the [CH₄], [CO], [H₂O] and [CO₂] with respect to k₁, k₂ and k₃ constants. These curves indicate the level of information intrinsic in the data with respect to each parameter. To compute the ∂C_j/∂k_i curves, the slope of a nearby secant line through the points (x - h, f(x - h)) and (x + h, f(x + h)) was determined, i.e. the second order numerical differentiation concept was adopted as . In this sense, a set of coupled differential equations, similar to equation 2, was solved. Notably, in the sensitivity curves with respect to k₁ constant, the [O₂] concentration is preferable compared to other data. However, for the k₂ and k₃ constants, the [CO] and [CO₂] data has greater sensitivity and can be used with the same efficiency in the inversion process.

Spark ignition engines working at 3600 rpm commonly present a biogas residence time of about 10^-3 s.^1,4,22 In this case, the rate constants retrieved in this work may reproduce this experimental condition and can be used for kinetic calculations by modifying only the initial chemical concentration. In this example, the initial molar fraction percent concentrations must be 9.04% [CH₄], 18.08% [O₂] and 1.84% [CO₂], which reproduces fuel characteristics. Since rate constants are the same as those shown in Table 1, the time concentration behavior is similar to that depicted in Figure 2. For experimental calculations, only the [CO] or [CO₂] temporal concentration must be input to the network, further illustrating the applicability of the method.

The WD model yields valuable information for a future more elaborate analysis. For example, a more detailed model will not change the importance of [CO₂] and [CO] concentrations, as predicted by the practical WD model. Also, the Hopfield neural network has shown to be an adequate method for retrieving rate constants from experimental data and not restricted to a specific model.

CONCLUSION

Kinetic rate constants of methane combustion were obtained in the present study from [CO] molar fraction percent concentration data using an artificial neural network. The modified WD kinetic model was adopted for its simplicity and accuracy in [CH₄], [CO] and [CO₂] concentration predictions. Simulated data were first calculated considering real conditions including biogas concentration, excess air ratio and temperature in combustion engines. From these data and a well-known WD kinetic model, the Levenberg-Marquardt and Simplex algorithms were compared against the Hopfield network methodology. The efficiency of these algorithms was investigated using random noises in the simulated data and error in the initial guesses. A neural network is a superior method of retrieving rate constants from simulated data with random noise, as evident from the above considerations.

The HNN is non-restrictive with regard to the system and can be used for any set of differential equations which describe a consecutive reaction mechanism. The algorithm was numerically stable in relation to different initial conditions and also more robust than Levenberg-Maquardt and Simplex algorithms, common optimization methods used to solve this kind of problem, for all the noise levels tested. The numerical stability of the method allows treatment of data with large experimental noise and guarantees efficient convergence of neurons even for initial conditions that deviate from the real estimated parameters.

The [CO] concentration was chosen because of its greater sensitivity values with respect to the three constants. However, the [CO₂] concentration can also be used with the same efficiency. The application of results obtained to spark ignition engines provides a useful prediction method in combustion chamber design and can be adjusted to obtain knock effects, activation energy and exhaust emission concentrations for different pressures, temperatures and excess air ratios.

ACKNOWLEDGEMENTS

The authors would thanks CNPq, Fapemig and PRPq for the financial support.

Recebido em 11/6/12; aceito em 10/9/12; publicado na web em 1/2/2013

• 1. Porpatham, E.; Ramesh, A.; Nagalingam, B.; Fuel 2008,87,1651.
• 2. Bendtsen, A. B.; Glarborg, P.; Dam-Johansen, K.; Combust. Sci. Technol. 2000,151,31.
• 3. Westbrook, C. K.; Dryer, F. L.; Combust. Sci. Technol. 1981,27,31.
• 4. Yin, J.; Su, S.; Yu, X. X.; Weng, Y.; Appl. Energy 2010,87,2102.
• 5. Hautman, D. J.; Dryer, F. L.; Schug, K. P.; Glassman, I.; Combust. Sci. Technol. 1981,25,219.
• 6. Westbrook, C. K.; Dryer, F.; Prog. Energy Combust. Sci. 1984,10,1.
• 7. Westbrook, C. K.; Creighton, J.; Lund, C.; Dryer, F. L.; Chem. Phys. 1977,61,2542.
• 8. Andersen, J.; Rasmussen, C. L.; Giselsson, T.; Glarborg, P.; Energy Fuels 2009,23,1379.
• 9. Gosiewski, K.; Pawlaczyk, A.; Warmuzinski, K.; Jaschik, M.; Chem. Eng. J. 2009,154,9.
• 10. Wang, Y.; Liu, Y.; Cao, Q.; Wang, C.; Che, D.; Energy Fuels 2011,25,3437.
• 11. Orbegoso, E. M. M.; Silva, L. F. F.; Novgorodcev Jr., A. R.; J. Braz. Soc. Mech. Sci. Eng. 2011,33,492.
• 12. Lemes, N. H. T.; Borges, E.; Braga, J. P.; Chemom. Intell. Lab. Syst. 2009,96,84.
• 13. Lemes, N. H. T.; Borges, E.; Braga, J. P.; J. Braz. Chem. Soc. 2007,18,1342.
• 14. Sebastião, R. C. O.; Braga, J. P.; J. Magn. Reson. 2005,177,146.
• 15. Hopfield, J. J.; Proc. Natl. Acad. Sci. U. S. A. 1984,81,3088.
• 16. Hopfield, J. J.; Tank, D. W.; Biol. Cybern 1985,52,141.
• 17. Lemes, N. H. T.; Sebastião, R. C. O.; Braga, J. P.; Inverse Prob. Sci. Eng. 2006,14,581.
• 18. Sebastião, R. C. O.; Pacheco, C. N.; Braga, J. P.; Veloso, D. P.; J. Magn. Reson. 2006,182,22.
• 19. Sebastião, R. C. O.; Lemes, N. H. T.; Virtuoso, L. S.; Braga, J. P.; Chem. Phys. Lett. 2003,378,406.
• 20. Vemuri, V.; Jang, G. S.; J. Franklin Inst. 1992,329,2417.
• 21. Forsythe, G. E.; Malcolm, M. A.; Moler, C. B.; Computer Methods for Mathematical Computations, Englewood Clifs: New Jersey, 1977.
• 22. Karagiannidis, S.; Mantzaras, J.; Boulouchos, K.; Proc. Comb. Inst. 2011,33,3241.

Publication Dates

• Publication in this collection
  18 Mar 2013
• Date of issue
  2013

History

• Received
  11 June 2012
• Accepted
  10 Sept 2012