Development of a reduced chemical-kinetic combustion model for practical fuels

(1)

Development of a Reduced Chemical-Kinetic

Combustion Model for Practical Fuels

A thesis accepted by the

Faculty of Aerospace Engineering and Geodesy of the Universit¨

at Stuttgart

in partial fulfillment of the requirements for the degree of

Doctor of Engineering Sciences (Dr.-Ing.)

by

Anton Zizin

born in Voronezh, Russia

Committee chair:

Herr Prof. Dr. rer. nat. Uwe Riedel

Committee member: Herr Prof. Dr.-Ing. Andreas Kronenburg

Date of defence:

12.06.2013

Institute of Combustion Technology for Aerospace Engineering

University of Stuttgart

(2)

(3)

List of Figures

2.1 Distillation curves for 3 gasoline samples. . . 21

2.2 Distillation curves for kerosene samples. . . 22

2.3 Distillation curves for diesel samples. . . 22

2.4 Phase diagram for Jet-A and diesel fuels. . . 23

2.5 Schematic flowchart for the calculation procedure of the phase diagram and critical properties of a mixture. . . 32

2.6 Phase diagrams for an H2S-CO2-N2-CH4-C2H6-C3H8-nC4H10-isoC4H10-C5H12 mixture. Comparison of the experimental data of Boukouvalas et al. and calculation using the proposed algorithm and the REFPROP© program. . . . 33

2.7 Equilibrium pressure vs. H2S mole fraction in an H2S-C6H14-C15H32 mixture. Comparison of the experiment and calculation. k12=0.07, k13= k23= 0. . . 34

2.8 Calculated and experimental dependencies of critical pressure and temperature for n-C4H10-CO2 system. Mixing rule parameter k12 = 0.06. . . 34

2.9 Experimental and calculated distillation curves for 75/25 and 50/50 mole frac-tion mixtures of C10H22-C14H30 at P = 83kPa and P = 70 kPa. . . 35

2.10 Phase diagrams calculated using developed algorithm for fuel surrogates com-pared to experimental data for Jet-A kerosene. . . 37

2.11 Experimental distillation curves for gasoline and calculated boiling points curve for RD 387 and its surrogate fuel. The near to constant boiling temper-ature of the surrogate is caused by almost equal boiling tempertemper-atures of the surrogate components. . . 38

2.12 Calculated two-phase diagram for RD 387 and its reference fuel. Dew-point curve and bubble-point curve for surrogate are indistinguishable, because of very close properties of the surrogate components. . . 39

2.13 Experimental distillation curves for kerosene samples and calculated fuel boil-ing points curve for the proposed fuel surrogate. . . 39

2.14 Calculated phase diagram for the proposed fuel surrogate, Table 9, compared to experimental data for Jet-A. . . 40

3.1 Hierarchical structure of the mechanisms of hydrocarbons combustion. . . 46

3.2 Main paths of the paraffins combustion model. . . 48

5.1 Schematic flowchart for the calculation procedure of the thermodynamical properties of molecules with the Benson group additivity method. Low tem-peratures are 300 - 1000 K, high temtem-peratures are 1000-5000 K. . . 61

5.2 Comparison of the calculated normalized heat capacity of C8H18(2-methylheptane) with the experimental data. . . 62

5.3 Comparison of the calculated normalized heat capacity of C12H23OO, calcu-lated with DLR-Therm and calcucalcu-lated using NASA polynomials. . . 63

5.4 Comparison of the normalized heat capacity of the C12H25OO molecule, cal-culated with DLR-Therm and calcal-culated using NASA polynomials from the JetSurf mechanism. . . 63

5.5 Comparison of the normalized heat capacity of C12H23OOH, calculated with DLR-Therm and calculated using NASA polynomials from the JetSurf mech-anism. . . 64

(6)

5.6 Dependency of the normalized heat capacity on temperature for molecules from the C12H26 low temperature oxidation submodel. Calculation is made with DLR-Therm. . . 64 5.7 Dependency of the normalized heat capacity on temperature for the C16H34

molecule calculated with DLR-Therm and compared with other works. . . 65 5.8 Dependency of the normalized heat capacity on temperature for the C16H33

molecule calculated with DLR-Therm and compared with other works. . . 65 5.9 Dependency of calculated normalized heat capacity on temperature for C16H32

and C16H31 molecules calculated with DLR-Therm and compared with other works. . . 66 5.10 Dependency of the normalized heat capacity on temperature for molecules from

the C16H34 low temperature oxidation submodel calculated with DLR-Therm. 66 5.11 Comparison of experimental and calculated ignition delay times for an i–C8H18/

air mixture. Stoichiometric mixture. Pressure approximately 13.5 bar. . . 67 5.12 Comparison of experimental and calculated ignition delay times for an i–C8H18/

air mixture. Stoichiometric mixture. Pressure approximately 40 bar. . . 68 5.13 Comparison of experimental and calculated ignition delay times for an i–C8H18/

air mixtures. Rich mixture, φ=2. Pressure approximately 40 bar. . . 68 5.14 Comparison of experimental and calculated ignition delay times for an i–C8H18/

air mixture. Lean mixture, φ=0.5. Pressure approximately 40 bar. . . 69 5.15 Comparison of experimental and calculated ignition delay times for a

methyl-naphthalene/air mixture. Lean mixture, φ=0.5. Pressures approximately 11 and 36 bar. . . 69 5.16 Comparison of experimental and calculated ignition delay times for a

methyl-naphthalene/air mixture. Stoichiometric mixture. Pressures approximately 10 and 40 bar. . . 70 5.17 Comparison of experimental and calculated ignition delay times for a n-decane/

air mixture. Stoichiometric mixture. Pressures approximately 8.9, 11 and 13.5 bar. . . 72 5.18 Comparison of experimental and calculated ignition delay times for a n-decane/

air mixture. Stoichiometric mixture. Pressures approximately 40 and 50 bar. . 73 5.19 Comparison of experimental and calculated ignition delay times for n-decane/

air mixtures. Lean mixtures, φ=0.5 and φ=0.67. Pressure approximately 50 bar. . . 73 5.20 Comparison of experimental and calculated ignition delay times for n-decane/

air mixture. Lean mixture, φ=0.5. Pressures approximately 12 and 13 bar. . . 74 5.21 Comparison of experimental and calculated ignition delay times for n-decane/

air mixture. Lean mixture, φ=0.25. Pressures approximately 10 and 33 bar. . 74 5.22 Comparison of experimental and calculated ignition delay times for a n-decane/

air mixture. Rich mixture, φ=2. Pressures approximately 13 and 50 bar. . . . 75 5.23 Comparison of experimental and calculated ignition delay times for a 70%

α-methylnaphthalene/30% n-decane/air mixture. Stoichiometric mixture. Pres-sures approximately 12 and 50 bar. . . 75 5.24 Comparison of experimental and calculated ignition delay times for a 30%

α-methylnaphthalene/70% n-decane/air mixture. Stoichiometric mixture. Pres-sures approximately 19 and 50 bar. . . 76

(7)

5.25 Comparison of experimental and calculated ignition delay times for a n-decane/air mixture. Stoichiometric mixture. Pressures approximately 14 and 41 bar. . . 78 5.26 Comparison of experimental and calculated ignition delay times for a n-decane/

air mixture. Lean mixture, φ=0.5. Pressures approximately 15 and 42 bar. . . 78 5.27 Comparison of experimental and calculated ignition delay times for a n-dodecane/

air stoichiometric mixture. Pressure approximately 23 bar. . . 79 5.28 Comparison of experimental and calculated ignition delay times for a lean

n-dodecane/air mixture. φ=0.5. Pressure approximately 22 bar. . . 79 5.29 Comparison of experimental and calculated ignition delay times for a lean

n-dodecane/O2 mixture. φ=0.5. Dilution with argon 96-97%. Pressure ap-proximately 16 bar. . . 80 5.30 Comparison of experimental ignition delay times for n-dodecane/air and

calcu-lated ignition delay time for hexadecane/air mixtures. Lean mixtures, φ=0.5. Pressure approximately 22 bar. . . 81 5.31 Comparison of experimental ignition delay times for n-dodecane/air and

