Integrated Procedure, Using Differential Evolution Optimization of Rate

Feb 20, 2018 - The unaffordable computational costs of detailed mechanisms in full three-dimensional applications have resulted in development of the ...
0 downloads 6 Views 878KB Size
Subscriber access provided by - Access paid by the | UCSB Libraries

Article

An Integrated Procedure, using Differential Evolution Optimization of Rate Parameters, for Design of Small and Accurate Multi-step Global Chemical Mechanisms Alireza Shakeri, and Karim Mazaheri Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.8b00144 • Publication Date (Web): 20 Feb 2018 Downloaded from http://pubs.acs.org on February 26, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

An Integrated Procedure, using Differential Evolution Optimization of Rate Parameters, for Design of Small and Accurate Multi-step Global Chemical Mechanisms Alireza Shakeri1*, Karim Mazaheri1 1

Center of Excellence of Aerospace Systems, Sharif University of Technology, Tehran, Azadi Avenue, Iran *Corresponding

author at: P.O.B 1458833351, Tehran, Iran.

E-mail address: [email protected]

ABSTRACT Three dimensional analysis of combustion chambers in industrial gas turbines suffers from lack of simple and accurate reduced mechanisms for oxidation of hydrocarbon fuels. Here an integrated procedure is introduced based on a differential evolution optimization technique. The procedure is flexible and modular and allows optimization of many rate parameters of a multi-step global mechanism based on many different combustion criteria and inlet or operational conditions. The procedure uses any selected chemical reactor model and any reference combustion mechanism provided. Sample design criteria used here are flame temperature, ignition delay time, and concentration of selected species, especially pollutants. Different low atmospheric pressures and high pressure combustion conditions are considered. Reactor models used here are zero and one dimensional reactors of PSR, PFR, IGNITION and PREMIX models. To show applicability and performance of the procedure, the algorithm is analyzed in detail during four different case studies. The first two cases consider application to methane, using a five-step global mechanism. Case 3 considers application of a 3-sep mechanism for simulation of propane oxidation, and finally kerosene oxidation via a 3-step mechanism to predict its transient behavior is studied in case 4. Many global, temporal and spatial combustion properties are studied, and it is shown that for most off-design conditions studied here the relative error of the most significant performance parameters is below one percent. Key Words: Multi-step Global Mechanisms, Optimization of Rate Parameters, Hydrocarbon Fuel Oxidation.

1. INTRODUCTION The unaffordable computational costs of detailed mechanisms in full three-dimensional applications have resulted in development of the so called reduced mechanisms. There are many proposals for reduction of detailed mechanisms to less expensive reduced mechanisms. This has resulted in introduction of many different mechanisms with different levels of details and comprehensiveness. Selection of insignificant mechanisms or missing the significant reactions will result in non-accurate or even erroneous results. Furthermore, mechanisms which are realistically large and comprehensive include a very large number of huge equations, and these equations have very different time scales, resulting in a stiff system of equations, and the time and space complexity of computations become intolerable for almost all industrial applications. This problem is more severe

1 ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 28

as the hydrocarbon fuel molecule gets larger, since the number of involved species and reactions grows exponentially with the molecule size. Unfortunately almost all liquid fuels used in engines include very large molecules. In fact even 1D simulations of these systems of equations are very expensive. One may sort these reduction mechanisms into five major categories, i.e., global multi-step, lumping, time scale analysis, stiffness reduction and skeletal reduction. Each category include a variety of different methods, and they are partly reviewed recently by authors [1]. The global multi-step mechanisms are probably the most reduced mechanism, have the highest cost advantage and obviously are viewed to have the least accuracy. The other methods usually are based on complex derivations, and result in still huge system of equations to be solved simultaneously [2-5]. Therefore their usage in full industrial applications are still very challenging. In global mechanisms, the rate of each reaction is a function of Arrhenius relation coefficients, i.e., preexponential factor, species concentration exponents, activation temperature, and concentration of each reactant with a power of an unknown parameter. Simplicity of these global multi-step mechanisms has motivated many researchers to optimize their parameters to reach to an acceptable level of accuracy, at least for a selected range of applications [6-9]. This is the main goal and contribution of this article. Many different methods are proposed to find parameters of these rate coefficients. Most previous global mechanisms are designed for accurate production of a single criterion, e.g., laminar flame speed, therefore their prediction of other parameters as species concentration in the products and flame temperature is very weak and unreliable. Many of these global mechanisms are designed using experimental data to accurately predict one of species concentration in the combustion products, e.g., CO2, H2, H2O, CO [26-28]. Initially most people used to use experimental data to find coefficients with the best match [10, 11]. Later, regression analysis [12] and optimization schemes like genetic algorithms [13] were applied to find more suitable coefficients. Many studies have been devoted to make these coefficients a function of combustion parameters to make these mechanisms more suitable for a wider range of conditions. The main combustion parameters used for these improvements were equivalence ratio [12] and the reactants input conditions [13]. Many researchers have developed these global mechanisms exclusively for a certain application or to predict certain combustion properties, e.g., the laminar flame speed [10], the ignition delay time [14, 15], CO concentration and the heat release or the flame temperature [16-18], and NO and CO pollutants concentration [7, 12]. Polifk et al. [13] used a premixed CRM and applied the genetic algorithm to find a two-step and a three-step global mechanism for methane oxidation. Harris et al. [19] used a PSR CRM and the genetic algorithm to develop a 19-step global mechanism for modeling Hydrogen and Air combustion. Later Elliot et al. [20] used multi-objective methods to improve this 19-step global mechanism. Elliot et al. [21] reviewed application of different genetic 2 ACS Paragon Plus Environment

Page 3 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

algorithms in optimization of global mechanism introduced for oxidation of hydrogen and methane. In another work, Elliot et al. [22] also used the genetic algorithm and the PREMIX and PSR CRMs to optimize two global mechanisms developed for aviation fuels. They also used hybrid genetic algorithms to develop a skeletal mechanism for oxidation of aviation fuels [23]. Amorim et al. [24] extensively compared genetic and differential evolution (DE) optimization algorithms to show that DE algorithms are more stable and have a faster convergence. Authors have recently used the DE algorithm to optimize the geometry of an experimental jet stirred reactor to produce less NO and CO pollutants [25]. One may introduce a numerical tool (combination of a hybrid CRM and an optimization algorithm) which may produce different global mechanisms for different regions inside a combustion field, and simultaneously considering many different parameters for accuracy criteria, as part of the optimization cost function (e.g., flame speed and temperature, and ignition delay time). In this way, one only needs to check the accuracy of a global mechanism in the vicinity of its reference point or inside the reference region. To implement this strategy, we need to develop a fast numerical tool to design simple multi-step global mechanisms with optimal accuracy considering many different combustion conditions and accuracy criteria. This target is sought in this article. A few advantages of the new algorithm “Differential Evolution Optimization of Rate Parameters” DORP follows: • Considering a range of combustion conditions and reactants conditions. • Using DE optimization algorithm that converges faster towards the global optimum. • Considering a range of different accuracy criteria as the cost function or performance criteria. • Modular structure which allows simple and easy usage for different fuels, different global mechanisms, different accuracy criteria, and different chemical reactor models (CRMs). Aspects of novelty of this work include: • Different from previous works, the new optimization algorithm, combined with different appropriate CRM models, has no restriction on the number of combustion conditions and parameters and their range of applications. Using combination of the PSR-PFR model, as its reactor model, which is more appropriate for gas turbine and many other industrial combustors. • The scheme is so fast and easily applicable that instead of designing a mechanism for a wide range of parameters, one may produce many different mechanisms each suitable for a small region close to the reference point.

3 ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 28

• The algorithm allows the user to select simultaneously many different performance or accuracy criteria for the cost function, including flame temperature, pollutants concentration, ignition delay time, and etc. The second section of the article describes CRMs used in these simulations. In the third section the DE optimization algorithm is introduced. Details of the DORP procedure and its implementation are explained in the fourth section. In the fifth section the algorithm is applied to many different cases to generate a valid model, and the accuracy of the model is assessed.

