ARTICLE pubs.acs.org/EF
Kinetics Oxidation of Heavy Oil. 1. Compositional and Full Equation of State Model Alexandre Lapene,*,†,‡ Gerald Debenest,§,^ Michel Quintard,§,^ Louis M. Castanier,† Margot G. Gerritsen,† and Anthony R. Kovscek† †
Energy Resources Engineering, Stanford University, Stanford, California 94305, United States Universite de Toulouse; INPT, UPS; IMFT (Institut de Mecanique des Fluides de Toulouse); Allee Camille Soula, F-31400 Toulouse, France ^ CNRS; IMFT; F-31400 Toulouse, France §
ABSTRACT: In-situ combustion (ISC) is an enhanced oil recovery method that involves coupled mass transport, fluid phase equilibria, and chemical reactions. Because of its complexity, a full understanding of ISC process phenomena and their modeling remains an open research topic. Ramped temperature oxidation (RTO) of crude oil is one of the methods to obtain the kinetics of oxidation. This paper provides a new method to interpret such RTO experiments and then gives an improved physical understanding compared with classical analytical methods. The reactor model is based on a compositional and full equation of state formulation. The numerical model is applied to a synthetic test case. Although generally not taken into account, this work shows that spatial transport effects are present despite the plug flow design of typical RTO reactors. Moreover, coupling between thermodynamics and chemistry is important even for extra-heavy oil that in a first approximation could be seen as totally nonvolatile.
’ INTRODUCTION In-situ Combustion. In-situ combustion (ISC) is a thermally enhanced oil recovery method. Enhanced oil recovery is broadly described as a group of techniques used to extract crude oil from the subsurface by the injection of substances that are not originally present in the reservoir with or without the introduction of extraneous energy.1 During ISC, a thin combustion front is propagated through the reservoir driven by injected air. The heat generated results in greater temperatures leading to a reduction in oil viscosity and an increase in oil mobility. In-situ combustion is a complex reactive transport process. Under the combined effects of phase change, chemical reaction, andheat and mass transfer, distinct regions are identified. Burger et al.2 distinguished four regions in the case of regular ISC, Figure 1, for a one-dimensional case. In zone 1, Figure 1, combustion has already taken place and all the fuel is consumed. The cold injected gas heats up upon contact with the porous medium and recovers a small amount of the heat released by the combustion. Zone 2, Figure 1, is the combustion zone. Heterogeneous and homogeneous energetic bond-breaking reactions occur between oxygen and residual hydrocarbon. The physical dimensions of the combustion front are of order centimeters. In zone 3, Figure 1, the heavy fractions that have been neither displaced nor vaporized are cracked. Finally, in zone 4, Figure 1, the temperature is too low for further chemical reaction to take place. Oil and water are displaced by the effluent gases. In the region nearest the reaction zones, successive vaporizations and recondensations of the lighter fractions of the oil take place. The same is observed for water originally present in the matrix. In the region where the temperature is lower than the condensation temperature of water, a water bank is observed. The water saturation is greater than the original water saturation. Ahead r 2011 American Chemical Society
Figure 1. Temperature and saturation profiles according to Burger et al.2 This figure has been adapted from reference 2.
of the cold water, an oil bank appears, with an oil saturation greater than the original oil saturation. Mechanistic ISC modeling requires a high-quality description of all the phenomena above. It is apparent that many of the phenomena occur near the combustion front. Thus, accurate reaction modeling is critical. Poor reaction modeling affects the modeling of all other physics. One of the challenging features in the development of reaction models for an ISC process is the complexity of the reactive system. Crude oil contains a large number of components, and consequently, reactions cannot be characterized fully. The oil Received: March 8, 2011 Revised: September 21, 2011 Published: October 07, 2011 4886
dx.doi.org/10.1021/ef200365y | Energy Fuels 2011, 25, 4886–4895
Energy & Fuels
ARTICLE
Figure 2. Schematic diagram of a ramped temperature oxidation reaction kinetic apparatus from the work of Lapene et al.11 This figure has been adapted from reference 11.
description and the reaction model have to be simplified while retaining the critical physics: the reaction model should give accurate reaction rates under a range of various conditions, such as variable initial water saturation and variable heating rate. To improve our understanding of the chemistry and aid the modeling, laboratory experiments are essential. The choice of experiments and their interpretation remain a key point in a modeling strategy. Kinetic Experiments. Various experimental devices are used for chemical studies of oil. We choose to classify these experiments depending on the scale. The combustion tube is a meter scale experiment. Prasad and Slater3 described it as a portion of the reservoir, simulated at the lab scale. Its objective is to mimic the full process. It is especially used as a global tool for ISC study in order to design and optimize field projects. Many experiments have been performed.46 Combustion tubes can also be used to study oil oxidation kinetics but have limited value in this context. It is generally too difficult to distinguish, analyze, and model chemical reactions among all the phenomena present in the tube. Instead, specialized smaller scale experiments are used. One such experiment is the thermogravimetric experiment. Thermogravimetric experiments are small scale experiments. Their main objective is to study the mass variation as a function of the time. A sample of the test material (oil or oil and rock) is placed into a small cup that is supported on, or suspended to, an analytical balance located outside the furnace chamber. The balance is zeroed, and the sample cup is heated according to a predetermined thermal cycle. Gas flows around the sample. The experimental apparatus can be coupled with a gas analyzer in order to examine the composition of the exhausted gases. Unfortunately, such experiments are not representative of porous media processes. The size of the sample does not make it possible to study reactions and take into account surface and catalytic effects. As a consequence, kinetic results obtained from these experiments are hard to upscale to combustion tube or field scale. Kinetic cell experiments are typically used to dinstinguish the critical chemical reactions that take place in porous media.4,711 A kinetic cell is a small metallic cylinder filled by crude oil, clay, water, and sand to be studied. It is located in a temperature
controlled furnace. Gas is injected at constant flux and continuously analyzed at the outlet. Figure 2 gives a schematic diagram of the reaction kinetic apparatus used by Lapene et al.11 It is a typical RTO setup. Kinetic Cell Model. Although the kinetic cell experiment is simplified compared to full ISC, modeling of the observed processes remains a challenge. The simplest modeling approach is analytical. In 1984, Fassihi et al.9 proposed a method, initially introduced by Weijdema.7 They assumed an Arrhenius model for the rate constant:
dCf ¼ AeEa =RT P02 a Cf b dt
ð1Þ
where b and a represent, respectively, the reaction orders for carbon and oxygen partial pressure, Ea is the activation energy (J/mol), T is the temperature (K), R is the gas constant (J/(mol K)), and A is a pre-exponential factor sometimes referred to as the Arrhenius constant. Similar forms are used by other authors.10 To enable analytical solutions, strong simplifying assumptions are made: • The reactions are homogeneous in the cell, i.e., there is no gradient effect. • Mass variation in the cell is caused only by chemical reactions, i.e., evaporation and condensation are not taken into account. • Variations of oxygen, carbon dioxide, and carbon monoxide are only due to chemical reactions, i.e., the presence of other gases as cracking gases or evaporated components is neglected. • The inlet and outlet air fluxes are the same, i.e., variation of mass within the gas phase is not taken into account. Kristensen et al.12 developed a numerical tool, called a “virtual kinetic cell”, especially designed for kinetic cell experiments. It rests on the solution of an ordinary differential equation system that describes chemical reaction and phase change modeling. The latter is computed by flash calculations.13,14 Spatial effects are not taken into account, and mass transfer is modeled through a simple volume integrated flux balance and a valve coefficient. When modeling a reaction kinetics apparatus, a numerical approach allows us to incorporate coupled nonlinear physics. Analytical methods, have been very useful while computer 4887
dx.doi.org/10.1021/ef200365y |Energy Fuels 2011, 25, 4886–4895
Energy & Fuels
ARTICLE
performance was poor, but they are limited for modeling of crude-oil kinetics cells because of the strong assumptions required as discussed previously. Accordingly, this paper presents a new fully compositional method for kinetic cell modeling. Because of the characteristic dimension of the kinetic cell (centimeters), we believe that the reaction is not homogeneous in the reactor. Transport gradient effects are taken into account by a full equation of state (EoS) compositional model. A full EoS also accurately accounts for phase behavior. This is important because of the coupling between chemistry and thermodynamics. Analytical methods that do not consider phase changes overestimate the mass of the fuel, and so, the kinetic parameters are not accurate. The fuel laid down is the organic matter remaining after oil evaporation and cracking, and not the total initial fuel. Finally, this paper is also part of a broader study from Lapene et al.15 on the interpretation of kinetic data from kinetic cell experiments to obtain reaction model parameters. These results, when put together, constitute a detailed methodology for kinetic cell modeling and data postprocessing.
’ SIMULATOR DESCRIPTION The kinetic cell experiment modeled is that illustrated in Figure 2. Mathematical Description. Usually, samples are composed of sand, clay, water, and oil. The water-rich phase is neglected because, in operating pressure conditions, the saturation temperature of water is much lower than the temperature at which first reactions are observed. Sand, clay, and solid components created by reactions constitute the solid phase. Oil and gas are mobile and miscible phases. The system is composed of Nc mobile components and Ns solid components. The temperature of the furnace is controlled during the experiment and recorded within the cell. We could solve the heat equation within the cell considering the furnace power as a boundary condition and then get the temperature field of the system. This calculation, however, requires many data to be accurate. First of all, we need to solve the heat equation over the whole system, including the furnace, the kinetic cell, and the sample. Many inputs are needed for computations, such as the thermal properties of the materials, the natural convection coefficient of the air, and so on. Moreover, solving the heat equation means that we have to deal with the heats of reactions and enthalpy of vaporization. Also, it is difficult to quantify and describe in the model how the temperature controller manages the furnace power. Because of all these reasons, we decide to take the recorded temperature history within the cell as an input. The temperature throughout the cell is taken as uniform. To verify the relevance of this assumption, several thermocouples were placed at various locations within the cell during a test. The thermocouples showed a near uniform temperature. The total mass conservation of the system is given by ∂ ½ϕFðMwo O þ Mwg V Þ þ Mws ξs þ ∇ðFo B vo ∂t
þ Fg B v gÞ ¼ 0
ð2Þ
where ϕ is the porosity, F is the global molar density, Mwβ is the molecular weight of the β phase, O is the oil molar fraction, V is the vapor (or gas) molar fraction, ξs is the molar density of the solid phase, Fβ is the density of the β phase, andBv β is the filtration
velocity of the β phase. The global molar density is defined as Fg F Sg þ o So ð3Þ F ¼ Mwg Mwo where Sβ is the saturation of the β phase. The mass conservation of each mobile component in the system, if molecular diffusion is neglected, is expressed as ∂ ðϕFMi zi Þ þ ∇ðFo B v o Yo, i þ Fg B v g Yg, i Þ ¼ q_ i , ∂t i ¼ 1, :::, Nc
ð4Þ
where Yβi is the mass fraction of the ith component in the β phase and zα is the global molar fraction of ith component in the pore space. The mass conservation of each solid component in the system is given by ∂ ðCs, i Mws, i Þ ¼ q_ s, i , ∂t
i ¼ 1, :::, Ns
ð5Þ
where Cs,i is the concentration of the ith solid component in the system. The macroscopic velocity of oil and gas phases is described by the generalized multiphase Darcy’s law krj w v j ¼ K ð∇Pj Fj gBÞ, B μj
j ¼ fo, gg
ð6Þ
where krβ is the relative permeability of the β phase, μβ is the w viscosity of the β phase, K is the absolute permeability tensor, Pβ the pressure in the β phase, and B g is the gravitational acceleration. Thermodynamic equilibrium is expressed by equality of fugacity of each mobile component in each mobile phase g
fio ðx o , T, Po Þ ¼ fi ðx g , T, Pg Þ,
i ¼ 1, :::, Nc
ð7Þ
where fiβ is the fugacity of the ith component in the β phase, xβ are the molar fractions of all the components in the β phase, and T is the temperature. Fugacities are calculated by the Peng Robinson EoS16 using classic mixing rules.17 The molar volumes of the oil and gas phases are evaluated as ~vj ¼
Zj RT , Pj
j ¼ fo, gg
ð8Þ
with ~vj as the molar volume of the phase j and Zj as the compressibility factor of phase j. Because of the poor predictability of the PengRobison EoS for molar volumes of liquid phases, the last equation is corrected (See the Appendix for details). The capillary pressure drop between the oil and gas phase is neglected, so Po ¼ Pg
ð9Þ
The set of Nr reactions occurring in the system is
∑
i∈R
j
ν i Ci f j
∑
i∈P
j
ν i Ci
ð10Þ
j
where Nr is the number of reactions, R j is the set of the index of reactant for the jth reaction, including mobile and solid species, P j is the set of the index of products for the jth reaction, including mobile and solid species, νji is the stoichiometric coefficient for reaction j and component i (negative for reactants and positive for products), and Ci is the ith component (mobile or solid). We note that "j, R j ∩ P j = L. 4888
dx.doi.org/10.1021/ef200365y |Energy Fuels 2011, 25, 4886–4895
Energy & Fuels
ARTICLE
We assume that the oxidation reactions are heterogeneous, which means the reaction occurs at the surface of the interface (or within a thin, compared to the pore scale, mixing layer). We assume a reaction dominant system, i.e., a high Damkholer number. At the microscale, the reaction rate of the oxidation reaction can be expressed for j as τj ¼ kj PO2
Y i ∈ R =O2 j
½Ci ni, j
ð11Þ
where kj is the reaction-rate constant of the jth reaction, [Ci] is the total concentration of the ith component in the system (including pore volume and solid volume) and ni,j is the order of the jth reaction and ith component. For the nonoxidation reaction, the reaction rate is τj ¼ kj
Y i ∈ Rj
½Ci ni, j
ð12Þ
The reaction rate constant is given by the Arrhenius law j
kj ¼ kj eðEa =RTÞ j
ð13Þ
Eja
where k and are, respectively, the pre-exponential factor and activation energy of the jth reaction. At the macroscale, which corresponds to the scale we are treating this problem, the expression of the reaction rate is the result of an upscaling stage. We will not introduce here the different steps because it is not the goal of this paper, but we want to underline that in the case of a reaction-dominated pore-scale problem, the macroscale reaction term is the product of the pore-scale reaction rate times the specific surface. Mathematically speaking, it is equivalent to assume that the pre-exponential factor in eq 13 is a function of the specific surface between the reacting phases. While the reactions are occurring and the system evolving, the surface area is changing because of phase change, chemistry, and so on. For such a complex system as the kinetic cell, because of the lack of data, it is impossible to characterize the surface area evolution, which is data that depends on the microscale configuration and phase repartitions. Mamora et al.4 proposed a model to take into account specific surface area variations. They considered the sand grain as a sphere. Unfortunately, this simple approach is not expandable to our problem. If the specific surface area cannot be characterized and then the pre-exponential is taken as a constant, we think it is more relevant to consider the total concentration instead of the partial pressure of oxygen in eq 11. Indeed, for a given pressure and temperature, while the gas quantity is decreasing, the specific surface area is likely to decrease as well and, then, the reaction rate too. In the case of the pre-exponential factor being taken as a constant, the expression of the reaction rate involving the partial pressure of oxygen does not necessarily decrease. Finally, we decided to express the macroscopic reaction rate for oxidation reaction as τj ¼ kj
Y i∈R
j
½Ci ni, j
ð14Þ
The source terms in eqs 4 and 5 are written as q_ k ¼ Mwk
Nr
j νk τ j , ∑ j¼1
k ∈ f1, :::, Nc g or k ∈ f1, :::, Ns g ð15Þ
Numerical Solution. The mathematical problem contains Nc + Ns independent variables. Thus, Nc + Ns equations need
to be solved. We choose to use the global variable formulation as opposed to the natural variable formulation (see the work of Cao18 for details). In this formulation, the primary unknowns are the following: • Pressure P • Nc 1 global mole fractions zi • Ns global concentrations Cs,iThe primary governing equations are the following: • Equation 2 • Equation 4 • Equation 5 Spatial discretization is performed by the finite volume method. Second-order schemes are used for diffusive terms, and first-order upwind schemes are fused or advective terms. Because of the symmetric geometry, injection, and heating conditions, the problem is a one-dimensional problem. Temporal discretization is done by a first-order scheme. The fully implicit method is adopted in order to preserve coupling and unconditional numerical stability. The nonlinear system resulting from the discretization stage is solved by Newton’s method. An automatic time step method is implemented. The time step is allowed to vary between a given range. If the Newton solutions fails, the time step decreases, and if n successive Newton solves are successful, the time step increases. Time step variation is executed following a given ratio. In this special problem, an important stage is the phase equilibrium calculation. Equation 7, complemented by mass fraction constraint equations, exactly defines the phase equilibrium problem. The primary variable formulation makes it possible to solve the problem within each Newton iteration inner loop, as a pressure/temperature negative flash.19,20 Thus, the coupling between phase changes, reactions, and mass transport is fully preserved. A postprocessing stage is added in order to take into account the specificity of the gas separation at the outlet of the cell. The gas composition at the outlet is flashed at 1 atm and 273.15 K, which are also the conditions of the experiment outlet. Test Case: Experimental Situation. We describe the test case used to evaluate the model and code. It is a virtual experiment designed to evaluate the code. Numerical parameters are taken as close as possible to experimental conditions. Numerical Parameters. The oil description characterization arises from an extensive PVT study presented in the Appendix. Table 3 shows the oil composition and properties used for this study. In addition to the sand, the only solid considered is coke. It is assumed to be pure carbon. Other components are water, oxygen, nitrogen, carbon dioxide, and carbon monoxide. The domain is a 0.04 m long and 0.255 cm in diameter cylinder. The outlet pressure is fixed at 6.98 bar, and the initial temperature is 300 K. The air injection rate at standard conditions is 1000 cm3/min. The injected gas is composed of oxygen (20.7%) and nitrogen. There are 2 g of oil in the sample. This corresponds to a saturation of 0.29 at initial conditions. For this test case, the temperature sustains a 0.05 K/s linear heating rate. We propose a reaction mechanism. It is simple in character but reproduces the global behavior observed during experiments. We include low temperature oxidation (LTO) reactions and high temperature oxidation (HTO) reactions. Three kinds of 4889
dx.doi.org/10.1021/ef200365y |Energy Fuels 2011, 25, 4886–4895
Energy & Fuels
ARTICLE
reactions occur. The first is a cracking reaction, in which a part of the oil is transformed into coke and methane. It is given by y y Cx Hy f x coke þ CH4 ð16Þ 4 4 The others are LTO and HTO reactions. Both follow this reaction mechanism 1 xð2m þ 1Þ y x þ z O2 f CO Cx Hy Oz þ 2 m þ 1 2 m þ 1 mx y þ CO2 þ H2 O ð17Þ m þ 1 2 where m is the carbon dioxide and carbon monoxide ratio. LTO occurs at a lower temperature than HTO. Kinetic data related to these reactions are provided in Table 1. Ten Block Simulation. A simulation with 10 blocks is performed. Results are plotted in Figure 3. The reaction model reproduces the classic behavior of heavy-oil combustion. Two peaks are observed in the oxygen consumption oxygen profile. Carbon monoxide and carbon dioxide are produced while oxygen is consumed. Coke and methane are created by the cracking reaction. This reaction is in competition with C50+ combustion. The results are classic, and a sensitivity analysis is performed to determine grid-size effects. Spatial Effects. Figure 4 gives oxygen mole fraction profiles in the gas phase at different times, for the 10 block simulation.
Oxygen mole fraction is not homogeneous in the kinetic cell, especially while strong reactions occur, for example during HTO reactions at 8000 s. Because of the strong reactivity of the system, oxygen coming into the cell reacts with oil. Oxygen is consumed, as it flows, through the cell. As a consequence, the reaction rate, depending on the oxygen concentration, also depends of the location, and the reaction is not homogeneous. Obviously, a 0-dimensional model cannot capture this spatial effect because the concentration is assumed to be the same throughout the domain. Note that the sample size within the cell is only 4 cm long, and the short length prohibits stronger gradients in oxygen concentration along the cell. Mesh Sensitivity Study. A mesh sensitivity study has been performed. Six mesh systems have been tested: 1, 2, 3, 5, 8, and 10 blocks. The time step is fixed at 100 s, in order to make certain that it remains constant for all the runs. According to a sensitivity test not presented here, this time step also insures time convergence, i.e., reducing the time step would not change the results. For comparison, we consider the oxygen consumption as a function of time because it is the main information provided by the experiment. Results are shown in Figure 5. The mesh influences the amplitude of the oxygen consumption peaks and their positions. Spatial effects observed in Figure 4 impact oxygen
Table 1. Kinetic Data for Test Case prereaction
organic
activation
exponential
m
reactant
energy (J)
factor (SI)
ratio
cracking
C50+
90000
2 103
LTO
C2C11
82000
2 104
2
85000 120000
2 102 1 105
4 5
C12C16 C17C21 C22C27 C28C35 C36C49 HTO
C50+ coke
Figure 4. Oxygen mole fraction in the gas phase.
Figure 3. Evolution of oxygen, carbon monoxide, and carbon dioxide (a) and oxygen, methane, and coke (b). 4890
dx.doi.org/10.1021/ef200365y |Energy Fuels 2011, 25, 4886–4895
Energy & Fuels
Figure 5. Oxygen consumption for different mesh sizes.
Figure 6. Oxygen mole fraction in the gas phase and oil composition for the nonvolatile and the flash model.
data observed at the outlet of the cell. Qualitatively, the case with an 8 block mesh seems to be the coarsest mesh to reach a fair convergence because plots for the 8 and 10 block meshes merge. However, for some practical applications, the 0D simulation could be considered as satisfying. A deeper investigation is needed with a more realistic reaction mechanism to conclude. Phase Equilibria Sensitivity. In our model, we use a full EoS formulation. Phase equilibria are calculated according to equality of fugacities as in eq 7. A two-phase flash problem is solved. As previously discussed, analytical models used for kinetic cell experiments do not include phase changes, so oil is considered as nonvolatile. In Figure 6, we compare simulations results for a “nonvolatile” approach and our approach using a flash calculation. In order to get a realistic nonvolatile assumption, oil components from C2C11 to C50+ are assumed to be inert for vaporization and condensation. The comparison clearly shows that oxygen consumption is drastically increased for LTO reactions but remains almost unchanged for HTO reactions. In the nonvolatile model, oil is only consumed by reactions. Components react at the same characteristic temperature. This temperature only depends on the reactions that are involved. The quantity of burned fuel for each reaction is exactly the initial
ARTICLE
Figure 7. Oxygen mole fraction in the gas phase and oil composition for the K-value and the flash model.
amount in the cell. However, because of thermodynamic effects, the quantity of fuel burning depends on the history of the system. If a component evaporates, it is advected out of the cell and its capacity to react is directly linked to its residence time in the system. On the contrary, if it remains liquid, it stays available and the reaction only depends on the reaction rate. Thus, the system involves a strong coupling between chemistry and thermodynamics. Our model handles this coupling naturally. Figure 6 exhibits three relevant situations. The first is characterized by a C2C11 component: it is evaporated and flushed from the reactor by the flow before chemical reactions. In the second, the component consumption is the result of both chemistry and thermodynamics. The C17C21 component illustrates this situation. Finally, the C50+ component is the perfect example for the third situation. It is completely consumed by the reaction before it is evaporated. Its evaporation only slightly affects the end of the reaction. That is why oxygen consumption during HTO reactions is almost the same for both models. Definitely, phase changes impact the simulation results and especially modify the oxygen consumption. The nonvolatile assumption is not suitable for our modeling purposes. For comparison, equilibrium data are also calculated using the K-value approach; see the work of Lapene et al.,20 in which the K-value method is presented. There are generated by Wilson’s correlation,21 which assumed that equilibrium relations only depend on temperature and pressure. Results are shown in Figure 7. The K-value method leads to reduced oxygen consumption for both reactions because the evaporation occurs earlier for the K-value method. Hence, less fuel is available for reactions. Results are sensitive to the approach used for thermodynamic modeling. Coupling between thermodynamics and chemistry is strong. In order to evaluate this coupling and quantify the effect on oxygen consumption, we write the equation of mass conservation of the oxygen component in one dimension. Only one fuel is considered, and oxygen is assumed to be only in the gas phase. Then, the equation is written as ∂ ∂ ½ϕSg CO2 þ ðCO2 vg Þ ¼ νO2 k½Sl CF þ Sg KF CF ϕ2 CO2 ∂t ∂x
ð18Þ
where CO2 is the oxygen concentration in the gas phase, vg is the gas Darcy’s velocity, νO2 is the stoichiometry coefficient in the reaction occurring between oxygen and fuel, CF is the 4891
dx.doi.org/10.1021/ef200365y |Energy Fuels 2011, 25, 4886–4895
Energy & Fuels
ARTICLE
Table 2. Oil Composition by GC Analysis boiling point (C)
concentration of fuel in the oil phase, KF is defined as the ratio between the concentration of fuel in the liquid phase and the one in the gas phase, k is the reaction rate (m3/(mol s)), ϕ is the porosity, Sg is the gas saturation, and Sl is the liquid saturation. Equation 18 is integrated over the cell volume V ðCO2 vg Þ Q_ O2 ∂ ½ϕSg CO2 þ ¼ νO2 k½Sl CF þ Sg KF CF ϕ2 CO2 ∂t V L
ð19Þ
where L is the length of the kinetic cell and Q_ O2 is the molar flux. Variables are nondimensionalized as vg ; vog k k ¼ o ; k CO CO 2 ¼ o 2 ; CO2
vg ¼
t ¼
tvog
L CF CF ¼ o CF Q_ O L Q_ O2 ¼ o 2 o CO2 v g
ð20Þ
The oxygen balance is written as
∂C ϕSg O 2 þ CO2 vg Q_ O2 ¼ kCF CO2 DaKF ∂t
ð21Þ
where DaKF = (νO2ϕ2koCoFL[SgKF + Sl])/(vog ) is a Damkohler number quantifying effects of reaction over advection. If KF f ∞, i.e., all the fuel remains in the liquid phase, DaKF f ∞ and the variation of oxygen concentration is only controlled by the reaction. If KF f 0, i.e., the fuel is only in the gas phase, DaKF f (νO2ϕ2koCoFLSl)/(vog ); in other words, considering experimental parameters, the oxygen concentration is controlled by the residence time L/vog of the fuel in the cell. Figure 8 illustrates the evaluation of log(DaKF) as a function of log(ko) and log(Kf) for experimental conditions: νO2 = 1, ϕ = 0.36, Sg = 0.30, CoF = 138 mol/m3, L = 4 102 m, and vog = 2 102 m/s. Two regions are distinct. The blue region corresponds to DaKF < 1, and the red region, to DaKF > 1. If the reaction is weak, the advection dominates except for high KF, but if the reaction is strong, the reaction always dominates. It is interesting to notice that DaKF can be close to one, or lower, only if the reaction is not strong. Accordingly, the sensitivity to the phase equilibria in test cases is because evaporation occurs before and while reactions occur.
mass
fraction
fraction 0
320.44
320.44
nitrogen N2
0.0004
109.32
109.32
carbon CO2
0.0132
0.0014
76.56
76.56
hydrogen H2S
0
0
259.06
259.06
methane C1
0.1621
0.0063
128.02
128.02
ethane C2
0.0027
0.0002
43.96
43.96
propane C3
0.0038
0.0004
10.94 30.92
10.94 30.92
i-butane i-C4 n-butane n-C4
0.002 0.0035
0.0003 0.0005
82.04
82.04
i-pentane i-C5
0.0015
0.0003
96.98
96.98
n-pentane n-C5
0.0013
0.0002
97
Figure 8. Evaluation of log(DaKF) as a function of ko and Kf.
mole fraction
156
hexanes C6
0.0062
0.0013
156
208.9
heptanes C7
0.0043
0.0011
208.9
258.1
octanes C8
0.0073
0.002
258.1
303.1
nonanes C9
0.0069
0.0022
303.1 345
345 385
decanes C10 undecanes C11
0.0142 0.0197
0.0049 0.007
385
419
dodecanes C12
0.0262
0.0103
419
455
tridecanes C13
0.0296
0.0126
455
486
tetradecanes C14
0.0328
0.0152
486
519.1
pentadecanes C15
0.0299
0.015
519.1
550
hexadecanes C16
0.0305
0.0165
550
557
heptadecanes C17
0.0222
0.0128
557 603
603 626
octadecanes C18 nonadecanes C19
0.0231 0.0248
0.0141 0.0159
626
651.9
eicosanes C20
0.0239
0.016
651.9
675
heneicosanes C21
0.0219
0.0155
675
696.9
docosanes C22
0.0208
0.0154
696.9
716
tricosanes C23
0.0178
0.0138
716
736
tetracosanes C24
0.0165
0.0133
736
755.1
pentacosanes C25
0.016
0.0134
755.1 774.1
774 792
hexacosanes C26 heptacosanes C27
0.0153 0.015
0.0134 0.0137
792.1
809.1
octacosanes C28
0.0138
0.013
809.1
826
nonacosanes C29
0.0123
0.012
C30+
0.3533
0.7187
826
∞
’ CONCLUSION A new kinetic cell simulator for crude-oil oxidation based on a compositional and full equation of state model has been developed. Compared to the classic approaches (analytical or 0D cell), it takes into account phase equilibria and spatial mass transport effects. The fully implicit solution allows display of all the possible couplings. Using this new numerical tool, some important remarks can be made: 1 Spatial heterogeneity of the reaction within the cell was clearly emphasized. This is an important point, because one can not assume that the same reaction occurs throughout the cell even for the convection-dominated regime with an excess of oxygen. It is therefore important to address the problem of transport in this cell and not to consider it as a 0-dimensional object. 2 Our dimensionless study (from eqs 19 to 22) has shown that coupling between chemistry and thermodynamics 4892
dx.doi.org/10.1021/ef200365y |Energy Fuels 2011, 25, 4886–4895
Energy & Fuels
ARTICLE
Figure 9. Extension of oil composition.
can be strong, especially when reaction occurs while evaporation is in progress. Moreover, results of the test case show that considering no evaporation—which is the usual approach—overestimates the quantity of reacting fuel. An approximate method such as the K-value method also exhibits discrepancies with the more-accurate flash approach. Finally, a rough modeling of thermodynamics/chemistry coupling will generate bad evaluations of the reaction scheme because a wrong approximation of the quantity of fuel leads to a wrong estimation of oxygen consumed. Then, in terms of reactive parameters, a first approach could be to misjudge kinetic parameters in Arrhenius functions. This topic is discussed in ref 15. Finally, this tool is a part of a toolbox useful to determine kinetic parameters when heavy-oil combustion takes place in a porous medium. If we couple this tool with an optimization algorithm and if we are able to build an objective function based on real experimental data, we are able to find solutions in terms of activation energies, pre-exponential factors, reaction orders, and other terms for each reaction. This is the idea of the forthcoming paper.15 A reaction scheme and associated reaction parameters will be given.
’ APPENDIX: OIL CHARACTERIZATION Experimental Data. The oil used for this modeling is an extraheavy oil. Gas chromatography (GC) data are provided in Table 2. Data are for a live oil. Oil Composition Extension. The first step of the oil characterization is to extend the composition from C30+ to C50+. This limit corresponds to the validity of correlation used. It is assumed that the composition distribution follows zðnÞ ¼ zmax qn nmax
where zmax is the greatest mole fraction, nmax is the carbon number of the component, and q is a real number. This last one is the fitted parameter to match the shape of the curve. The composition is extended using eq 22. Results are plotted in Figure 9. Evaluation of Properties. The molecular weights of analyzed components are calculated considering the mole fraction and mass fraction. For extended components, these properties are extrapolated using Mwi ðnÞ ¼ 13:82n
ð23Þ
where n is the number of carbon atoms. LeeKesler correlations22,23 are used to express critical temperature, Tc, critical pressure, Pc, and acentric factor ω. Tuning of C50+ Fraction. In a first step, properties of C50+ are calculated by using the previous correlations. Their range of predictibility does not allow property predictions for components heavier than C50. A common way is to tune properties of the last fraction in order to match the phase envelope by using flash calculation.24 Some saturation pressure data are available from differential liberation experiments. Critical properties of the C50 fraction are tuned to match experimental data. Flash simulations are performed in order to get the optimal values. Comparisons between experimental data and flash data for tuned and nontuned properties are shown in Figure 10. Volume Translation Correction. The PengRobinson EoS is used for volume calculation of the gas and oil phase, although it is not accurate for liquid volume prediction. One can apply a volume correction to avoid this issue. Peneloux et al.25 proposed to correct the volume of a component by a constant translation as shown in following vPEN, i ¼ vPR, i ci
ð22Þ 4893
ð24Þ
dx.doi.org/10.1021/ef200365y |Energy Fuels 2011, 25, 4886–4895
Energy & Fuels
ARTICLE
where, for the ith component, vPEN,i is the volume corrected by the Peneloux translation, vPR,i is the volume predicted by the PengRobinson EoS, and ci is the correction. This correction is not quite accurate for heavy components and depends on the temperature. Many works have been done on this topic.26,27 It is impossible to experimentaly get volume (or density) of each components or cuts. So, as a reference, we propose to use the Rackett equation28 which is known to be very accurate29 for volume prediction at saturation. Rackett equation is
Finally, because the oil used in the experiment is a dead oil (stored at laboratory conditions), the mixture is flashed many times at laboratory conditions and the flashed composition is conserved.
The regression is done for a temperature between 50 and 250 C. Out of this range, the Rackett equation is not valid.
Lumping and Final Composition. The lumping is an essential step in the oil characterization because, obviously, numerical simulation cannot be performed with too many components (52 components in this case). At the same time, the lumping, which defines the step of reducing the number of fractions by lumping them in a given numbers of other fractions, does not have to be too coarse at the expense of the accuracy. Whitson30 suggested to lump fractions according to their carbon number and molecular weight. This approach only works for light oils. This method focuses on light fractions whereas for heavy oil, the description should be more accurate for heavy components, which are predominant. As far as we know, there is not an acceptable approach to lump components for thermodynamics purposes. Also, because of the aim of this work, the lumping should be also acceptable for chemistry purposes. Therefore, because of the lack of theory, we decided to adopt a simple approach. Components are lumped into eight fractions: C1, C2C11, C12C16, C17C21, C22C27, C28C35, C36C49, and C50+. This lumping makes possible to preserve the shape of the distribution shown in Figure 9. Equivalent properties are calculated by linear mixing rules for extensive properties and using Leibovici’s rules for other properties.31 The compositional description of the extra-heavy oil is introduced in Table 3. NC and NH are respectively the number of
Figure 10. Phase envelopes for flash and experimental results.
Figure 11. Evolution of oil saturation and density for full and lumped compositions.
vRA, i ¼
! RTc, i n ZRA , Pc, i
n¼1 þ
T 1 Tc, i
!2=7
ð25Þ
where Z RA is the Rackett coefficient. If the specific gravity is known, It can be estimated by this expression Mw P c ZRA ¼ ð26Þ 999RTc SG By assuming the volume correction is a linear correction as a function of the temperature, ci = αiT + βi, and by assuming the Rackett equation as a reference, αi and βi coefficients are deduced from a linear regression as ci ¼ vPR, i vRA, i ≈ αi T þ βi
ð27Þ
Table 3. Extra Heavy Oil Compositionnal Description zi
coupe
Mw
Tc (K)
Pc (bar)
ω
α
β 7
NC
NH
C1
0
215.15
51.99
0.008
0
1.5402 10
C2C11
0.0195
143.6
635
24
0.464
2.6559 108
9.6624 106
10.5422187
17.1516878
C12C16
0.181
193.8
701
19
0.608
8.614 108
4.2135 105
14.1865432
23.5837245
7.9091 105
19.0034845
35.3584878 44.4786478 56.7048882
16
7
1
C17C21
0.154
263.4
772
15
0.787
1.2707 10
C22C27 C28C35
0.135 0.115
336.2 430.4
826 879
12 9
0.946 1.104
1.6002 107 1.9655 107
0.00012037 0.0001755
24.3176319 31.1478656
C36C49
0.101
573.0
936
7
1.272
2.463 107
0.00025697
41.4653185
C50+
0.291
1033.9
1260
6
1.65
1.6075 108
0.000465
74.8
4894
4
75.4668797 135.4
dx.doi.org/10.1021/ef200365y |Energy Fuels 2011, 25, 4886–4895
Energy & Fuels carbon and hydrogen atoms. This data is useful to define stoichiometry coefficients for chemical reactions. Finally, the lumping stage is validated by a comparison using the kinetic cell simulator for the full description and the lumped description. Experimental geometry conditions are respected, but no reactions are input. Only the evaporation effects are observed. This is an extreme situation where only phase equilibria and mass transport are coupled. Results are plotted in Figure 11. The lumped composition gives similar results as the full composition. This lumped composition is considered as acceptable. Model results should not be sensitive to this lumping stage.
’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected]. Present Addresses ‡
TOTAL, avenue Larribau, 64018 Pau, France.
’ ACKNOWLEDGMENT The authors would like to thanks TOTAL SA for its financial support and Schlumberger for the postdoctoral fellowship of A.L. ’ REFERENCES
ARTICLE
(18) Cao, H. Development of techniques for general purpose simulators, Ph.D. thesis, Stanford University, 2002. (19) Whitson, C.; Michelsen, M. Fluid Phase Equilib. 1989, 53, 51–71. (20) Lapene, A.; Nichita, D.; Debenest, G.; Quintard, M. Fluid Phase Equilib. 2010, 297 (1), 121–128. (21) Wilson, G., A modified Redlich-Kwong equation of state applicable to general physical data calculations, Paper No15C, 65th AIChE National meeting, May, 1968. (22) Lee, B.; Kesler, M. AIChE J. 1975, 21, 510–527. (23) Kesler, M.; Lee, B. Hydrocarbon Process. 1976, 55, 153–158. (24) Zurita, R. A.; McCain, W. An Efficient Tuning Strategy to Calibrate Cubic EOS for Compositional Simulation; SPE Annual Technical Conference and Exhibition, San Antonio, TX, Sept 29Oct 2, 2002 (25) Peneloux, A.; Rauzy, E.; Freze, R. Fluid Phase Equilib. 1982, 8, 7–23. (26) Soreide, I. Improved Phase Behavior Predictions of Petroleum Reservoir Fluids from a Cubic Equation of State, Doctor of Engineering Dissertation, Norwegian Institute of Technology: Trondheim, Norway, 1989. (27) Magoulas, K.; Tassios, D. Fluid Phase Equilib. 1990, 56, 119–140. (28) Rackett, H. J. Chem. Eng. Data 1970, 15 (4), 514–517. (29) Riazi, M. Characterization and Properties of Petroleum Fractions; ASTM, 2005. (30) Whitson, C. SPE J. 1983, 683–694. (31) Leibovici, C. Fluid Phase Equilib. 1993, 87, 189–197.
(1) Lake, L. Enhanced Oil Recovery; Prentice Hall, 1989. (2) Burger, J.; Sourieau, P.; Combarnous, M. In Recuperation Assistee du Petrole: Les Methodes Thermiques; Publications de l’Institut Franc-ais du Petrole, 1984. (3) Prasad, R.; Slater, J. High Pressure Combustion Tube Tests. 5th SPE/DOE Symposium on Enhanced Oil Recovery, Tulsa, Oklahoma, April 2023, 1986. (4) Mamora, D.; Ramey, H.; Brigham, W.; Castanier, L. Kinetics of In Situ Combustion; 1993. (5) Sarathi, P. In situ combustion handbook - principles and practices; United States Department of Energy: 1999. (6) Lapene, A.; Castanier, L.; Debenest, G.; Quintard, M.; Kamp, A.; Corre, B. Effects of Water on Kinetics of Wet In-Situ Combustion. 2009 SPE Western Regional Meeting, San Jose, CA, March 2426, 2009. (7) Weijdema, J. Studies on the Oxidation Kinetics of Liquid Hydrocarbons in Porous Media With Regard to Subterranean Combustion. Erdol Kohle. Erdgas Pet. 1968, No. Sept, 520–526. (8) Fassihi, M.; Brigham, W.; Ramey, H. Reaction Kinetics of In Situ Combustion-Part 1 Observations. SPE California Regional Meeting, SPE 8907, August 1984. (9) Fassihi, M.; Brigham, W.; Ramey, H. SPE J. 1984 (Aug); SPE 9454. (10) Cinar, M.; Castanier, L.; Kovscek, A. Improved Analysis of the Kinetics of Crude-Oil In-Situ Combustion. SPE Western Regional and Pacific Section AAPG Joint Meeting, Bakersfield, California, USA, March 29April 2, 2008. (11) Lapene, A.; Castanier, L.; Debenest, G.; Quintard, M.; Kamp, A.; Corre, B. SPE Reserv. Eval. Eng. 2009, 12 (4), 508–517. (12) Kristensen, M.; Gerritsen, M.; Thomsen, P.; Michelsen, M.; Stenby, E. Transport Porous Media 2007, 693, 383–409. (13) Michelsen, M. Fluid Phase Equilib. 1982, 9, 1–19. (14) Michelsen, M. Fluid Phase Equilib. 1982, 9, 21–40. (15) Lapene, A.; Debenest, G.; Quintard, M.; Castanier, L.; Kovscek, A.; Gerritsen, M. Energy Fuels, to be submitted. (16) Peng, D. Y.; Robinson, D. B. Ind. Eng. Chem. Fundam. 1976, 15, 59–64. (17) Prautniz, J.; Lichtenthaler, R.; Azevedo, E. Molecular Thermodynamics of Fluid-Phase Equilibria (3rd Edition). Prausnitz, jm/ Lichtenthale , Prentice Hall , 1998-10-22; Wiley, 1964. 4895
dx.doi.org/10.1021/ef200365y |Energy Fuels 2011, 25, 4886–4895