cal-culated ignition delay time for hexadecane/air mixtures. Stoichiometric mix-tures. Pressure approximately 23 bar. . . 82 5.32 Comparison of experimental ignition delay times for n-tetradecane/air and

calculated ignition delay time for hexadecane/air mixtures. Lean mixtures, φ=0.5. Pressures approximately 12 bar and 38 bar. . . 82 5.33 Comparison of experimental ignition delay times for n-tetradecane/air and

calculated ignition delay time for hexadecane/air mixtures. Stoichiometric mixtures. Pressures approximately 15 and 40 bar. . . 83 5.34 Comparison of experimental ignition delay times and calculated ignition delay

times for a Jet-A/air mixture. Stoichiometric mixture. Pressures approxi-mately 25 and 50 bar. . . 84 5.35 Comparison of experimental ignition delay times and calculated ignition delay

times for a Jet-A/O2/N2 mixture. Stoichiometric mixture. Pressure approxi-mately 18 bar. . . 84 5.36 Comparison of experimental ignition delay times and calculated ignition delay

times for a Jet-A/air mixture. Lean mixture, φ=0.5. Pressure approximately 20 bar. . . 85 5.37 Comparison of experimental ignition delay times and calculated ignition delay

times for a Jet-A/air mixture. Stoichiometric mixture. Pressures approxi-mately 22 and 25 bar. . . 85 5.38 Comparison of experimental ignition delay times and calculated ignition delay

times for a Jet-A/air mixture. Stoichiometric mixture. Pressure approximately 9 bar. . . 87 5.40 Comparison of experimental ignition delay times and calculated ignition delay

times for a Jet-A/air mixture. Rich mixture, φ=2. Pressure approximately 9 bar. . . 87 5.41 Rate of production analysis of the paraffins combustion mechanism. Initial

temperature 1192 K. Pressure 27 bar. Calculated ignition delay 93 µs. Analysis shown at t=50 µs. . . 88

(8)

5.42 Rate of production analysis of the methylnaphthalene combustion mechanism. Initial temperature 1192 K. Pressure 27 bar. Calculated ignition delay time

93 µs. Analysis shown at t=50 µs. . . 89

5.43 Rate of production analysis of the propyl-benzene combustion mechanism. Initial temperature 1192 K. Pressure 27 bar. Calculated ignition delay 93 µs. Analysis shown at t=50 µs. . . 90

5.44 Sensitivity analysis for the H atom. T=1531.8 K, p=9.251 bar, φ=1. Calcu-lated ignition delay 4.5 µs. Analysis shown at t=4.4 µs. . . 91

5.45 Sensitivity analysis for C16H34. T=1531.8 K, p=9.251 bar, φ=1.Calculated ignition delay 4.5 µs. Analysis shown at t=4.4 µs. . . 92

5.46 Sensitivity analysis for the H atom. T=880 K, p=25 bar, φ=1. Calculated ignition delay 2350 µs. Analysis shown at t=2300 µs. . . 92

5.47 Sensitivity analysis for C16H34. T=880 K, p=25 bar, φ=1. Calculated ignition delay 2350 µs. Analysis shown at t=2300 µs. . . 93

6.1 Mechanism reduction strategy. . . 94

6.2 The flow chart of the RedMaster-R program. . . 95

6.3 The flow chart of the RedMaster-R program. . . 96

6.4 Comparison of ignition delay times calculated using the full and set of reduced mechanisms 1.1-1.5 for a Jet-A/air mixture. Stoichiometric mixture. Pressure approximately 9 bar. . . 99

6.11 Comparison of ignition delay times calculated using the full and set of reduced mechanisms 4.2.1-4.2.5 for a Jet-A/air mixture. Stoichiometric mixture. Pres-sure approximately 25 bar. . . 103

(9)

6.14 Comparison of ignition delay times calculated using the full and set of reduced mechanisms 4.3.1-4.3.5 for a Jet-A/air mixture. Stoichiometric mixture. Pres-sure approximately 9 bar. . . 104 6.15 Comparison of experimental ignition delay times and calculated ignition delay

times for a Jet-A/O2/N2 mixture. Stoichiometric mixture. Pressure approxi-mately 18 bar. . . 105 6.16 Comparison of experimental ignition delay times and calculated ignition delay

times for Jet-A/air mixture. Stoichiometric mixture. Pressure approximately 25 bar. . . 105 6.17 Comparison of experimental ignition delay times and calculated ignition delay

times for a Jet-A/air mixture. Rich mixture, φ=2. Pressure approximately 9 bar. . . 106

(10)

List of Tables

1 Hydrocarbons used in the surrogate fuel models. . . 19

2 Thermodynamical properties of hydrocarbons used in surrogate fuel models. . 20

3 Types of fuels. . . 21

4 Comparison of main types of equations of state. . . 26

5 Equilibrium concentrations in the propane-cyclohexane mixture; experimental data are compared with calculations using the presented algorithm and the NIST REFPROP© package. . . 35

6 Composition of 7 hydrocarbons mixture, experimental and calculated critical temperature and pressure. . . 36

7 Composition of 5 hydrocarbons mixture, experimental and calculated critical temperature and pressure. . . 36

8 Compositions of tested mixtures. . . 38

9 Composition of the proposed blend for kerosene combustion and its basic prop-erties. . . 40

10 List of molecules for which thermodynamical properties are calculated. . . 62

11 C10H22 submechanism. Units are mole, s, cm3, Kelvin. . . 70

14 Aromatic compounds . . . 86

15 Reduced mechanisms obtained with RedMaster-S. . . 98

(11)

List of symbols

———————————————– Roman symbols ———————————————– A preexponential factor in the arrhenius function

a coefficients in the NASA polynomials

A Helmholtz free energy

b coefficient in the cubic equation of state ¯

B coefficient in the virial equation of state

C cubic term

¯

C coefficient in the virial equation of state

[C] concentration

¯

D coefficient in the virial equation of state Ea activation energy

f production rate

F coefficient in the Troe approximation

G Gibbs energy

H enthalpy

k reaction rate constant

k parameter in the mixing rules

l parameter in the mixing rules

M number of reactions

n number of moles

nT total number of moles in the system ¯

n power of T in the reaction rate

N number of substances or species in the system NT number of molecules in the system

P pressure pr reduced presure Q amount of heat Q stoichiometric matrix R gas constant S entropy T temperature U internal energy v molar volume V volume W work

x species volume fraction (in liquid phase) y species volume fraction (in gas phase)

Z compresibility

————————————————– Greek symbols ————————————————– α parameter of the cubic equation of state

γ parameter of the cubic equation of state ε parameter of the cubic equation of state Λ matrix which defines the limit of stability

µ chemical potential

(12)

ν stoichiometric coeficient

νij stoichiometric coefficient of the i-th species in the j-th reaction ξ parameter of the cubic equation of state

φ stoichiometric ratio ϕ fugacity coefficient Ψ arbitrary parameter ————————————————— Subscripts ————————————————— crit critical f formation i number of substance mix mixture ————————————————– Superscripts ————————————————— i number of substance L liquid phase r residual

V vapor or gas phase

° reference state

————————————————– Abbreviations ————————————————–

BWR Benedict-Webb-Rubin equation of state

CFD Computational Fluid Dynamics

EoS Equation of State

LES Large Eddy Simulations

NASA National Aeronautics and Space Administration NIST National Institute of Standards and Technology

PCA Principal Components Analysis

QSS Quasi-Steady-State

RANS Reynolds-Averaged-Navier-Stokes

SRK Soave-Redlich-Kwong EoS

TBP Temperature Boiling Point

(13)

Kurzfassung