2. CHEMICAL REACTION MODELS (CRM) To assess accuracy and reliability of a global mechanism, either experimental data or reliable computational tools shall be used. Here we use several zero and one dimensional chemical reaction models. The reaction models used here are functions PSR, PFR, IGNITION and PREMIX of CHEMKIN [29-31]. PSR and PFR are used when prediction of NO and CO concentrations are sought, and IGNITION is used to predict ignition delay time. Details of these models are described below [32]. For a better assessment of accuracy, the PREMIX model is also used. 2.1.

PSR and PSR-PFR Models To predict species concentration and flame temperature, both PSR and PFR models are extensively used by

different people [2-7]. The combination of these two models is motivated by our understanding that the flame region is well simulated by the PSR model, and the region after the flame is well simulated by the PFR model [12]. The combustion conditions, i.e., the equivalence ratio, and reactants initial temperature and the reactor pressure are part of known input parameters. To describe how the reaction rates of the base elementary reactions are calculated, consider the general form of an elementary reaction: 

 ,



́ ∙  ,  = 1, … ,   ́ ∙    



(1)



in which r is the reaction number, and R is the total number of elementary reactions.  is the mole fraction of

species k, and ́  and ́  are the stoichiometric coefficients of species k in the respectively forward and

backward directions of reaction r. These coefficients are all given by a file of reactions data. The information regarding the third body efficiency to find the contribution of each species in the corresponding reaction is also given in this file. This file also includes values of all coefficients of the Arrhenius relation (9)

4 ACS Paragon Plus Environment

Page 5 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

In this study we have neglected heat loss, therefore we do not need to estimate it. To find the stoichiometric coefficients in relation (2) we consider lean oxidation of a hydrocarbon fuel with air as given by:   +

1 1−∅ 3.76  + " #$% + 3.76*% + → $% + % $ + . 0  + " $% +  + " *% ∅ 4 2 ∅ 4 ∅ 4

(2)

Outputs of the PSR model are species mass fraction 2 , the adiabatic flame temperature T, and mixture density ρ. The known (input) data for the PFR model are value of variable ρ,4, 5 and 2 at x=0, which are produced

by PSR function, in addition to combustor diameter, mass flux, and combustor length. The function produces the spatial distribution of ρ, 5 and 2 , and also pressure p and mixture velocity, u. Other inputs are similar to the PSR model. 2.2.

IGNITION Model Many people have used the IGNITION function of CHEMKIN to analyze transient characteristics of a

reaction, e.g., the ignition delay time, and this function is validated by experimental data [14-15]. We use this function to design a three-step global mechanism for oxidation of kerosene. The IGNITION model input data include the combustor volume V, the initial temperature T0, the pressure P, and the reactants equivalence ratio ɸ. The model predicts the temporal variations of temperature T(t), and each

specie mass fraction 2 #6+.

2.3.

PREMIX Model This reactor models a combustion channel with a variable cross sectional area A(x) and length of L. The

reactants enter at x=0, and products leave the combustor at x=L. Input data for this model are the reactants temperature, pressure, and equivalence ratio. The model computes the spatial variation of species concentration and temperature.

3. DIFFERENTIAL EVOLUTION (DE) OPTIMIZATION Differential Evolution (DE) is one of the most efficient evolutionary algorithm [35, 36], developed based on perturbing the population member of each generation using randomly selected, and the scaled difference of a distinct pair of members. The alteration of the population is achieved through Mutation and Crossover operators. An optimization procedure includes three main elements: design variables, a cost function, and constraints. DE optimization process starts from random generation of Np members in a D-dimensional space (D: number of 5 ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 28

design variables). Each generation consists of Np members, and each member represents a vector with D components or gens. The number of members of each generation (Np) is selected by us. Each member of a generation is considered as a vector:

% : 7,8 = 97,8 , 7,8 , … . , 7,8 ;, < = 1,2, … *=

(3)

where G is the generation number. The allowable maximum and minimum values for design variables are set as optimization constraints: x?@A ≤ x@, ≤ x?DE , j = 1,2, … D B

B

B

(4)

Np randomized members are assessed by the cost function, to find the cost value of each member. The cost function describes what we want to optimize and is a D variable function. The next step is the mutation step which will be repeated Np times for each generation. Following Eq.1, the Mutation operator generates an intermediate

population of vectors or members (HI7,8 ) based on addition of a scaled difference of two randomly-selected vectors, I%.8 and IJ,8 , over the base vector, I ,8 . Here, F is a scalar which is used for scale differencing and is usually a real number in the range of [0, 2].

HI7,8 = I ,8 + K#I%,8 -IJ.8 +

(5)

Next, according to Eq. 28, the mutant vector is discretely crossed with the selected base vector I ,8 in the generation to generate the trial vector (M NI7,8 ). The Crossover operator is characterized by a parameter,

current

termed crossover rate ( ), which governs how many components from the mutant vector are inherited by the trial vector.

M NI7,8 #O+ = P

HI7,8 #O+ CO2 CO2=>CO+0.5O2 N2+O2 =>2NO N2+O2 =>2NO

A1[CH4]B1[O2]C1exp(-D1/Ru/T) A2 [CO]B2[O2]C2exp(-D2/Ru/T) A3[CO2]B3exp(-D3/Ru/T) A4 [O2]B4[CO]C4exp(-D4/Ru/T) A5 T-0.5 [N2] [O2] exp(-D5/Ru/T)

* This mechanism contains seventeen unknowns

4.5. Cost Functions The target function or cost function is the main property that we like to be minimized. Since the main disadvantage of global mechanisms is their lack of acceptable accuracy, here we use different accuracy criteria which are appropriate for different industrial applications. The DORP procedure allows very easily to select different accuracy criteria. For each case studied in section 5 and based on the CRM used, we select one of the cost functions described in equation (10-13). ∅Š

vwx xˆ‰ #∅+ − 2s‡ #∅+†+ T∅ {€,‚ƒ„ = … #†2s‡

(10)

∅‹

9 ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

{€,‚ƒ„Œ‚„ =

xˆ‰ xˆ‰ 1  †2s‡ #+ − 2s‡ #+† †2{‡ #+ − 2{‡ #+† †4’ … # + + xˆ‰,ˆ xˆ‰,ˆ Ž ‘ 2s‡ 2{‡

{“€”,‚ƒ„Œ‚„ =

vwx

vwx

xˆ‰

xˆ‰ vwx xˆ‰ 1  †4’ #+ − 4’ #+† †2{‡ #+ − 2{‡ #+† … # + +T xˆ‰,ˆ Ž ‘ 4’xˆ‰,ˆ 2{‡ vwx

#+ − 4’

vwx

4’xˆ‰,ˆ

#+† +T

xˆ‰ {‹•€Š•,–8s–—–‡s = †57v˜ − 57v˜ † vwx

Page 10 of 28

(11)

(12) (13)

Here, ɸ is the equivalence ratio, x is the axial coordinate along the combustor length, bas denotes the base or reference mechanism, glb denotes the global mechanism, and τign is the ignition delay time. More elaborate cost functions could be designed to involve non-uniform weight function for different points of range of combustion condition, or even multi-point or multi-objective and robust optimization techniques could easily be applied, but we found these fairly simple cost functions effective enough.

4.6. Constraints Each design variable in the DORP procedure need to a specific given range. This range is selected based on our experience, and by trial and error, so that the optimum point does not belong to the boundaries. These ranges are given in each case study in section 5.

4.7. DE Inputs and Outputs The DE routine has many input data. The main input are the cost functions (given in equations (10-13)), the initial value of the design variables (given in tables 2-3), and their possible range of variation. Each time the DE routine generates a new member of a new generation of people, a new set of reaction rate coefficients is proposed, which is written in a text file called GlblMech.inp, and is sent as an input for CHEMKIN. 4.8. Implementation to CFD computations The above procedure may be used in 3D CFD computations very easily. Authors have used a reduced multistep global mechanism to analyze NO production in a 3D JSR reactor [25, 53]. Two proposals are given here. Many people (e.g. [51]) have tried decomposing of a full industrial reactor to a network of CRMs, and using global mechanisms with fixed coefficient for each CRM. In fact there are semi-automatic procedures (e.g. [52]) which start with a rough CFD computation to be used in design of an appropriate network of CRMs. The contribution of this work is that now we may use those initial CFD results to find range of variables and select appropriate cost functions in different regions and produce a distinctive mechanism for each CRM. This will obviously enhance significantly the accuracy, while will not affect the computational cost. The optimized network could be used for 10 ACS Paragon Plus Environment

Page 11 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

many computations. Another idea is that we use results of a rough 3D CFD computation to disintegrate the computational domain to several zones, each matched with a different optimized global mechanism. This is in many ways more efficient than the common practice of using one global mechanism in the whole domain.

5. RESULTS To study applicability and accuracy of the DORP procedure, four cases are studied here. In each case, a hydrocarbon fuel oxidation is modeled by both a reference mechanism and an optimized global mechanism, and the difference of these two and many other reduced mechanisms are compared. Table 4 shows details of these cases, and results of analysis of each case are given in a following separate section.

Table 4 Details of the different case studies for implementing of DORP Approach. 1 2 3

Case Fuel Base mechanism Global mechanism CRM Range/Inlet operating Conditions

ɸ T0(K) P (atm) Combustion optimized parameters Cost function relationship

Operating conditions

4

CH4

CH4

C3H8

C10H20

GRI-3.0[38] 5 step[12] PSR

GRI-3.0[38] 5 step[12] PSR-PFR

Range of ɸ

Inlet Conditions

0.4-1.0 600 1.0, 30.0 XNO Eq. 10

0.70, 0.66 368, 573 6.12, 6.28 XNO, XCO, Tf Eq. 11

SanDiego[40] 3 step PSR-PFR Inlet Conditions 1.0, 0.5, 0.5 300,600,600 1, 10., 30. XCO, Tf Eq. 12

Wang[2] 3 step IGNITION Inlet Conditions 1.0 1500 1.0 τign Eq. 13

5.1. Case 1: M5 mechanism for methane oxidation with air for a range of ɸ Since methane is probably the most important hydrocarbon fuel and is extensively used in the industry, many reduced and global mechanisms are proposed for simulation of its oxidation with air. Regarding environmental concerns, those reduced and global mechanisms with capability of predicting pollutant concentrations, i.e. NO and CO, are of main concern [7,12]. Here we use the DORP procedure to find its capability to optimize some global mechanisms proposed by other people. In the first case of our study, we use the global mechanism M5 (see table 3) proposed by Nicol [12] for air oxidation of methane. The base conditions used by Nicol are (τ=0.5 ms, T0=650 K, P=1 atm, ɸ=0.55). The seventeen design variables (rate coefficients) are given in table 3. The CRM used here is a PSR. There are two other wellknown reduced mechanisms. One is a five-step global mechanism introduced by Mallampalli [3] which is based on 11 ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

time scale analysis, under pressure of 1 and 30 atms, and for a range of equivalence ratios between 0.4 and 1.0. The other one is a skeletal mechanism proposed by KF [39]. Mallampalli compared his results with those of the full mechanism of GRI-2.11 [3]. The DORP procedure is used to optimize the M5 global mechanism, and we will compare its performance for a range of input data and combustion conditions with the other two reduced mechanisms and two detailed mechanism of GRi-2.0 and GRI-3.0. Using DE optimization engine, all seventeen design variables are optimized. The range for all design parameters are given in table 5. These are found based on our experience and trial and error. The convergence history versus generation number for two pressure conditions of one and 30 atms (case 1 in table 4) is shown in figure 2. The cost function is based on accuracy of CO and NO concentration and the flame temperature, in a PSR reactor. In less than 60 generations convergence is achieved. Each generation involves about 30-40 persons, i.e., the CRM is called about 2000 times. In a machine with an AMD eight-core 4.0GHz CPU it takes about 30 minutes to achieve a converged solution. The resulted optimum values for all seventeen design variables for both operating pressures are given in table 6. 1

1 atm 30 atm

0.8

Convergence history of case 1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 28

0.6

0.4

0.2

0

20

40 Generations

60

80

Figure2 Convergence history of the optimized rate coefficients for M5 mechanism for case 1 of table 5 obtained from DORP approach using a PSR configuration.

Table 5 Ranges of all 17 design variables used in DORP approach for M5 mechanism for case 1 (units are moles, cm3, sec, Kelvin). Reaction A B C D number 1 2 3 4 5

min 1016.9 1020.5 1013.8 1034.8 1016.0

max 1018.8 1021.9 1014.9 1035.9 1016.9

min 1.10 1.10 1.00 3.80 ---

max 1.90 1.90 1.01 4.90 ---

min 0.10 1.10 --0.20 ---

max 0.90 1.90 --0.90 ---

12 ACS Paragon Plus Environment

min 0.1 e+5 0.1 e+5 1.0 e+5 1.0 e+5 1.1 e+5

max 0.9e+5 0.9e+5 1.9e+5 1.9e+5 1.6e+5

Page 13 of 28

Table 6 Optimized values of all 17 rate coefficients for M5 mechanism under operating conditions of case 1 (table 5) using DORP algorithm (units are moles, cm3, sec, Kelvin). No. Reaction Description Reaction rates with optimized coefficients 1 2 3 4 5

CH4+1.5O2 =>CO+2H2O CO+0.5O2 =>CO2 CO2=>CO+0.5O2 N2+O2 =>2NO N2+O2 =>2NO

P=1 atm (T0=600 K, ɸ=0.4-1.0)

P=30 atm (T0=600 K, ɸ=0.6-1.0)

1018.066[CH4]1.476[O2]0.5exp(-43918/Ru/T)

1018.392[CH4]1.710[O2]0.596exp(-12277/Ru/T) 1020.975[CO]1.572[O2]1.787exp(-30594/Ru/T) 1014.746 [CO2]1.004exp(-163184/Ru/T) 1034.985[O2]3.821[CO]0.246exp(-189571/Ru/T) 1016.353 T-0.5 [N2] [O2] exp(-141727/Ru/T)

1021.434[CO]1.791[O2]1.4exp(-24924/Ru/T) 1014.253[CO2]1.01exp(-118103/Ru/T) 1035.179[O2]4.185[CO]0.751exp(-93012/Ru/T) 1016.353T-0.5 [N2] [O2] exp(-129334/Ru/T)

Comparison of results with other mechanisms

The performance of the DORP optimization procedure applied to M5 global mechanism at conditions of case 1 in table 4 for a range of different equivalence ratios 0.4-1.0 and two reduced mechanisms of Mallampalli [3] and KF [39], and two detailed mechanisms of GRI-2.11 [41] and GRI-3.0 [38] are compared in this section. Note that all four mechanisms are computationally much more expensive than M5 mechanism, and the optimized M5 is at most expected to approximate to their accuracy. We also assume that GRI-3.0, as our base or reference mechanism, is the most accurate one. Figure 3 and 4 compares respectively NO and CO concentration for different equivalence ratios for all mechanisms. The NO predictions are very well, and even better than Mallampalli [3] and GRI-2.11. Errors for CO concentrations are slightly higher.

500

1500 Methane-air combustion T0=600 K, P= 30 atm, τ=2 ms

Methane-air combustion T0=600 K, P= 1 atm, τ=2 ms

400

1000

300 5 step present GRI-2 [41] GRI-3 [38] KF [39] 5 step Mall. [3]

200

5 step present GRI-2 [41] GRI-3 [38] KF [39] 5 step Mall.[3]

NO(ppmv)

NO (ppmv)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

500

100

0

0.4

0.5

0.6 0.7 0.8 Equivalence ratio (φ)

0.9

0

1

0.4

0.5

0.6 0.7 0.8 Equivalence ratio(φ)

0.9

1

(a) (b) Figure 3 NO concentration versus ɸ based on the optimized rate parameters for the M5 mechanism obtained using DORP approach by a

13 ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

PSR reactor, in comparison with GRI-2.11[41], GRI-3.0 [38], KF [39], and Mallampalli [3], τ= 2.0 ms, T0=600 K, (a) P= 1.0 atm and (b) P= 30.0 atms

Another important combustion performance parameter is the flame temperature. The flame temperature and concentration of other species CO2, H2O and O2 are also predicted by these five mechanisms and are compared in figures 5-6. One notes that although we have only used NO concentration in our cost function, still most performance parameters and combustion products are predicted by an acceptable level of accuracy. This confirms that even a simple cost function is very effective in simulation optimization, since most combustion parameters are inter-related in many different ways, and matching NO concentration makes most other performance criteria fairly satisfied. 1

3

Methane-air combustion T0=600 K, P= 30 atm, τ=2 ms

Methane-air combustion T0=600 K, P= 1 atm, τ=2 ms

0.8

2

CO(%mole fraction)

CO (%mole fraction)

2.5

5 step present GRI-2 [41] GRI-3 [38] KF [39] 5 step Mall. [3]

1.5

1

0 0.4

0.6 5 step present GRI-2 [41] GRI-3 [38] KF [39] 5 step Mall.[3]

0.4

0.2

0.5

0.5

0.6 0.7 0.8 Equivalence ratio (φ)

0.9

0 0.6

1

0.7 0.8 Equivalence ratio(φ)

(a)

0.9

(b)

Figure 4 CO concentration versus ɸ based on the optimized rate parameters for M5 mechanism obtained using DORP approach by a PSR reactor, in comparison with GRI-2.11[41], GRI-3.0 [38], Skeletal KF[39], and Mallampalli [3], τ= 2.0 ms, T0=600 K, (a) P= 1.00 atm and (b) P= 30.0 atms.

2600

0.1 Methane-air combustion T0=600 K, P= 30 atm, τ=2 ms

0.09

Methane-air combustion T0=600 K, P= 30 atm, τ=2 ms

CO2(mole fraction)

2400 Flame temperature(K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 28

2200

5 step present GRI-2 [41] GRI-3 [38] KF [39] 5 step Mall.[3]

2000

1800 0.6

0.7

0.8 Equivalence ratio(φ)

0.9

0.08

0.07

5 step present GRI-2 [41] GRI-3 [38] KF [39]

0.06

0.05

0.04 0.6

1

14 ACS Paragon Plus Environment

0.7

0.8 Equivalence ratio(φ)

0.9

1

Page 15 of 28

(a)

(b)

Figure 5 (a) Flame temperature and (b) CO2 concentration versus ɸ based on the optimized rate parameters for the M5 mechanism obtained using DORP approach in a PSR reactor in comparison with GRI-2.11 [41], GRI-3.0 [38], and KF[39], τ= 2.0 ms, T0=600 K, P= 30.0 atm.

0.2

0.1 Methane-air combustion T0=600 K, P= 30 atm, τ=2 ms

Methane-air combustion T0=600 K, P= 30 atm, τ=2 ms

0.18

0.08

O2(mole fraction)

H2O(mole fraction)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

0.16

0.14

5 step present GRI-2 [41] GRI-3 [38] KF [39]

0.12

0.1 0.6

0.7

0.8 Equivalence ratio(φ)

0.9

5 step present GRI-2 [41] GRI-3 [38] KF [39]

0.06

0.04

0.02

0 0.6

1

0.7

0.8 Equivalence ratio(φ)

0.9

1

(a) (b) Figure 6 (a) H2O and (b) O2 concentration versus ɸ based on the optimized rate parameters for the M5 mechanism obtained using DORP approach in a PSR reactor, in comparison with GRI-2.11 [41], GRI-3.0 [38], and KF[39], τ= 2.0 ms, T0=600 K, P= 30.0 atm.

Accuracy assessment

To have a quantitative assessment of accuracy for these comparisons, the L2 norm of error as defined by equation (14) is used here. s

1 7 − 7 % } = ™ . 0 š * ˆ

‘.›

(14)

7

N is the number of data in comparison, yi is the ith data produced by the reference mechanism, and xi is the data from the optimized global mechanism. Here we use comparison points selected by Mallampalli [3], i.e., four points of equivalence ratios of 0.4, 0.6, 0.85, and 1.0 at a pressure of one atmosphere, and three points of 0.6, 0.85 and 1.0 at 30.0 atmospheres. The performance of M5 optimized mechanism and the reduced mechanism of KF [39] and the detailed mechanism of GRI-3.0 [38] are quantitatively compared. Also, results of KF [39] and GRI-2 [41] are partially compared with each other. Summary of results are given in table 7. One observes that the average relative error for NO predictions is about 12% at pressure of one atmosphere, and about 3% at 30 atmospheres. Similar relative errors for CO are respectively about 6% and 4%. The relative errors for flame temperature, and concentration of CO2 and H2O species is much better, and in fact the difference of results with those of GRI-3.0 in these parameters 15 ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 28

are surprisingly negligible. This shows that the optimized M5 mechanism is a very simple and accurate model for simulation of methane oxidation with air. Therefore the DORP procedure is successful to satisfy our accuracy requirements.

Table 7. Quantitative comparison of results of optimized M5 with results of GRI-3.0 [38] and KF [39] of case 1. (All numbers are in percent) NO

Mechanisms

1 atm 3.12 3.19 11.3

M5 v.s. GRI-3 M5 v.s. KF Mall [3] v.s. GRI-2

30 atm 3.12 3.23 1.41

CO

1 atm 3.33 3.34 3.0

30 atm 4.1 4.1 4.7

T

0.53 0.55

CO2

H2O

30 atm 0.76 0.85 0.77 0.87

O2

0.62 0.62

5.2. Case2: M5 mechanism for methane oxidation with air based on inlet operating condition Here we use DORP procedure to optimize seventeen rate coefficients of M5 mechanism introduced in table 3. This case is very similar to the first case. The differences are: • The chemical reactor is a one-dimensional PSR-PFR. • The base or reference mechanism is GRI-3.0 [38] (our main reference for error assessment). • There are two inlet or combustion conditions considered here (based on results given in [7]): 

Inlet condition 1: P=6.12 atms, T0=368 K, ɸ0=0.7, τ=2.5 ms.



Inlet condition 2: P=6.28 atms, T0=573 K, ɸ0=0.66, τ=4.0 ms.

• The cost function (equation (11)) is based on relative errors, and include flame temperature as well as NO and CO concentrations. • Results of optimized M5 are compared with GRI-3.0 [38], KF [39] and eight-step global

mechanism of Novosselov [7]. The possible range of all seventeen design variables are selected based on our experience, and trial and error. The reference mechanism is again GRI-3.0. Each generation has about 30-40 members, and complete convergence is achieved in less than 80 generations. The convergence history for the cost function, based on relative error of NO and CO concentrations and flame temperature indicates that the PSR+PFR chemical reactor is called about 3000 times, which corresponds to a CPU time of about 40 minutes on an AMD eight-core 4.0 GHz CPU.

Comparison of results with other mechanisms

Figures 7-9 compare results of the optimized M5 mechanism for both inlet conditions of case 2, with those of detailed mechanism of GRI-3.0 [38], the skeletal mechanism of KF [39], and the eight step global mechanism of 16 ACS Paragon Plus Environment

Page 17 of 28

Novosselov [7]. In these figures, the spatial distribution of NO and CO concentration and the flame temperature along the PFR reactor are shown. One observes that the optimized M5 mechanism follows very closely the GRI-3.0 full mechanism for all parameters and conditions. The enhanced accuracy respect to the first case is probably due to that we do not deal with a range of different inlet parameters, and also the better cost function has helped to uniformly decrease the error in performance parameters. The higher relative error of Novosselov [7] with respect to GRI-3.0 is due to the fact that this mechanism, although more complex than our M5 mechanism, is designed for a range of inlet parameters, while the optimized M5 is designed for only the inlet conditions in case 2. 0.15

18

Methane-air combustion T0=368 K,P0=6.12 atm, φ0=0.7

Methane-air combustion T0=368 K,P0=6.12 atm, φ0=0.7

CO mole fraction(%)

NO mole fraction (ppmv)

16

5 step present GRI-3.0 [38] KF [39] 8 step Novosselov [7]

14

12

0.1 5 step present GRI-3.0 [38] KF [39] 8 step Novosselov [7]

0.05

10

8

0 0

0.2

0.4 0.6 0.8 Distance along the PSR-PFR (cm)

1

0

0.1 0.2 0.3 Distance along the PSR-PFR (cm)

0.4

(a) (b) Figure 7 (a) NO (b) CO concentration along the PFR reactor based on the optimized M5 mechanisms obtained from DORP, compared with GRI-3.0 [38], skeletal KF[39], and Novosselov [7.], at condition 1: T0=368 K, P0= 6.12 atm, ɸ0=0.7. 44

0.1 Methane-air combustion T0=573 K, P0=6.28 atm, φ0=0.66

42

Methane-air combustion T0=573 K, P0=6.28 atm, φ0=0.66

0.08 40

CO mole fraction (%)

NO mole fraction (ppmv)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

38

36

34

5 step present GRI-3.0 [38] KF [39] 8 step Novosselov [7]

32

30

0

0.2

0.4 0.6 0.8 Distance along the PSR-PFR (cm)

0.06 5 step present GRI-3.0 [38] KF [39] 8 step Novosselov [7]

0.04

0.02

0

1

0

0.2

0.4 0.6 0.8 Distance along the PSR-PFR (cm)

1

(a) (b) Figure 8 (a) NO (b) CO concentration along the PFR reactor based on the optimized M5 mechanisms obtained from DORP, compared with GRI-3.0 [38], skeletal KF39], and Novosselov [7], at condition 2: T0=573 K, P0= 6.28 atm, ɸ0=0.66.

17 ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

2020

1940

Methane-air combustion T0=573 K, P0=6.28 atm, φ0=0.66

Methane-air combustion T0=368 K,P0=6.12 atm, φ0=0.7

1920 Temperature (K)

2000

Temperature (K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 28

1900

1980

5 step present GRI-3.0 [38] KF [39] 8 step Novosselov [7]

1880

1860

0

0.1

0.2 0.3 0.4 Distance along the PSR-PFR (cm)

5 step present GRI-3.0 [38] KF [39] 8 step Novosselov [7]

0.5

1960

0.6

0

0.2

0.4 0.6 0.8 Distance along the PSR-PFR (cm)

1

(a) (b) Figure 9 Temperature profiles along the PFR reactor based on the optimized M5 mechanisms obtained from DORP, compared with GRI-3.0 [38], skeletal KF[39], and Novosselov [7], for (a) T0=368 K, P0= 6.12 atm, ɸ0=0.70, (b) T0=573 K, P0= 6.28 atm, ɸ0=0.66 .

Accuracy assessment

To have a quantitative comparison, again we use the L2 norm of the relative errors given by equation (14). Results of the optimized M5 mechanism are compared with those of GRI-3.0 [38] and KF [39], for both inlet conditions of case 2, and results are shown in table 8. Results show that the DORP procedure has optimized the M5 mechanism easily for both inlet conditions of case 2, and although the M5 mechanism is much simpler and computationally less expensive than GRI-3.0 or KF mechanisms, in terms of performance in prediction of NO concentration or flame temperature, its difference with GRI-3.0 is not significant. To show the advantage of M5 to Novosselov [7], the relative difference of Novosselov mechanism with GRI-3.0 is also given in table 8. We have shown in one of our previous works [25] that the resulting 5-step mechanism designed here using a combined PSR-PFR CRM, is applicable in analysis and optimization of a three-dimensional JSR combustor. In that article, combustion of a turbulent premixed methane-air mixture is analyzed to find concentration of many major and minor species, and the flame structure. Of course, more validation is necessary for better assessment of this procedure. Table 8. Quantitative comparison of results of optimized M5 with results of GRI-3.0 [38] and KF [39] for both conditions of case 2 in table 5. (All numbers are in percent.) Mechanisms

NO 0.11 0.36 7.3

M5 in comparison with GRI-3[38] M5 in comparison with KF [39] Novosselov [7] in comparison with GRI-3 [38]

Condition 1 CO T 9.7 0.08 11.4 0.08 5.5 0.17

18 ACS Paragon Plus Environment

NO 0.16 0.53 1.98

Condition 2 CO T 1.95 0.14 1.96 0.14 5.3 0.17

Page 19 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

5.3. Case3: P3 mechanism for propane oxidation with air based on inlet operating conditions Propane is one of the other very popular fuels both for home and industrial applications. In this case we study how the DORP procedure may be implemented for optimization of rate coefficients of a three-step mechanism introduced in table 2. This mechanism is proposed by us, based on Nicol [12]. •

The chemical reactor for optimization procedure is a one-dimensional PSR-PFR.



The base or reference mechanism is San Diego [40].



There are three inlet or combustion conditions considered here:



-

Inlet condition 1: P=1.0 atms, T0=300 K, ɸ0=1.0, τ=2.0 ms.

-

Inlet condition 2: P=10. atms, T0=600 K, ɸ0=0.5, τ=2.0 ms.

-

Inlet condition 3: P=30. atms, T0=600 K, ɸ0=0.5, τ=2.0 ms.

The cost function (equation (34)) is based on relative errors, and include flame temperature as well as CO concentration.



Results of optimized P3 are compared with San Diego [40].

We generate three different P3 mechanisms, each optimized for one inlet condition. Each optimized mechanism is independently compared with the reference mechanism in different CRMs, i.e., PSR, PSR-PFR, and PREMIX reactors, and for different residence times. The premixed reactor model includes non-homogenous diffusive effects.

The possible range of all eleven design variables are selected based on our experience, and trial and error. The reference mechanism is San Diego [40]. Details of the DE optimization routine is similar to previous cases. Complete convergence is achieved in less than 30 generations. In this optimization the PSR+PFR chemical reactor is called about 800 times, which corresponds to a CPU time of about 10 minutes on an AMD eight-core 4.0 GHz CPU. Results

Figure 10 compares results of the optimized P3 mechanism for inlet condition 1 of case 4, with those of detailed mechanism of San Diego [40]. In these figures, the spatial distribution of flame temperature and CO and CO2 concentrations along the PFR reactor are shown. One observes that the optimized P3 mechanism follows very closely the detailed mechanism of San Diego [40]. The relative error for CO and CO2 concentrations and the flame temperature are respectively 1.85%, 0.62% and 0.5%.

19 ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

2300

0.04

0.11

P=1.0 atm, T0=300 K, φ0=1.0 Propane-air combustion

Propane-air combustion P=1 atm, T0=300 K, φ=1

2200

2150

0.03

0.1

0.02

0.09

CO2 mole fraction

CO mole fraction

2250

Flame temperature(K)

2100 Line: Dan Diego mechanism [40] Symbol: 3 step mechanism

2050

0

0.2

0.4 0.6 Distance along PSR-PFR(cm)

Lines: Sandiego mechanism Symbols: 3-step mechanism

0.8

0.01

1

0.2

0.4

0.6

0.8

1

0.08

Distance along PSR-PFR(cm)

(a)

(b)

Figure 10 (a) Flame temperature (b) CO and CO2 concentration along the PSR-PFR based on optimized P3 mechanism (symbols) for propane oxidation in the PSR-PFR configuration, compared with San Diego mechanism [40] (line), τ= 2.0 ms, T0=300 K, P= 1.0 atm,ɸ=1.0

To see how the above optimized P3 mechanism performs in other reactor models in the same operational conditions (i.e., inlet condition 1), we use the one dimensional non-homogenous PREMIX model of CHEMKIN. Figure 11 shows the spatial distribution of four species: O2, CO, CO2 and H2O along the reactor. Data (symbols) follow the reference San Diego mechanism [40] (lines) very closely, and the average relative errors are below 0.25% for all four species. 0.25 P=1.0 atm, T0=300 K, φ=1.0 Propane-air combustion

Lines: 247-reaction (SanDiego) Symbols: 3-rection (Present)

0.2 O2

Species mole fraction

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 28

H2O

0.15

0.1

CO2

0.05 CO

0 0

0.1

0.2

0.3

Distance along the PREMIXED(cm)

Figure 11 Species concentration along the PREMIX configuration based on optimized P3 mechanism (symbols) for propane oxidation in the PSR-PFR configuration, compared with San Diego mechanism (Line), T0=300 K, P= 1.0 atm,ɸ=1.0.

To study the validity of the multi-step global mechanism optimized for a special combustion condition for other residence times, we use the P3 mechanism produced for inlet condition 2 and 3 by a PSR-PFR reactor. The residence time for this procedure is 2 ms. The mechanism is then applied to a PSR or a PSR-PFR reactor, for different residence times (1-10 ms for PSR and 0.-0.4 ms for PSR-PSR), and results are compared with those of San

20 ACS Paragon Plus Environment

Page 21 of 28

Diego [40]. Figures 12-14 respectively show variations of flame temperature, CO and CO2 concentrations for different residence times. The left figures correspond to results from a PSR reactor, while the right figures are results for a PSR-PFR reactor. The reference data (lines) correspond to San Diego detailed mechanism [40]. Results show that the optimized P3 mechanism well follows the San Diego mechanism for all conditions, and the relative error for flame temperature and CO and CO2 concentrations are always well below 1%. Generally speaking, these results show that the DORP procedure is both efficient and reliable for most applications in propane oxidation.

1760

1780 Propane-air combustion T0=600 K, φ=0.5

Propane-air combustion T0=600 K,φ=0.5

1750

Flame temperature(K)

Flame temperature(K)

1760

1740

1730

1720 SanDiego mech, P=10 atm 3step DORP, P=10atm SanDiego mech, P=30 atm 3step DORP, P=30 atm

1710

1700

1

2

3

4

5

6

7

8

1740 SanDiego mech, P=10 atm 3 step, DORP, P=10 atm SanDiego mech, P=30 atm 3 step DORP, P=30 atm

1720

9

1700

10

0

0.1

Residence time in PSR(ms)

0.2 0.3 Residence time in PFR(ms)

0.4

(a) (b) Figure 12 Flame temperature versus residence times in (a) PSR and (b) PSR-PFR configurations based on the optimized P3 mechanism for propane oxidation, compared with the SanDiego mechanism [40], T0=600 K, P= 10, 30 atm, ɸ=0.5. 0.2

0.3

Propane-air combustion T0=600 K,φ=0.5

Propane-air combustion T0=600 K, φ=0.5

0.25

Mole fraction of CO(%)

0.15

Mole fraction of CO(%)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

0.2 SanDiego mech, P=10 atm 3step DORP, P=10 atm SanDiego mech, P= 30 atm 3step DORP, P= 30 atm

0.15

0.1

SanDiego mech, P=10 atm 3 step, DORP, P=10 atm SanDiego mech, P=30 atm 3 step DORP, P=30 atm

0.1

0.05 0.05

0

0 1

2

3

4

5

6

7

8

9

10

0

Residence time in PSR(ms)

0.1

0.2 0.3 Residence time in PFR(ms)

0.4

(a) (b) Figure 13 CO concentration versus residence times in (a) PSR and (b) PSR-PFR configurations based on the optimized P3 mechanism for propane oxidation, compared with the SanDiego mechanism [40], T0=600 K, P= 10, 30 atm,ɸ=0.5.

21 ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

0.062

0.061

Propane-air combustion T0=600 K,φ=0.5

Propane-air combustion T0= 600K,φ=0.5

Mole fraction of CO2(%)

0.06

Mole fraction of CO2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 28

0.059 SanDiego mech, P=10 atm 3step DORP, P= 10 atm SanDiego mech, p=30 atm 3step DORP, P= 30 atm

0.058

0.057

1

2

3

4

5

6

7

0.06

SanDiego mech, P=10 atm 3 step, DORP, P=10 atm SanDiego mech, P=30 atm 3 step DORP, P=30 atm

0.058

8

9

0.056

10

0

Residence time in PSR(ms)

0.1

0.2 0.3 Residence time in PFR(ms)

0. 4

(a) (b) Figure 14 CO2 concentration versus residence times in (a) PSR and (b) PSR-PFR configurations based on the optimized P3 mechanism for propane oxidation, compared with the SanDiego mechanism [40], T0=600 K, P= 10, 30 atm,ɸ=0.5.

5.4. Case 4: P3 mechanism for kerosene oxidation with air for ignition characteristics predictions

One of the most important properties of a fuel is its ignition characteristics. To use a reduced mechanism in design of an internal combustion engine, it need to be accurate enough in prediction of combustion properties during the ignition phase of flame formation. In this case, we study capabilities of the DORP procedure in optimization of a multi-step global mechanism to be well accurate for ignition simulation. For this case we optimize oxidation of kerosene, which is a heavy hydrocarbon, using the same very simple three-step of table 2, called K3 mechanism, which was also used in case 4 for oxidation of propane. We analyze how this optimized K3 mechanism may predict the ignition delay time of oxidation of kerosene. Kerosene, a heavy hydrocarbon fuel, is a mixture of alkanes, cycloalkanes, and aromatics. Its kinetics has been studied at different range of operating conditions with experimental tools, such as jet stirred reactors, shock tubes, and flat flame burners. In these experimental studies, many different mixtures of kerosene are used. Some of these studies are listed here: a mixture including 79% n-undecane, 11% 1,2,4-trimethylbenzene, 10% n-propyl cyclohexane at the atmospheric pressure and temperature range of 873-1033 K [42], mixture of 100%-n-decane at pressure range of 1.0-40.0 atm and temperature range of 750-1150 K [43], mixture of 100%-n-decane at temperature range of 550-1600 K at the atmospheric pressure [44], four mixtures including 100% n-decane, 74% n-decane, 26% n-propylbenzene, 74% n-decane, 26% n-propylcyclohexane, and 74%n-decane, 15%npropylbenzene, 11% n-propylcyclohexane under atmospheric pressure and temperature range of 900-1300 K [45], again these four mixtures under pressure of 1-40 atm [46] and many other studies such as [47-50].

22 ACS Paragon Plus Environment

Page 23 of 28

To produce K3 mechanism, we use the mixture considered by Wang et al. [2], i.e., 74%n-decane, 15%npropylbenze, 11% propyl cyclohexane. Details of implementation of the DORP procedure is very similar to case 4, and are not repeated here. This time the IGNITION reactor of CHEMKIN is used as the CRM engine. During this procedure, appropriate ranges for all eleven rate coefficients are selected by us, and using a cost function defined by equation (13) all eleven rate coefficients are optimized. Figure 15 shows results of the K3 mechanism (symbols) in oxidation of kerosene in an IGNITION reactor model, and compares these results with results of 382-step skeletal mechanism of Wang [2]. The inlet condition is: T0=1500 K, P= 1.0 atm, and ɸ=1.0. Temporal variation of O2, CO and CO2 species mole fractions predicted by both mechanisms are shown in figure 15. One observes that the very simple mechanism of K3 reproduces the temporal variation of major reactants and products with a high degree of accuracy. Again the relative errors are always much less than 1%. The ignition delay time may be defined as the time of maximum slope for CO mole fraction (time of maximum rate of production). This is about 195 µs for both mechanisms, and any difference between

them is hardly observed.

100 P=1.0 atm, T0=1500 K, φ=1.0 C10H20-air combustion O2

Species mole fraction

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

10

-1

CO

10

-2

CO2 Lines: 382-step (Wang 2012) Symbols: 3-step(Present)

10

-3

150

200

250

300

Time(µs)

Figure 15 Temporal variation of species concentration based on IGNITION reactor for oxidation of kerosene (C10H20) based on the optimized K3 mechanism, compared with Wang skeletal mechanism [2], T0=1500 K, P= 1.0 atm,ɸ=1.0.

6. CONCLUSIONS Access to simple chemical models is vital for three dimensional analysis of combustion chambers in all industrial applications. In this article a new approach, called DORP, is presented for fast and accurate optimization of multi-step global mechanisms for oxidation of hydrocarbon fuels. The differential evolution optimization procedure and chemical reactor models used in this procedure are described in detail. It is shown that modular application of cost functions and CRMs, makes application of the DORP procedure very easy, fast and effective. It is 23 ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 28

shown that any ideal CRM (or their combinations), and any simple or complex cost function may be easily used for this procedure. Therefore we claim that anyone may design global mechanisms just for application of his own interest. There is almost no restriction in selection of the cost function, or a combination of them. To evaluate how the procedure is applied and how accurate is the off-design performance of the resulted mechanism, five cases are studied. The first two cases study a five-step mechanism M5 produced to simulate oxidation of methane at different conditions. Here accuracy in prediction of concentration of pollutant species in the products of combustion and the flame temperature were design criteria. We observed that for a wide range of off-design equivalence ratios the relative errors were always below 1% for high pressure conditions, while it was still acceptable for low pressure simulations. In the second case the spatial accuracy in a PFR reactor was analyzed. Most remarkable observation is that the optimized M5 mechanism closely follows the reference mechanism of GRI3.0 in most off-design conditions considered here. In cases 3 and 4 the procedure is applied to other fuels, i.e., propane and kerosene. In case 3 it is shown that the P3 global mechanism, although very simple, still works very well in different off-design conditions, i.e., different residence times, and for reactor models different from its parent CRM, e.g., PREMIX reactor, and for different combustion criteria, e.g., spatial variation of flame temperature and pollutants concentration. Errors are surprisingly below 1% in most conditions. The reference mechanism for P3 is the detailed mechanism of San Diego. Finally case 4 is devoted to temporal assessment of the K3 mechanism, an optimized 3-step mechanism for simulation of kerosene oxidation. The reference mechanism is that of Wang, and the optimized K3 follows within a small margin the reference mechanism for most conditions studied here. One concludes that the DORP procedure is a simple, and fast procedure and is well suited for design of small multi-step global mechanisms, optimized for a range of inlet and operating conditions. The procedure is very flexible in usage of base optimization reactor models, and performance criteria. For an industrial combustor with different regions and combustion conditions, it is also easy to design different mechanisms, each appropriate for a zone of computation. Acknowledgment This work was supported by Sharif University of Technology. Interested researchers may ask to receive a copy of codes used here by sending a request to our email address.

REFERENCES 1. Shakeri, A.; Mazaheri, K. Using Sensitivity Analysis and Gradual Evaluation of Ignition Delay Error to Produce Accurate Low-Cost Skeletal Mechanisms for Oxidation of Hydrocarbon Fuels under High-Temperature Conditions. Energy. Fuels. 2017, 31, 11234. 24 ACS Paragon Plus Environment

Page 25 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Wang, Q. D.; Fang, Y. M.; Wang, F.; Li X.Y. Skeletal Mechanism Generation for High Temperature Oxidation of 2. Kerosene Surrogates. Combust. Flame. 2012, 159, 91. 3. Mallampalli, H.P.; Fletcher, T. H.; Chen, J.Y. Evaluation of CH4/NOx Reduced Mechanisms Used for Modeling Lean Premixed Turbulent Combustion of Natural Gas. J. Eng. Gas. Turbine. Power. 1998, 120, 703. 4. Lam, S. H.; Goussis, D. A. The CSP Method for Simplifying Kinetics. Int. J. Chem. Kinet. 1994, 26, 461. 5. Sung, C. J.; Law, C. K.; Chen, J. Y. Augmented Reduced Mechanisms for NO Emission in Methane Oxidation. Combust. Flame. 2001, 125, 906. 6. Kim, J. P.; Schnell, U.; Scheffknecht, G. Comparison of Different Global Reaction Mechanisms for MILD Combustion of Natural Gas. Combust. Sci. Technol. 2008, 180, 565. 7. Novosselov, I. V.; Malte, P. C. Development and Application of an Eight Step Global Mechanism for CFD and CRN Simulations of Lean Premixed Combustors. J. Eng. Gas. Turbine. Power. 2008, 130, 1. 8. Anderson, J.; Rasmussen, C. L.; Giselsson, T.; Glarborg, P. Global Combustion Mechanisms for Use in CFD Modeling Under Oxy-Fuel Conditions, Energy. Fuels. 2009, 23, 1379. 9. Franzelli, B.; Riber, E.; Sanjose, M.; Poinsot. T. A Two-Step Chemical Scheme for Kerosene–Air Premixed Fames. Combust. Flame. 2010, 157, 1364. 10. Westbrook, C. K.; Dryer, F. L. Simplified Reaction Mechanisms for the Oxidation of Hydrocarbon Fuels in Flames. Combust. Sci. Technol. 1981, 27, 31. 11. Jones, W. P.; Lindstedt, R. P. Global Reaction Schemes for Hydrocarbon Combustion. Combust. Flame. 1988, 73, 233. 12. Nicol, D. G.; Malte, P. C.; Hamer, A. J.; Roby, R. J.; Steele, R. C. Development of a Five-step Global Methane Oxidation-NO Formation Mechanism for Lean Premixed Gas Turbine Combustion. Eng. Gas. Turbine. Power. 1999, 121, 272. 13. Polifk, W.; Geng, W.; Dobbeling, K. Optimization of Rate Coefficients for Simplified Reaction Mechanisms with Genetic Algorithm. Combust. Flame. 1998, 113, 119. 14. Muller, U. C.; Peters, N.; Linan, A. Global Kinetics for n-heptane Ignition at High Pressures. Int. Symp. Combust. 1992, 24, 777. 15. Choi, J. Y. A Quasi Global Mechanism of Kerosene Combustion for Propulsion Applications., 47th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, San Diego, California, 2011. 16. Ruckert, F. U.; Sabel, T.; Schnelll, U.; Hein, K. R. G.; Risio, B. Comparison of Different Global Reaction Mechanisms for Coal Fired Utility Boilers. Prog. Comput. Fluid. Dy. 2003, 3,130. 17. Hurt, R. H.; Calo, J. M. Semi Global Intrinsic Kinetics for Char Combustion Modeling. Combust. Flame. 2001, 125, 1138. 18. Hautman, D. J.; Dryer, F. L.; Schug, K. P.; Glassman, I. A Multi-Step Overall Kinetic Mechanism for the Oxidation of Hydrocarbons. Combust. Sci. Technol. 1981, 25, 219. 19. Harris, S. D.; Elliott, L.; Ingham, D. B.; Pourkashanian, M.; Wilson, C. W. Optimization of Reaction Rate Parameters for Chemical Kinetic Modeling of Combustion Using Genetic Algorithms. Comput. Methods. App. Mech. Eng. 2000, 190, 1065. 20. Elliot, L.; Ingham, D. B.; Kyne, A. G.; Mera, N. S.; Pourkashanian, M.; Wilson, C. W. Multi Objective Genetic Algorithm Optimization for Calculating the Reaction Rate Coefficients for Hydrogen Combustion. Ind. Eng. Chem. Res. 2003, 42, 1215. 21. Elliot, L.; Ingham, D. B.; Kyne, A. G.; Mera, N. S.; Pourkashanian, M.; Wilson, C. W. Genetic Algorithms for Optimization of Chemical Kinetics Reaction Mechanisms. Prog. Energy. Combust. Sci. 2004, 30, 297. 22. Elliot, L.; Ingham, D. B.; Kyne, A. G.; Mera, N. S.; Pourkashanian, M.; Wilson, C. W. A Novel Approach to Mechanism Reduction Optimization for an Aviation Fuel/Air Reaction Mechanism Using Genetic Algorithm. J. Eng. Gas. Turbine. Power. 2004, 128, 255.

25 ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 28

23. Elliot, L.; Ingham, D. B.; Kyne, A. G.; Mera, N. S.; Pourkashanian, M.; Whittaker, S. Reaction Mechanism Reduction and Optimization for Modeling Aviation Fuel Oxidation Using Standard and Hybrid Genetic Algorithms. Comput. Chem. Eng. 2006, 30, 889. 24. Amorim, E. P. S.; Xavier, C. R.; Campos, R. S.; Santos, R. W. Comparison Between Genetic Algorithms and Differential Evolution for Solving the History Matching Problem. Proc. Comput. Sci. Appl. 2012, 12, 635. 25. Mazaheri, K.; Shakeri, A. Numerical Optimization of Laboratory Combustor Geometry for NO Suppression. Appl. Therm. Eng. 2016, 102, 1328. 26. Duterque, J.; Borghi, R.; Tichtinsky, H. Study of Quasi Global Schemes for Hydrocarbon Combustion. Combust. Sci. Technol. 1981, 26, 1. 27. Dupont, V.; Pourkashanian, M.; Williams, A. Modeling of Process Heater fired by Natural Gas. J. Inst. Energy. 1993, 66, 20. 28. Dryer, F. L.; Glassman, I. High Temperature Oxidation of CO and CH4. Int. Symp. Combust. Inst. 1973, 14, 987. 29. Kee, R. J.; Rupley, F. M.; Miller, J. A. CHEMKIN II: A FORTRAN chemical kinetics package for the analysis of gasphase chemical cinetics. Sandia Report 1989, SAND89-8009, Sandia National Laboratories. 30. Glarborg, P.; Kee, R. J.; Grcar, J. F.; Miller, J. A. A FORTRAN Program for Modeling Well-stirred Reactors. Sandia National Laboratories 1986, Livermore, CA, Report No. SAND86-8209. 31. Kee, R. J.; Grcar, J. F.; Smooke, M. D.; Miller, J. A. A FORTRAN program for modeling steady laminar onedimensional PREMIX flames 1985. Sandia Report, SAND85-8240, Sandia National Laboratories. 32. Turns, S. R. An introduction to combustion; Third Ed. Mc Graw hill Company: New York, 2012. 33. Stewart, P. H.; Larson, C. W.; Golden, D. M. Pressure and Temperature Dependence Reactions Proceeding Via a Bound Complex. 2. Application to 2CH3=>C2H5+H*. Combust. Flame. 1989, 75, 25. 34. Lindemann, F. A.; Arrhenius, S.; Langmuir, I.; Dhar, V; Perrin, V; Lewis, W. C. M. Discussion on the Radiation Theory of Chemical Action. Transactions of The Faraday Society. 1992, 17, 598. 35. Storn, R.; Price, K.V. J. Global Optimization. 1997, 11, 341. 36. Price, K. V.; Storn, R. M., Lampinen, J. A. Differential Evolution - A Practical Approach to Global Optimization. Natural Computing Series, Springer-Verlag, Berlin, Germany, 2005. 37. Mishra, S. Some New Test Functions for Global Optimization and Performance of Repulsive Particle Swarm Method, North-Eastern Hill University, Shillong (India), 23 August 2006. 38. Bowman, C.; Hanson, R.; Davidson, D.; Gardiner, W. J.; Lissianski, V.; Smith, G.; Golden, D.; Frenklach, M.; Goldenberrg, M. GRI 3.0 Detailed Mechanism, Accessed on 30 July 1999; http//www.me.berkeley.edu/gri3_mech/. 39. Karalus, M. F.; Fackler, K. B.; Novosselov, I. V.; Kramlich, J. C.; Malte, P. C. A Skeletal Mechanism for the Reactive Flow Simulation of Methane Combustion. Proceedings of The ASME Turbo Expo 2013: Tutbine Technical Conference and Exposition, San Antonio, Texas, U.S.A., June 3-7, 2013. 40. http://www-mae.ucsd.edu/~combustion/cermech/ 41. Bowman, C. T.; Hanson, R. K.; Gardiner, W. C; Lissianski, V.; Frenklach, M.; Goldenberg, M.; Smith, G. P.; Crosley, D. R.; Golden, D. M. GRI Mech 2.11, http://www.gri.org, 1995. 42. Dagaut, P.; Reuillon, M.; Boettner, J. C.; Cathonnet, M. Kerosene Combustion at Pressures Up to 40 atm: Experimental Study and Detailed Chemical Kinetic Modeling. Int. symp. Combust. Inst. 1994, 25, 919. 43. Battin-Leclerc, F.; Fournet, R.; Glaude, P. A.; Judenherc, B.; Warth, V.; Come, C. M.; Scacchi, G. Modeling of the Gas Phase Oxidation of n-Decane from 550 to 1600 K. Proc. Combust. Ins. 2002, 28, 1597. 44. Dagaut, P.; Ristori, A. The Combustion of Kerosene: Experimental Results and Kinetic Modeling Using 1-to 3Component Surrogate Model Fuels. Fuel. 2006, 85, 944. 45. Dagaut, P. Kinetics of Jet Fuel Combustion Over Extended Conditions: Experimental and Modeling. J. Eng. Gas. Turb. Power. 2007, 129, 394. 46. Cathonnet, M.; Gueret, C. B.; Chakir, A.; Dagaut, P.; Boettner, J. C.; Schultz, J. L. On the Use of Detailed Chemical Kinetics to Model Aeronautical Combustors Performances. Proceeding of the third European Propulsion Forum 1992, EPF91, ONERA 26 ACS Paragon Plus Environment

Page 27 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Paris, AAAF.

47. Patterson, P. M.; Kyne, A. G.; Pourkashanian, M.; Williams, A.; Wilson, C. W. Combustion of Kerosene in Counter-flow Diffusion Flame. J. Propul. Power. 2000, 16, 453.

48. Dagaut, P.; Ristori A. E. ; Bakali, A.; Cathonnet, M. Experimental and Kinetic Modeling Study of the Oxidation of n-propyl Benzene. Fuel. 2002, 81, 73.

49. Starik, A. M.; Titova, N. S.; Torokhov, S. A. Kinetics of Oxidation and Combustion of Complex Hydrocarbon Fuels: Aviation Kerosene. Combust. Explos. Shock waves. 2013, 49, 392.

50. Cathonnet, M.; Voisin, D.; Etsouli, A.; Sferdean, C.; Reuillon, M.; Boettner, J. C. Kerosene Combustion Modeling Using Detailed and Reduced Kinetic Mechanisms. Symposium applied vehicle technology panel on gas turbine engine combustion 1999; 14: 1-9. 51. Lee, D.; Park, J.; Jin, J. A Simulation for Prediction of Nitrogen Oxide Emission in Lean Premixed Combustor. J. Mech. Sci. Technol. 2011, 25, 1871. 52. Novosselov, I. V. Chemical Reactor Network Modeling of Combustion Systems, Ph.D. Dissertation, University of Washington, 2006. 53. Shakeri, A.; Mazaheri, K. Numerical Study of Pollutant Emissions in a Jet Stirred Reactor under Elevated Pressure Lean Premixed Conditions. Math. Probl. Eng. 2016, 2016, 1.

27 ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

DORP Procedure for Design of Small and Accurate Multi-step Global Chemical Mechanisms

Fuel Type Select Full or Skeletal Mechanism as Benchmark Define Operating Conditions (Point or range)

Chemical Reaction

Models (CRMs)

NO Concentration (ppmv)

500 Methane-air combustion T0=600 K, P= 1 atm, =2 ms

400 300 200 100 0

0.4

0.5

0.6 0.7 0.8 0.9 Equivalence ratio ()

1

500 Define Operating Conditions (Point or Range)

DE Optimization Code

Select a Multi-step Global Mechanism from Literature

Initialize Design Variable Define Design Variables of the Selected Multi-step Global Mechanism

NO Concentration (ppmv)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

Page 28 of 28

Methane-air combustion T0=600 K, P= 1 atm, =2 ms

400 300 200 100 0

ACS Paragon Plus Environment

Optimized 5 step Global Mech. GRI-2 GRI-3 Skeletal 5 step Mallampalli

0.4

0.5

0.6 0.7 0.8 0.9 Equivalence ratio ()

1