Diese Arbeit befasst sich mit der Modellierung von Ersatzbrennstoffen sowie deren physikalischen und chemischen Eigenschaften. Aufgrund der großen Zahl an unterschiedlichen Komponenten in herkömmlichen Kraftstoffen ist es nicht möglich, die Eigenschaften des Kraftstoffes mit Hilfe von Gleichungen, die für alle vorhandenen Spezies aufgestellt werden, zu berechnen. Noch schwieriger ist es, für einen solchen Kraftstoff einen kinetischen Mechanismus zu entwickeln. Stattdessen wird ein Gemisch weniger, ausgewählter Kohlenwasserstoffe verwendet, das den tatsächlichen Kraftstoff bei der Modellierung ersetzt. Ein solches Gemisch wird als Surrogat bezeichnet. Das Ziel dieser Arbeit ist es, eine Methode zu entwickeln, die einen reduzierten Verbrennungsmechanismus für ein Surrogat für Jet-A aufbaut, angefangen von der Entwicklung des Modells für den Treibstoff bis zum Algorithmus für die Reduktion des kinetischen Mechanismus.

Vor der Entwicklung des numerischen Codes wird ein kurzer ¨Uberblick ¨uber vorhandene Zustandsgleichungen gegeben und deren Vor- und Nachteile sowie Genauigkeit besprochen.

Die entwickelten numerischen Werkzeuge erlauben eine schnelle Berechnung des

Phasendiagramms des Kraftstoffs, der Destillationskurve und des kritischen Punktes. Das physikalische Modell und die entwickelten Werkzeuge wurden anhand mehrerer Experimente f¨ur verschiedene Mischungen von Kohlenwasserstoffen erfolgreich validiert. Das erhaltene Surrogat-Modell wird durch die Optimierung der genannten Eigenschaften sowie der Verbrennungsenthalpie, der Bildungsenthalpie, des molekularen Gewichts und des Ruß-Index erhalten. Das verwendete Surrogat f¨ur Jet-A besteht aus 11% Propylcyclohexan, 14% iso-Oktan, 22% n-Dodecan, 28% 1-Methylnaphthalin und 25% n-Hexadecan.

Für eine Vielzahl großer Kohlenwasserstoffe sind bereits kinetische Modelle entwickelt worden, unter anderem für n-Decan, n-Dodecan und n-Hexadecan. Diese Modelle sind für ihren jeweiligen Geltungsbereich durch zahlreiche Stoßrohrexperimente validiert. Die Mechanismen für n-Dodecan und n-Hexadekan werden zu den Mechanismen der anderen Komponenten des Surrogats, die am Institut für Verbrennungstechnik des Deutschen Zentrums für Luft- und Raumfahrt entwickelt wurden, hinzugefügt. Der so entstandene

Gesamtmechanismus besteht aus 183 Spezies und 1239 Reaktionen. F¨ur diesen

Gesamtmechanismus werden numerische Simulationen durchgeführt und die Ergebnisse mit Stoßrohrexperimenten für Jet-A Kerosin verglichen. Von besonderem Interesse sind der Niedertemperaturbereich und der Bereich des negativen Temperaturkoeffizienten. Deshalb werden Sensitivitätsanalysen und Analysen der Produktionsraten für jede Surrogat-Komponente durchgeführt.

Auf Grund des Mangels an thermodynamischen Daten f¨ur Spezies im

Niedertemperaturbereich wird ein Code entwickelt, der die thermodynamischen Eigenschaften unter Verwendung der Gruppenadditivitätsmethode von Benson berechnet. Für die Extrapolation der Wärmekapazität zu höheren Temperaturen wurden, aufgrund ihrer guten Präzision und Stabilität, Willhoit- Polynome verwendet. Die Validierung der mit dem neuen Tool entwickelten Werte erfolgt für eine Reihe von Kohlenwasserstoffen durch den Vergleich mit Werten aus Datenbanken oder Berechnungen aus anderen Publikationen.

Für den Reduktionsalgorithmus werden mehrere verbreitete Reduktionsmethoden beschrieben und verglichen. Anhand ausgewählter Methoden wird der kinetische Mechanismus analysiert und reduziert. Schließlich werden für die Reduktion des gesamten Mechanismus zwei zusätzliche Programme entwickelt. Anhand der KINALC Programme werden zwei Analysen zu mehreren Zeitpunkten für ausgewählte numerische Experimente

(14)

durchgeführt: eine Analyse der Zusammenhänge zwischen den einzelnen Spezies und eine Analyse der Wichtigkeit der Reaktionen. Durch einen Vergleich der erstellten Listen werden automatisch die unwichtigen Spezies und Reaktionen erkannt und aus dem Mechanismus entfernt. Nach dieser Methode lassen sich hervorragende Ergebnisse, auch bei hohen Reduktionsgraden, erhalten. Allerdings existiert eine Untergrenze für die Anzahl der Spezies, unter welcher eine weitere Reduktion zu sehr schnell wachsenden Fehlern führt. Bei diesen Reduktionsgraden tendiert der Algorithmus dazu, einige Spezies zu entfernen, die für die Qualität des Mechanismus von grundlegender Bedeutung sind. Dies begrenzt die minimale Größe des Mechanismus auf 70-80 Spezies. Der Einfluss unterschiedlicher Reduktionsparameter auf die Qualität des resultierenden reduzierten Mechanismus wird untersucht und dargestellt.

(15)

Abstract

This work deals with modeling of the surrogate fuels including their physical as well as chemical properties. The huge number of components of real fuels makes it impossible to model its properties with equations for all of the constituents. Even less feasible is to create a kinetic mechanism for such a fuel. As an alternative a mixture of a few hydrocarbons is used to replace the real fuel. Such mixture is called a surrogate fuel. The goal of the present work is to construct a complete method for the creation of a reduced combustion model for a Jet-A surrogate fuel, starting from the development of the fuel model itself to the algorithms for a kinetic mechanism reduction.

Prior to the development of the numerical code a short overview of available equations of state is given and a comparison of their advantages, disadvantages and precision is made. The developed numerical tools allow fast calculation of the fuels phase diagram, distillation curve and critical point. The physical model and developed tools are successfully validated over a number of experiments for different mixtures of hydrocarbons. The elaborated sur-rogate model is obtained by optimization of the above mentioned characteristics as well as combustion enthalpy, formation enthalpy, molar weight and sooting tendency index. The Jet-A surrogate model consists of 11% propylcyclohexane, 14% iso-octane, 22% n-dodecane, 28% 1-methylnaphtalene and 25% n-hexadecane.

Kinetic models are developed for the large hydrocarbons: decane, dodecane and n-hexadecane. These models are thoroughly tested over a wide range of shock-tube experiments for each of these pure components as well as mixtures of n-decane and α-methylnaphthalene. The developed sub-mechanisms are added to the sub-mechanisms of other surrogate fuel components developed at the Institute of Combustion Technology at the German Aerospace Center (DLR) producing the full mechanism with 183 species and 1239 reactions. Numer-ical simulations made with the full mechanism for the surrogate blend are compared with shock-tube experiments for Jet-A kerosene. Of special interest is the low-temperature region and negative temperature coefficient. A rate of production analysis for each surrogate fuel component and sensitivity analysis will be given.

Due to the lack of the thermodynamical data for the low-temperature species another computational tool is developed. It calculates thermodynamical properties using the Benson group additivity method. For extrapolation of the heat capacity to higher temperatures Willhoit polynomials are used due to good precision and stability. The verification is made by comparison of the results for a range of hydrocarbon molecules with database values and calculations in other works.

For the reduction algorithm several common reduction methods are described and com-pared. Using some of these methods analysis and reduction of the full kinetic mechanism is made. Finally, for the reduction of the full mechanism two additional programs are devel-oped. Using the KINALC package they accomplish two analyses in several time points for chosen numerical experiments: importance of reactions and species. Comparing the obtained lists they automatically choose unimportant species and reactions and remove them from the mechanism. Excellent results are obtained even at high degrees of reduction. Nevertheless, there exists a minimal number of species, after which further reduction causes very fast grow-ing errors. At the same time the algorithm tends to remove some very basic species at high degrees of reduction. This limits the minimal size of the mechanism to 70-80 species. The influence of different reduction parameters on the performance of the resulting reduced mechanism is investigated and will be given.

(16)

1 Introduction

Modern energy infrastructure all over the world relies on combustion for energy conversion. Hydrocarbons are especially important for transportation and absolutely irreplaceable in air-borne transport due to their high energy density. For example, the heat of combustion of kerosene is 43M J

kg , whereas the specific energy of lithium-ion batteries is only 0.36-0.90 M J

kg . Hydrogen has an even higher specific energy 121M J_kg , but using it in practical applications entails many problems due to its high diffusivity and flammability. Moreover it has very low volumetric density and has to be highly compressed or liquefied, which causes technical problems and additional expenses. Renewable energy sources show capability for supple-menting and gradually replacing fossil fuels, but hydrocarbons will remain the main energy contributor for at least the next few decades.

Increasing environmental awareness and rising fuel costs require improvements in effi-ciency of energy infrastructure and reduction in pollution. Computational Fluid Dynamics (CFD) is a powerful tool that can help us achieve these goals through optimization of com-bustion processes. Recently several mathematical models and algorithms are implemented to ease large-scale simulations. Large Eddy Simulation (LES) and Reynolds-averaged Navier-Stokes equations (RANS) are widely used approximations of turbulent flow. Flamelet model and CFD look-up tables are used to reduce the computational cost of chemical kinetics cal-culations.

Despite these simplifications and significant increase in available computational power over the last several decades, exact modeling of combustion with all physical and chemical processes taken into consideration is still far beyond reach. This situation will not change for at least several decades. Computational costs (time and memory required) rise very fast with the growth of the number of variables (species) in the problem. Depending on the numerical method used in computation the cost of modeling is between the quadratic and cubic power of the number of variables in the problem. Currently for smaller molecules CFD simulations can be done with detailed or semi-detailed chemistry. For larger fuel molecules the only option are global mechanisms, which consist of a few global reactions only. Some authors also include reduced pollutants chemistry. Such calculations provide general information about flow, turbulence, and temperature. However, they do not give precise information about the underlying chemistry and pollutants, NOx and soot formation.

Regular fuels consist of hundreds of hydrocarbons each generating a large number of minor hydrocarbons during combustion. To avoid generation of bulky and unusable mechanisms one can use surrogate fuel models, which would consist of only a few components but still would be able to reproduce real fuel’s physical and chemical properties. Thorough selection of a few such components would be a first step in the creation of the reduced kinetic mechanism. Until now most of suggested surrogate fuels are proposed based on only a few physical or chemical properties: either spray formation, evaporation characteristics or combustion properties, such as ignition delay or a flame speed. This leads to a strong limitation of the applicability of corresponding surrogate.

The reduction of kinetic mechanisms is another problem for CFD simulations. The ac-cepted approach to reduction is to investigate the whole mechanism and then remove reactions and species, which are unimportant for the range of parameters considered (pressure, tem-perature, and equivalence ratio). This approach provides controllable errors in ignition delay time, flame speed or concentration profiles. Therefore, methods of mechanism analysis and reduction are tightly coupled. Many procedures are suggested since the first works of Semenov

(17)

and Zeldovich: quasi steady state assumption (QSSA), partial equilibrium assumption, sensi-tivity, atom flux analysis, principal component analysis (PCA), uncertainty analysis, intrinsic low dimensional manifold (ILDM), lumping (methods based on the merging of species con-centrations into groups), directed relation graph, and computational singular perturbation (CSP). Some of these methods allow one to reduce a detailed mechanism to a skeletal one, which usually is still too big for direct application in CFD. Others use complex expressions for reaction rates of these mechanisms, which can be prohibitive for some numerical methods used in commercial CFD codes.

Conversely, global mechanisms consist of a few reactions only. Such mechanisms cannot describe auto ignition and are used only in flame modeling. One way to build a global mechanism is to fold a number of elementary reactions into several global reactions. Another way is to choose the most important reactions in different flame zones and remove the rest.

The third approach is to start with 1-reaction mechanism and to build the mechanism bottom-up. Knowing what species concentration or physical properties need to be repro-duced, it is possible to construct several global reactions and then fit the reaction rates. The shortest possible mechanism is the following global reaction of paraffin combustion:

CxH2x +2+ (3x+1₂ )O2 → xCO2+ (x + 1)H2O

This reaction gives reasonable flame speeds for a considerable equivalence ratio range, but over-predicts the temperature. Therefore another reaction is added [1]:

CxHy+ (x₂ +₂y)O2 ↔ xCO + y₂H2O CO + 1₂O2 ↔ CO2

For an estimation of NOx production a third reaction is inserted [2]: CxHy+ (x₂ +₂y)O2 ↔ xCO + y₂H2O

CO + 1₂O2 ↔ CO2 N2+ O2 ↔ 2NO

This method works only for small mechanisms. It is impossible to produce such mech-anisms automatically and fit reaction rate constants for larger number of global reactions. This leads to the necessity of developing an intermediate reduction step between the skeletal mechanism and the global one.

The goal of the present work is to construct a complete method of creation of a reduced mechanism for a specific surrogate fuel, from development of the model fuel itself to the algorithms for reduction and elaboration of the global combustion mechanism. Therefore the central parts of the work are: (1) development of the surrogate fuel model, (2) development of the kinetic mechanism and calculation of the thermodynamical properties of the new molecules, and (3) development of reduction methods and elaboration of skeletal kinetic mechanisms for the proposed fuel model.

(18)

2 Surrogate Fuels

2.1 Surrogate Fuels

Common fuels distilled from crude oil, such as kerosene, petrol, diesel, as well as some modern bio-fuels consist of hundreds or thousands of individual molecular components. Each fuel must meet its certification standard but its chemical composition can vary significantly depending on the source, processing methods, and even the time of the year.

For that reason instead of constant adjustment of the fuel compound for the computer modeling it is sensible to use a model fuel or surrogate fuel, which consists of only a few different hydrocarbons but gives a good representation of the performance and emissions behavior of the real fuel. This makes it possible to create a chemical model, which can describe physical and chemical properties of the real fuel with sufficient precision and, at the same time, is simple enough for use in combustion and evaporation models. Some parts of Chapter 2 are published in [3].

2.2 Reference Fuel Evaporation Characteristics

The hydrocarbons in practical fuels are primarily paraffins, naphthenes (cycloparaffins) and aromatics. Table 1 summarizes the hydrocarbons used in different chemical models for gaso-line, kerosene and diesel [4–13]. These species have been used successfully to represent hun-dreds of hydrocarbons in surrogate mixtures in simulations of fuel oxidation. Let us consider the following two questions: How do these species represent the evaporation characteristics? What kind of criteria do we need to evaluate these characteristics?

First, the requirements for the fuel volatility are traditionally described by the distilla-tion curve (temperature boiling point, TBP curve - boiling temperature vs. the volume of evaporated fraction during distillation), which gives a distribution of species by boiling tem-peratures. Very often instead of the distillation curve a fuel is characterized by the percentage volume evaporated for a few boiling temperatures only. This defines maximal boiling temper-ature. The most important distillation characteristics are the temperatures at the top end of the range, as these provide a measure of the proportion of heavier components, which corre-late particularly closely with the aromatics content of the fuel. Heavier components in a fuel have more potential for incomplete vaporization and combustion, which results in increased production of smoke or soot. Predefinition of the lower end values of the distillation curve, which correspond to higher temperatures, reduces the proportion of the heavy components and produces a cleaner exhaust. It can also reduce the density and the viscosity of the fuel as these properties are closely connected. If one wants to describe the distillation curve for a mixture, one has to select a set of hydrocarbons with the boiling temperatures lying between the lowest and highest temperatures of the distillation curve. The boiling temperatures of individual hydrocarbons used in surrogate fuel models are provided in the Table 2.

Table 3 sums up the common types and properties of fuels, which will be used in the further investigation.

In most cases it is possible to generate a variety of different sets that will describe a given boiling range. The distillation curve for an individual fuel is usually obtained from experiments [15–17] and this data can be used for validation of the surrogate blend content [13, 17].

(19)

Table 1: Hydrocarbons used in the surrogate fuel models.

Type of fuel Primary reference fuel components

Gasoline [4, 11] toluene n-pentane iso-octane iso-pentane n-heptane cyclopentane ethanol cyclohexane iso-butane m-xylene benzene iso-octane toluene cyclooctane methylcyclohexane propylbenzene n-heptane 1,2,4- trimethylbenzene ethyl n-propylcyclohexane benzene Jet [4, 6–9, 11–13] n-nonane butylbenzene tetralin tetramethylbenzene decalin n-decane methylnaphthalene n-dodecane n-undecane n-tetradecane n-octylcyclopentane n-hexadecane n-pentadecane iso-octane heptamethylnonane decalin n-hexadecane tetraline naphthalene Diesel [14]

ethyl-/ propyl-/ butylcyclohexane n-decane

ethyl-/ propyl-/ butylbenzene n-dodecane

toluene tetradecane

methylcyclohexane n-decylbenzene

n-heptane iso-cetane

grades have different distillation curves. It should also be mentioned that there is not enough data describing distillation curves of various fuels. This can present difficulties in evaluating experimental errors in available measurements or finding data related to the selected fuel. The phase behavior (vapor-liquid equilibrium and accordingly surface tension and viscosity) over the full ranges of temperatures and pressure values is characterized by two-phase diagrams, Fig. 2.4. The situation with experimental two-phase diagrams is even worse than that with distillation curves – there are no other measurements of liquid-gas equilibrium for real fuels besides the data reported in [11] and for gasoline there is no data at all. The phase equilibrium defines the chemical composition of a surrogate mixture: blends and real fuel must have identical boiling range/phase behavior, critical temperature and pressure.

For mixtures, these properties are not linearly dependent on the percentage of the com-ponents and are not easy to fit. They are calculated from the equation of state, i.e. they are linked directly to the thermodynamic state functions. For the effective development of the optimal surrogate fuel model, one has to approximate the thermodynamic functions of the multi-component real fuel by those of the simplified surrogate. In what follows, we identify tools that can be used to calculate the thermodynamic equilibrium, i.e. two-phase diagrams

(20)

Table 2: Thermodynamical properties of hydrocarbons used in surrogate fuel models [18].

Formula Name Tboil Tcrit Pcrit

n- and iso- paraffins CH(CH3)3 C5H12 C5H12 C6H14 C7H16 C7H16 C8H18 C8H18 C9H20 C10H22 C11H24 C12H26 C13H28 C14H30 C15H32 C16H34 C16H34 iso-butane n-pentane iso-pentane hexane heptane 2-methylhexane n-octane iso-octane nonane n-decane undecane dodecane tridecane tetradecane pentadecane hexadecane heptamethylnonane, HMN 261.34 309.22 300.99 341.88 371.57 363.18 398.82 372.39 423.97 447.3 469.08 489.48 508.63 526.76 543.83 559.98 520.25 407.85 469.7 460.39 507.6 540.2 530.1 568.70 543.9 594.6 617.7 639 658 675 693 708 723 693 36.4 33.7 33.81 30.25 27.4 27.3 24.9 25.7 22.9 21.1 19.8 18.2 16.8 15.7 14.8 14 15.7 Cycloparaffins C6H12 C7H14 C8H10 C8H16 C8H16 C9H18 C10H18 cyclohexane methylcyclohexane m-xylene cyclooctane ethylcyclohexane propylcyclohexane trans-decalin 353.93 374.09 412.34 424.31 405.0 429.0 460.42 553.5 572.19 617 647.2 609 625.6 687 40.73 34.71 35.41 35.7 27.7 27.1 46 Aromatics C6H6 C7H8 C8H10 C9H12 C9H12 C10H12 C11H10 benzene toluene ethylbenzene propylbenzene 1,2,4-trimethylbenzene tetralin 1-methylnaphtalene 353.24 383.79 409.36 432.35 442.5 480.75 517.84 562.05 591.75 617.15 638.35 649.1 720 772 48.95 41.08 36.09 32 32.32 36.5 36

and critical points, for multi-component mixtures. Using these tools we calculate two–phase diagrams for a few surrogates from literature and optimized chemical contents of a blend, which satisfies the thermodynamic properties of practical fuels.

In calculation of vapor-liquid equilibrium for mixtures one faces two main problems: first, choosing an equation of state (EoS), which can provide reliable results, is simple, and would not require additional calculations to determine the equation parameters. The second prob-lem is choosing the mixing rules to evaluate EoS parameters for multi-component mixtures.

(21)

Table 3: Types of fuels.

Fuel Description

Jet-A

Aviation kerosene based fuel designed for use in aircraft powered by gas-turbine engines. Generally consists of hydrocarbons with

8-16 carbon atoms.

Jet-A1 Same as Jet-A but has lower freezing point (-47 °C vs. -40 °C). Jet-B Aviation kerosene with enhanced cold-weather performance due to

higher naphtha content (hydrocarbons with 5-15 carbon atoms). JP-8 (Jet Propellant-8) Military jet fuel, also used as a replacement of

diesel fuels.

S-8 Synthetic aliphatic jet fuel, which is free of aromatic hydrocarbons, an alternative to JP-8.

AI-91 Gasoline with an octane number 91.

RD 387 Gasoline with an octane number 87 .

For non-ideal cases (big difference in molecular size, high polarity) simple mixing rules are not applicable. -10 0 10 20 30 40 50 60 70 80 90 100 110 280 300 320 340 360 380 400 420 440 460 480 T b o il , K

Distillate volume fraction, %

Gasoline, experiment, de Oliveira et al. AI 91 gasoline, experiment, Smith et al. Gasoline, experiment, Anderson et al.

Figure 2.1: Distillation curves for 3 gasoline samples: de Oliveira et al. [15], Andersen et al. [19] and distillation curves of the gasoline with antiknock index (AI) 91 from Smith et al. [16].

(22)

0 20 40 60 80 100 420 440 460 480 500 520 540 560

580 _{Jet-A, experiment, Smith et al.}

Jet-A, experiment, Violi et al. Jet-A #1, experiment, Rachner et al. Jet-A #6, experiment, Rachner et al. Jet-A #7, experiment, Rachner et al. Jet-A #8, experiment, Rachner et al. Jet-A #11, experiment, Rachner et al.

T b

o

il

,

K

Distillate Volume Fraction, %

Figure 2.2: Distillation curves for kerosene samples: experimental distillation curve for Jet-A from Smith et al. [20], boiling-point curves for Jet-A fuel from Violi et al. [13] and 5 different Jet-A kerosenes from Rachner [21].

0 10 20 30 40 50 60 70 80 90 100 380 400 420 440 460 480 500 520 540 560 580 600 620 640 660 680 700

Diesel (Atas), experiment, Demirbas Diesel (Iran), experiment, Demirbas Diesel, experiment, Bugorkova et al.

T b

o

il

,

K

Distillate Volume Fraction, %

Figure 2.3: Distillation curves for diesel samples. Experimental distillation curve from Demir-bas [22] and Bugorkova et al. [23].

(23)

600 620 640 660 680 10 15 20 25 P , a tm T, K Vapor Luquid

Figure 2.4: Phase diagram for Jet-A and diesel fuels from [11]; symbol * used for critical point; symbols ◦ used for bubble points, and symbols × used for dew points.

2.3 Thermodynamical Properties of Gases and Liquids

Physical properties of all substances depend directly on the nature of the molecules they consists of. Theoretically one can calculate all physical properties of fluids and gases from that molecular behavior, but such calculations are very tedious for pure substances and should be repeated for each new combination of substances for mixtures. Therefore, from the times when molecular theory was not known until present simpler generalizations were developed and used. Though they have limits of applicability they are still very useful.

The first correlation of properties was stated by ´Emile Clapeyron in 1834 and called ideal gas equation of state (EoS):

P ν = RT, (2.1)

where P is a pressure, ν the molar volume, R the universal gas constant and T the temperature. Equation (2.1) can be rewritten as

P ν

RT = 1. (2.2)

Often instead of dependency of pressure on temperature and volume, compressibility Z is used: Z = P V_RT . For the ideal gas it is always equal 1. For the real gases compressibility shows a deviation from the ideal gas model.

The first law of thermodynamics connects the amount of heat supplied to the system Q, work performed by the system on its surroundings W , and the change in the internal energy of a system U :

dU = δQ − δW. (2.3)

The entropy is defined as dS = δQ_T .

For a system with a changing amount of particles equation (2.3) can be written as:

(24)

µ ≡ dU dNT

S,V

is a chemical potential, or an energy of insertion of one particle in the system if entropy and volume are held constant. NT is the number of particles in the system and V is the volume.

Enthalpy of a system is defined as H = U + P V, therefore

dH = dU + P dV + V dP. (2.5)

Gibbs energy G = H − T S = U + P V − T S, dG = V dP − SdT + µdNT, from what follows µ ≡

dG dNT P,T . Another thermodynamic function is the Helmholtz free energy.

A = U − T S. It shows the maximal work which can be done by the closed system at constant temperature. Its differential is

dA = −P dV − SdT and for systems with different number of particles

dA = −P dV − SdT − µdNT. (2.6)

From (2.6) follows that µ = dA_dn_{T ,V} .

Lewis relation connects chemical potential of the substance i µi and parameters of a system (pressure, temperature):

µi = µi° + RT ln(fi) (2.7)

where µi° is the reference chemical potential of the ideal system, and fi is the fugacity or an effective pressure of the real gas. Fugacity is equal to the pressure of an ideal gas which has the same chemical potential as the real gas. In other words, the ideal gas EoS can be used for real gas, if fugacity is used instead of pressure.

The fugacity coefficient ϕ is defined as a relation between fugacity and pressure:

ϕ = f

P. (2.8)

Obviously for an ideal gas ϕ ≡ 1.

For the i-th substance in the mixture equation (2.8) takes the form: ϕi =

fi yiP

. (2.9)

Calculation of the two-phase equilibrium can be done through fugacities. This is based on the fact, that when liquid and vapor are in the equilibrium, each component has equal fugacities of liquid and vapor phases.

The fugacity coefficient can also be calculated through other thermodynamic functions, which in turn can be found using EoS.

From (2.7) follows µi− µi° = RT ln f_i fi° . (2.10)

The superscript ° refers to the reference ideal-gas state. However from (2.9) follows that fi° = P °yi, thus RT ln ϕi = RT ln fi P°yi = ∂ ∂Ni (A − A°)T ,V,Nj6=i. (2.11)

(25)

From (2.6) follows that at constant temperature and composition dA = −P dv. Therefore residual Helmholtz energy Ar _{= A − A}_{° is an integral from the reference molar volume v° to} the real molar volume v:

A − A° = − ˆ v v° pdv = − ˆ v ∞ P −RT ν dv − RT ln v v°. (2.12)

Combining (2.12) and (2.11) we get RT ln ϕi = RT ln fi P°yi = − ˆ V ∞ ∂P ∂ni T,V,nj6=i − RT V ! dV − RT ln Z. (2.13)

2.4 Choosing of Equation of State

The ideal-gas EoS is based on the assumption, that all molecules move randomly and do not interact. This equation is simple but generally fails at lower temperatures or higher pressures. It fails for heavy molecules and does not allow phase transition. Since then dozens of more sophisticated equations of state have been proposed. They can be broadly divided into four groups: beginning from relatively simple ones (Peng-Robinson, Redlich-Kwong) to more complex ones (Patel, Zabaloy and Vera) and perturbation models. A detailed review and comparison of equations of state and mixing rules can be found in [24–28]. Table 4 gives a summary of all equations that have been considered in this work.

The virial equation of state has a theoretical basis in statistical mechanics, its coeffi-cients relate to the intermolecular forces and can be directly derived:

Z = 1 + ¯ B v + ¯ C v2 + ¯ D v3 + .... (2.14)

With appropriate assumptions about the mathematical form of intermolecular forces, theoretical expressions can be developed for each of the coefficients in second and other terms. In this case the second virial coefficient ¯B corresponds to the interactions between pairs of molecules, the third virial coefficient ¯C to triplets, and so on. Accuracy can be increased indefinitely by considering higher order terms.

The general form of the cubic equation of state is:

P = RT

v − b−

α

v2_{+ γv + ε}, (2.15)

where α, b, γ, ε are parameters.

Various modifications of cubic equations (Van der Waals, Soave-Redlich-Kwong, Peng-Robinson) are widely used because all equation parameters can be easily evaluated from critical temperatures and pressures of studied substances. Martin [29] and Trebble and Bishnoi [30] carried out a detailed comparison of different cubic equations of state.

The general form of the Rubin (BWR) and Modified Benedict-Webb-Rubin equations of state is:

Z = 1 + f1(T ) v + f2(T ) v2 + f3(T ) v3 + f4(T ) (γ1+ γ2/V2) vm exp(− γ2 V2), (2.16)

where the functions fi can contain up to 30 additional parameters. It provides a high level of equation flexibility and accuracy. This resulted in the development of many modifications,

(26)

Table 4: Comparison of main types of equations of state.

EoS Advantage Disadvantage

Virial equation of

state • Relatively simple

• Not precise at high pressure • Coefficients are very difficult to determine

• Usually can be applied only for single-phase gas system

Cubic

• Simple

• Analytic expressions for density or fugacity are not too

demanding

• Adequate, but does not provide perfect accuracy • For polar and associating substances 4 and more parameters should be used • Inapplicable for mixtures of fluids with big differences in molecular size

BWR and MBWR • Good precision

• Many parameters should be calculated

• Results depend strongly on errors in experimental data used for the parameter evaluation • Generalized modifications (with smaller number of parameters) are not accurate Based on

Perturbation theory (Perturbed hard chain, Statistical associated Fluid, Perturbed hard chain of spheres)

• Good precision • Applicable for mixtures of

substances with big difference in

molecular size

• Applicable for high pressures

• Very complicated • For mixtures of unlike

molecules 3 or 4 pure-component parameters are needed

which are widely used for engineering purposes. For example the Span-Wagner modification [31]: p ρRT = IP ol X i=1 niτtiξdi+ IP ol+IExp X i=IP ol niτtiξdiexp(−γξ2), (2.17) where values of IP ol ,IExp, ti , and di are defined in a “trial and error” procedure.

In the perturbation models a compressibility factor Z is found using thermodynamic partial derivation of residual Helmholtz energy Ar _[24]:

Z = 1 − v∂(A r_{/RT )}

∂v |τ . (2.18)

The molecules are considered as a collection of functional groups, the residual Helmholtz energy is written as a sum of reference and perturbation terms

(27)

where perturbation terms can be obtained from rigorous theoretical calculations. For example in Statistical Associated-Fluid Theory, the residual Helmholtz energy consists of Helmholtz energies of an ideal gas Aideal, the spherical segments Asegm, the chains Achain, to represent non-spherical molecule as a chain of spherical segments, and associations Aassoc [18, 32], to take into account all association sites of the molecule:

A = Aideal + Asegm + Achain+ Aassoc. (2.20)

Taking into consideration all advantages and disadvantages of the above-mentioned equa-tions of state, the Redlich-Kwong-Soave modification of cubic EoS is chosen for calculation of the phase equilibrium. It combines good accuracy for large hydrocarbon molecules and simplicity, which allows analytical expression for fugacity.

2.5 Choosing Mixing Rules

There is a large number of different mixing rules for the EoS parameters: beginning from simple Van der Waals rules to complicated relations containing Gibbs energy. For example Van der Waals mixing rule for an arbitrary parameter Ψ is:

Ψmix = n X i=1 n X j=1 yiyjΨij, (2.21)

where the cross parameter Ψij is defined by combining rules:

Ψij =

Ψii+ Ψjj

2 (1 − lij) (2.22)

if Ψ = b from (2.15), and

Ψij =pΨiiΨjj(1 − kij) (2.23)

if Ψ = α from (2.15). Parameters lij and kij are found by fitting to the experimental data, if they are not available they are usually considered to be zero. It is usually kij ≡ kji. This mixing rule also has many composition-dependent modifications [25–28, 33, 34]. If for Ψ = b all coefficients lij = 0 in (2.22), then the mixing rule (2.21) can be simplified:

Ψmix= n X

i=1

yiΨi (2.24)

For a more precise result or in the case of strong non-ideality of substances Excess Gibbs Energy (GE) mixing rules should be applied [35]. A detailed review and comparison of equations of state and mixing rules can be found in [18, 24, 32, 35].

Hydrocarbons used in the reference blends rarely have sizes that differ by a factor ex-ceeding two. Hydrocarbons are non-polar molecules with well established thermodynamical parameters. For these reasons the simple cubic equation of state (Soave modification of Redlich-Kwong EoS) and Van der Waals mixing rules (2.21) and (2.23) for α and (2.24) for b are used to calculate the two-phase diagram with sufficient accuracy.

(28)

2.6 Calculation of Reference Fuel Evaporation Characteristics

The method for developing a phase diagram calculation presented here is based on works of Michelsen [36, 37] and Asselineau et al. [38]. In this thesis the Soave-Redlich-Kwong (SRK) equation of state is used:

P = RT

v − b− α

v(v + b). (2.25)

The system is supposed to have N components. The system is in equilibrium, when the fugacities of the liquid and vapor phases are equal. Fugacities are found from fugacity coefficients, which are calculated from equation (2.11) and the Soave-Redlich-Kwong equation of state (2.25). The following expression for the fugacity coefficient ϕi is obtained:

ln(ϕi) = −ln z 1 − bmix v + bi bmix (z − 1)+ + αmix bmixRT bi bmix − 2 PN j=1 √ αiαjxj amix ! ln 1 + bmix v . (2.26)

The full derivation of this formula is given in the Appendix 1.

Parameters α and b from (2.25) for a mixture of several components are defined with following mixing rules:

αL_mix= N X i=1 N X j=1 xixjαij, bLmix= N X i=1 xibi, (2.27)

for the liquid phase, and

αV_mix= N X i=1 N X j=1 yiyjαij, bVmix= N X i=1 yibi, (2.28)

for the gas phase (vapor). Coefficients of the pure components [39] are: αii= 0.42747 · ξi· R2T_crit,i2 P2 crit,i , bii= 0.08664 RTcrit,i Pcrit,i , (2.29) where ξi = " 1 + 1 − κi r T Tcrit !#2 , (2.30) κi = 0.48508 + 1.5571ωi− 0.15613ωi2, (2.31)

and cross term in the mixing rule for parameter α in (2.27), (2.28): αij =

√

αiiαjj(1 − kij). (2.32)

The parameter kij is found by minimizing the discrepancy of experimental and calculated thermodynamic data. The strategy for the two-phase equilibrium determination can be explained as follows: from given composition of the liquid phase xi and system pressure P one has to calculate the composition of gas phase yi and the corresponding equilibrium

(29)

temperature T . This problem is called bubble point T problem. We have N + 3 unknowns: the gas phase composition yi, i = 1..N , the temperature T and the compressibilities of the liquid and gas phases - zL, zV, and N + 3 equations: two equations of state: for liquid and gas phase, N equilibrium conditions for each substance in the mixture:

f_iV = yiϕVi P = xiϕLi P = f L

i , i = 1...N (2.33)

and one equation for concentrations yi: N X

i=1

yi = 1. (2.34)

Now we can form a system of N + 3 equations to be solved:

F( ˜V) = 0, (2.35)

where the vector of variables ˜VT _{= (T, z}L_{, z}V_{, y}

1, y2, ..., yN), and vector of functions is F1 = P − _v−bRTL mix + αLmix v(v+bL mix) F2 = P − _v−bRTV mix + αVmix v(v+bV mix) F3 = ΦL1x1 − ΦV1y1 ... FN +1= ΦLN −1xN −1− ΦVN −1yN −1 FN +2= ΦLNxN − ΦVNyN FN +3= 1 − PN j=1yj. (2.36)

This system of equations is solved with the Newton-Raphson method and the equilibrium temperature, compressibilities and gas phase composition are obtained for the given pressure and composition of the liquid phase.

2.7 Calculation of Reference Fuel Critical Points

The calculation of critical points is based on the works of Heidemann and Khalil [40] and Michelsen and Heidemann [41]. The critical point is placed on the stability limit of the stable phase. This limit is determined by expanding the expression for the Helmholtz free energy in a Taylor series around a test point:

A − A0− N X i (∂A/∂ni)4ni = = 1 2 N X j N X i ∂2A ∂ni∂nj 4ni4nj+ 1 6 N X k N X j N X i ∂3A ∂ni∂nj∂nk 4ni4nj4nk+ O(4n4). (2.37)

The system is stable if the expression above is positive for all possible deviations of the concentration 4ni. This is assured if the quadratic form (the first term on the right hand side of the equation 2.37) is positive-definite. On the limit of stability the quadratic form is positive semidefinite and the stability is defined by cubic and further terms. The system is

(30)

stable if its stability is defined by a term of an even order. Therefore, the critical point is reached if the cubic term equals zero. Derivatives of the Helmholtz energy can be expressed through the fugacity

∂2_A ∂ni∂nj = RT∂ ln ϕi ∂nj (2.38) and ∂3A ∂ni∂nj∂nk = RT∂ 2_{ln ϕ} i ∂nj∂nk (2.39) Therefore, the system lies on the limit of stability when the quadratic term is positive semidefinite, or what is the same

det(Λ) = 0, (2.40)

where the matrix Λ is formed by elements Λij = ∂2A ∂ni∂nj |T,ν= RT ∂ϕi ∂nj |T ,ν (2.41)

This is equivalent to that the system of linear equations

Λ4n = 0 (2.42)

has a nonzero solution 4n = (4n1, 4n2, ..., 4nN)T.

The definition of the critical point requires that with 4n found from (2.42) and inserted into (2.37) the first non-vanishing term should be of even power, therefore the cubic term should be equal to zero.

C = N X i N X j N X k ( ∂ 3_A ∂ni∂nj∂nk )4ni4nj4nk = 0. (2.43)

Constants of pure components and mixing rules are similar to those used for the phase diagram construction. For the case of the Soave-Redlich-Kwong EoS, elements of matrix Λ become ΛijnT = RT δij xi + (βi+ βj)F1+ βiβjF12 + +amix bmix βiβjF3− aij amix F5+ (βiβj− αiβj − αjβi)F6 , (2.44) where the coefficients are

αi = 1 amix X xjaij, βi = bi bmix . (2.45)

And the cubic term (2.43) can be calculated using the following expression: Cn2_T = RT X∆n3_i x2 i + 3N (βF1)2+ 2(βF1)3 + +amix bmix 3β2(2α − β)(F3+ F6) − 2β 3 F4− 3βaF6 (2.46) where

(31)

α =P i∆niαi, β = P i∆niβi, a = _a1 mix P i P j∆ni∆njaij, N = P i∆ni, and F1 = 1 K − 1, F2 = 2 K + 1, F3 = 2 1 K + 1 2 , F4 = 2 1 K + 1 3 , F5 = 2 ln K + 1 K , F6 = F2− F5, K = v bmix . (2.47) The full derivation of the formulas (2.44) and (2.46) is given in the Appendix 2.

Fig. 2.5 (b) shows a schematic flowchart for the critical points calculation procedure. Such algorithm allows for the calculation of all critical points in the area of interest. For mixtures consisting of a few (up to 5) substances a calculation takes several seconds on an Intel Pentium 4 3.0 GHz, whereas for complicated mixtures (more than 20 substances) it can reach up to a minute.

2.8 Distillation Curve

A distillation curve shows the dependency of the boiling temperature on the evaporated fraction of the liquid. Since the more volatile substances boil at lower temperatures they evaporate first. Thus the boiling temperature for a complicated fuel grows with the growth of the distillate volume. Distillation curves can be easily predicted using a phase diagram algorithm. The equation for the material balance for species i in the liquid phase is:

d( ¯M xi) = − ¯V yi (2.48)

where ¯M is the total number of moles in the liquid phase and ¯V is the molar evaporation rate. The first boiling point and vapor composition y_i(0) are found for the initial mixture x(0)_i . Assuming that ¯V ¯M , the composition for the k-th iteration x(k)_i is recalculated using the equilibrium vapor phase composition y_i(k−1) from the previous iteration:

¯

M (x(k)_i − x(k−1)_i ) = − ¯V y_i(k−1).

2.9 Numerical Algorithms

Several programs are written for calculating the fuel properties described. The calculation of the phase diagram is performed with the Newton-Raphson method. The Fig. 2.5 shows the flow-chart diagrams of the programs.

First for the calculation of the phase diagram data are read: range of pressures at which equilibrium should be calculated, composition of the liquid phase xi, critical temperature Tcriti, pressure Pcriti, and acentric factor ωi of each pure component. The Newton-Raphson

method needs a first approximate solution. Mollerup [37] suggested a good approximation for the y

(0) i

xi ratio in hydrocarbons systems:

Ki = y(0)_i xi = Pcriti P exp 5.42 1 − Tcriti T (2.49)

(32)

Take P Read liquid phase

composition xi,

properties of pure substances Tcriti, Pcriti, ωi,

Assume y(0) Calculate F, J= Jacobian F Solve Divide by 2 No Take previous solution as initial point Print Results. Take next P Yes First iteration ? Yes

Are new values reasonable ? Yes No No 1 J− ∆ = −V ⋅F ur ur new old = + ∆ V V V ur ur ur toler toler v f < ∆ < V F r r V r ∆ (a) Calc b_mix v = 1.05·b_mix to 4

·

b_mix, step 0.001·b_mix Calculate Det(Λ) Calculate C No Find ∆n

Find points where |C| = 0 Yes |Det(Λ)|<ε1 Det( ) T T Det '( )T crit= crit− Λ Λ Calculate Pcrit Print Results Save T_j,P_j,v_j,C_j Take P Read liquid phase

composition x_i, properties of pure substances Tcrit_i, Pcrit_i, ω_i

(b)

Figure 2.5: Schematic flowchart for the calculation procedure of the phase diagram (a) and critical properties of a mixture (b).

Then follows the standard Newton-Raphson method for the solution of a nonlinear system of equations. The program calculates the Jacobian J and the step 4 ~V, if a value of the function F in the new ~Vnew = ~Vold+ 4−→V point is closer to zero than in the previous point, then the step is taken, otherwise 4−→V is divided by a factor of 2. When the function F equals zero with sufficient precision, equilibrium is calculated for the next pressure, and the previous solution is used as a starting solution for the Newton-Raphson iterations.

The calculation of the critical points is carried out by mathematical calculation of the quadratic term in the Taylor series decomposition of the Helmholtz energy at all reasonable values of the mixture molar volume ν. To cover all possible molar volumes the loop starts from v = 1.05 bmix to 4 bmix with a sufficiently small step 0.001 bmix, where bmix is the coefficient from EoS calculated for the mixture. From equation (2.42) the vector 4n is found. It is then used in the calculation of the cubic term C. Points where C = 0 are critical points. This method can provide several solutions Pcrit, Tcrit, that corresponds to the mathematical

(33)

0 20 40 60 80 100 120 140 120 150 180 210 240 270 300 330 T, K P , a tm

calculation, this work calculation, REFPROP experiment, Boukouvalas et al.

Figure 2.6: Phase diagrams for an H2S-CO2-N2-CH4-C2H6-C3H8-nC4H10-isoC4H10-C5H12 mixture. Comparison of the experimental data of Boukouvalas et al. [42] and calculations using the proposed algorithm and the REFPROP© [43] program.

solution of the system of equations. If present, false critical points can be easily identified by unphysical values of temperature or pressure (negative).

The program for the distillation curve calculation uses the subroutines of the phase dia-gram evaluation. First, the composition of the liquid phase, pressure and critical properties of the pure components are read. Given one mole of the liquid phase, a small amount of substance 4µ is evaporated and the composition of the vapor is calculated from the initial mixture composition. Then composition of the liquid phase is adjusted, according to the amount of substance evaporated. These iteration is repeated until no liquid remains.

2.10 Results of Surrogate Fuel Modeling

For the validation a variety of mixtures and experimental data [44–49] is used to verify the accuracy of the algorithms described above. In general, the SRK EoS shows good correla-tions for hydrocarbons. Significant discrepancies can be observed for highly polar substances (methanol-water mixture) and hydrocarbons with big size differences (methane-decane, methane-propylbenzene). Nevertheless, the calculation of an ethane-methylcyclohexane mix-ture provides satisfactory results.

Table 5 shows the correlation of calculations with experimental data [46] for propane-cyclohexane. The calculations performed for propane-toluene and propane-m-xylene mixtures [46] demonstrate the same conformity with the experimental data. For the calculations, the pressure and liquid phase composition are preset, equilibrium temperature and vapor phase composition are predicted. Mixing rules parameters kij are listed in Table 5 along with the results obtained using the REFPROP© code [50] provided by the National Institute of Standards and Technology (NIST). This code uses the Span-Wagner modification [31] of the BWR EoS (2.17). REFPROP© has an extensive but non-exhaustive list of species and

(34)

0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 10 20 30 40 50 60 70

80 _{experiment, Laugier et al.}

calculation, this work

P , a tm H 2S mole fraction

Figure 2.7: Equilibrium pressure vs. H2S mole fraction in an H2S-C6H14-C15H32 mixture. Comparison of the experiment of Laugier et al. [44] and calculation. k12=0.07, k13= k23 = 0.

0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 0 100 200 300 400 500 40 50 60 70 80 90 100 P c rit _{, a} tm T c ri t , K

Butane mole fraction T_crit, experiment, Olds et al. T_crit, calculation, this work

P_crit, experiment, Olds et al. P_crit, calculation, this work

Figure 2.8: Calculated and experimental [45] dependencies of critical pressure and tempera-ture for n-C4H10-CO2 system. Mixing rule parameter k12 = 0.06.

(35)

0 20 40 60 80 100 460 480 500 520 540

C10(75%)-C14(25%) experiment, 83 kPa, Huber et al. C10(75%)-C14(25%) calculation, 83 kPa, this work C10(50%)-C14(50%) experiment, 83 kPa Huber et al. C10(50%)-C14(50%) calculation, 83 kPa this work C10(50%)-C14(50%) experiment, 70 kPa Huber et al. C10(50%)-C14(50%) calculation, 70 kPa this work

T e m p e ra tu re , K

Distillate volume fraction, %

Figure 2.9: Experimental and calculated distillation curves for 75/25 and 50/50 mole fraction mixtures of C10H22-C14H30at P = 83kPa and P = 70 kPa. Experimental data and surrogate mixture are from Huber et al. [17].

Table 5: Equilibrium concentrations in the propane(1)-cyclohexane(2) mixture; experi-mental data are compared with calculations using the presented algorithm and the NIST REFPROP© [43] package. Mixing rule coefficients k12= k21= 0.016.

P(atm) Texp(K) Tcalc(K) TNIST(K) x1 exp y1 exp y1 calc y1 NIST 2.52 5.47 8.44 11.45 13.01 9.18 12.44 16.09 21.42 26.05 38.29 17.76 28.62 46.98 58.23 313 313 313 313 313 393 393 393 393 393 393 473 473 473 473 312.4 311.2 311.9 312.8 312.1 389.5 389.1 388.8 390.2 391.6 470.2 470.4 473.6 475.7 281 296 306 310.5 313 354 357 360.5 373 382 450 433 438 442 0.152 0.378 0.615 0.847 0.97 0.138 0.208 0.287 0.395 0.491 0.731 0.054 0.175 0.372 0.515 0.906 0.962 0.9813 0.9922 0.9983 0.68 0.757 0.812 0.855 0.881 0.908 0.232 0.49 0.619 0.608 0.8286 0.9427 0.9776 0.9934 0.9989 0.7016 0.7928 0.8538 0.932 0.975 0.349 0.671 0.852 0.913 0.96 0.97 0.98 0.99 0.999 0.805 0.84 0.886 0.88 0.91 0.3 0.58

Development of a reduced chemical-kinetic combustion model for practical fuels