Multicomponent effects in liquid–gas filtration combustion

Combustion and Flame 169 (2016) 51–62

Contents lists available at ScienceDirect

Combustion and Flame journal homepage: www.elsevier.com/locate/combustflame

Multicomponent effects in liquid–gas filtration combustion M.A. Endo Kokubun a,∗, N. Khoshnevis Gargar b, H. Bruinning b, A.A. Mailybaev a a b

Instituto Nacional de Matemática Pura e Aplicada – IMPA, Rio de Janeiro, Brazil TU Delft, Civil Engineering and Geosciences, The Netherlands

a r t i c l e

i n f o

Article history: Received 1 December 2015 Revised 5 April 2016 Accepted 6 April 2016

Keywords: Filtration combustion In situ combustion Medium temperature oxidation Multicomponent effects

a b s t r a c t This paper develops the theory of liquid–gas filtration combustion, when an oxidizer (air) is injected into porous rock containing two-component liquid fuel. We found a qualitatively new combustion mechanism controlled by the successive vaporization and condensation of the liquid phase sustained by the reaction. Motivated by the problem of recovery of light oil by air injection, as an enhanced oil recovery method, we consider a liquid composed of light and medium pseudo-components. The light part is allowed to oxidize and vaporize, while the medium part is non-volatile and only oxidizes. The liquid mobility depends strongly on its composition, with a small viscosity (high mobility) of the purely light component and a high viscosity for the purely medium (immobile) component. We show that the combined vaporization and condensation in the combustion wave leads to accumulation of the light component in the upstream part of the wave, considerably increasing mobility and, therefore, playing a crucial role in the mechanism of the combustion process. We describe physical implications of this effect, as well as its importance for applications. The results are confirmed by numerical simulations. © 2016 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

1. Introduction Filtration combustion, when an oxidizer is injected into a porous medium containing fuel, has numerous applications in technology and nature (enhanced oil recovery by in situ combustion, coal gasification, self-propagating high-temperature synthesis, smolder waves etc.). A large number of theoretical studies in this area consider solid (immobile) fuels [1,2]. Then the combustion process can be described using a single-phase flow model for an oxidizer. Examples of such studies are forced forward smoldering [3], downward buoyant filtration combustion [4] and the advance of in-situ combustion fronts in porous media [5]. These models can be extended to study the effects of gas–solid non-equilibrium in filtration combustion [6] and the transition from smoldering to flaming regimes [7], for example. When the fuel is a liquid with a relatively low viscosity and boiling temperature, the two-phase flow model becomes essential. In this case, combustion process couples both to the multi-phase flow and to the phase transition mechanism, increasing the complexity of the problem. Recent theoretical [8,9], numerical [10] and experimental [11] advances in this problem showed the fundamental difference of the combus-

∗

Corresponding author. E-mail addresses: [email protected], [email protected] (M.A. Endo Kokubun), [email protected] (N. Khoshnevis Gargar), j.bruining@ tudelft.nl (H. Bruinning), [email protected] (A.A. Mailybaev).

tion process with liquid fuels, compared to the classical filtration combustion for solid fuels. In this paper, our goal is to deepen the understanding of this process by considering multicomponent liquid fuels. We show that the presence of multiple components in a liquid fuel with different physical properties has a dramatic effect on the internal structure of the combustion wave, manifesting itself in a highly counterintuitive way. The model considered in this work is motivated by the petroleum engineering applications, where in situ combustion is considered as a technique for enhancing the recovery rate due to a lowering of the oil viscosity thus increasing its mobility [12–14]. This technique is usually applied for heavy oils. However, due to thermal expansion and gas drive promoted by the oxidation reaction and vaporization, it can improve the recovery of light oils [15–17]. When air is injected into the reservoir at medium pressures (10–90 bar), the oxidation mechanism is fundamentally different for light and heavy oils. The heavy oil undergoes cracking, due to the increase in temperature, which generates coke deposited on the rock surface. This coke reacts with the injected oxygen, leading to the formation of a high-temperature oxidation (HTO) wave that propagates in the reservoir [5,18]. On the other hand, the oxidation of light oil leads to the scission of liquid molecules, generating a gaseous product. In this case, the temperature is bounded by the boiling point of the oleic phase, such that it does not get large. This is termed the medium-temperature oxidation (MTO) mechanism and recent experiments support the

http://dx.doi.org/10.1016/j.combustflame.2016.04.008 0010-2180/© 2016 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

52

M.A. Endo Kokubun et al. / Combustion and Flame 169 (2016) 51–62

Nomenclature Arl , Arm cg Cm , Co fo , fg ko , kg n P Patm Qrl , Qrm Qv R so , sg t T Tb Tbn Tini Tlac , Tmac uo , ug , u ugl , ugk , ugr uol , uom uinj Wv Wrl , Wrm Xl , Xm ini Xm x Yl , Yk , Yr in j Yk

λ

μo , μ g μl , μ m ν ol , ν gl , ν om , ν gm ρo, ρg ρl, ρm ϕ

frequency factors for light (l) and medium (m) component reactions (1/s) heat capacity of gas (J/mol K) heat capacities of rock and oil (J/m3 K) oil and gas fractional flow functions oil and gas relative permeabilities (m2 ) reaction order prevailing gas pressure (Pa) atmospheric pressure (Pa) enthalpies of light and medium component reactions per mole of oxygen (J/mol) light component vaporization heat (J/mol) ideal gas constant (J/mol K) oil and gas saturations time (s) temperature (K) light component boiling temperature (K) boiling temperature at Patm (K) initial temperature (K) activation temperatures for reactions of light and medium components (K) oil, gas and total Darcy velocities (m/s) Darcy velocities of light (l), oxygen (k) and remaining (r) gas components (m/s) Darcy velocities of light and medium component (m/s) air injection Darcy velocity (m/s) vaporization rate (mol/m3 s) reaction rates of light and medium components (mol/m3 s) volumetric fractions in oleic phase medium component fraction in initial oil spatial coordinate (m) molar fractions in gas phase (mol/mol) oxygen fraction in injected gas thermal conductivity of porous rock (W/m K) viscosities of oil and gas (Pa s) viscosities of purely light and purely medium oleic phases (Pa s) stoichiometric coefficients for light and medium component reactions molar densities of oil and gas (mol/m3 ) molar densities of purely light and purely medium oleic phases (mol/m3 ) rock porosity

existence of reaction occurring at low temperatures [11]. Pascual et al. [19] performed a high pressure tube test using light crude oil to simulate the LTO process, with a stable reaction front at a temperature of 250 °C. They showed that reservoir oil had excellent burning characteristics, which made the process technically feasible. When thermal losses through the rock are high, or if the heat released by the reaction is not sufficient to increase the temperature significantly, the oxidation reaction occurs not far from the initial reservoir temperature [20]. In this case, oxidation is slow and can be incomplete, and oxygen consumption occurs in a larger reservoir zone [21]. The MTO process for light oil requires the two-phase (liquid– gas) flow model. The combustion (MTO) wave was analyzed theoretically [8,9] for a simplified two-phase model with one pseudo-component oil (single component liquid fuel) and the

solution to this problem was obtained as a series of traveling waves. Inside the combustion wave the reaction region separates the vaporization zone upstream and the condensation zone downstream. Note that the position of the vaporization zone in the upstream part of the combustion wave is opposite to the HTO process, where it travels on the downstream side [18]. The speed of the wave was shown to be equal to the Buckley–Leverett characteristic speed evaluated at a point separating the vaporization and condensation regions (called a resonance condition), and this condition allowed determining all macroscopic parameters of the combustion wave. The MTO process within a two-component oil model was studied in [22] by numerical simulations, where an oleic mixture composed of light (low viscosity) and medium (high viscosity) fractions was considered. The light component was allowed to vaporize and oxidize, while the medium non-volatile component only oxidized. Numerical simulations showed that the light component has a tendency to accumulate in the upstream part of the combustion wave despite of vaporization and the much higher mobility, indicating the importance of multicomponent considerations. However, for mixtures composed predominately of a medium non-volatile fraction, the transient behavior featuring the HTO process was reported. In the present work we provide the analytical theory for liquid–gas filtration combustion with a two-component oil as a liquid fuel. In our formulation, following the model of [22], we consider a two-component oil mixture, consisting of light and medium pseudo-components, where the medium component is non-volatile and immobile. The latter property is represented by the liquid mobility depending on the composition, which is high for a purely light oil and vanishes for purely medium oil. Such model is the simplest possible that considers multicomponent effects and it mimics a wide variety of oil types from very light to very heavy ones. The result of this paper is the detailed study of the combustion wave profile providing analytic formulas describing its macroscopic properties: limiting states and speed. We show that, due to the phase transition mechanism, the medium component is almost completely expelled from the reaction region to the downstream side of the wave and, thus, it does not react. This important property, which is valid for any initial oil composition cannot be captured with one pseudo-component models, e.g., [8,9]. Furthermore, the described effect is very counterintuitive: the medium component decreases the oil mobility and, thus, would be expected to accumulate on the opposite upstream side. The consequences are dramatic: the combustion wave speed turns out to be independent on the initial oil composition even if the medium component dominates in the initial oil. For such heavy initial mixtures, the coupled reaction/flow/phase-transition mechanism in the combustion wave yields a strong increase of gas drive and provides high downstream oil saturations, both favorable for oil recovery applications. Finally, we perform numerical simulations supporting our theoretical conclusions. 2. Two pseudo-component oil model We study a combustion front in two-phase flow in which a gaseous oxidizer (air) is injected into a porous rock filled with oil composed of a light and a medium fraction. The light oil (volatile component) can both vaporize and oxidize, whereas the medium oil (non-volatile component) can only oxidize. In our applications we disregard gaseous phase reactions, as annihilation of free radicals at pore walls reduces drastically the corresponding reaction rates [23]. We summarize the reaction process of liquid components in the following reaction equations:

νol (light component) + O2 → νgl (gaseous products), νom (medium component) + O2 → νgm (gaseous products),

(2.1)

M.A. Endo Kokubun et al. / Combustion and Flame 169 (2016) 51–62

i.e., one mole of oxygen reacts with ν oi moles of oleic (liquid) hydrocarbons (where i = l, m denotes light and medium components, respectively) generating ν gi moles of gaseous products (H2 O, CO2 , etc.). This reaction mechanism is used in order to consider the medium temperature oxidation, when oxygen reacts with oil, leading to a scission of liquid molecules and generating gaseous products. Such mechanism is fundamentally different from the low-temperature oxidation (LTO), when the oxygen molecule is incorporated into the liquid oil, generating another liquid compound, and the high-temperature oxidation, where the oxygen reacts with the coke generated by the cracking due to high temperatures. Water originally present or condensed from the reaction products is neglected. Effect of water in the one-component case was studied numerically in [10] and it was shown that although recovery rates may change, the overall wave structure does not change. The changes in our solution due to the presence of water will be discussed in future works. We study two-phase (liquid and gas) one-dimensional flow in the positive spatial direction x. The oil has saturation so , describing the oil occupied fraction of pore volume. The saturation of gas is, therefore, equal to sg = 1 − so. In the gaseous phase, we distinguish the molar fraction of the light hydrocarbon Yl and the molar fraction of oxygen Yk . The medium oil component does not exist in the gaseous phase. The remaining components with fraction Yr = 1 − Yk − Yl consist of reaction products and inert components of the injected gas. For the ideal gas, the molar fractions Yl , Yk and Yr are equal to the volume fractions. We assume that the two-component oil has neither heat nor volume effects due to mixing. Let Xm and Xl be the volume fractions of the medium and light components (Xm + Xl = 1), and ρ m and ρ l be the pure oleic molar densities of the medium and light components, respectively. We remark that the model in a previous paper [22] was based on molar fractions. The molar mass balance equations for liquid and gas components can be written as [22]

∂ ∂ ϕ X ρ s + ρ u = −νomWrm , ∂ t m m o ∂ x m om

(2.2)

53

It is convenient to express the liquid and gas velocities as

uo = u f o ,

ug = u f g ,

(2.9)

where u is the total Darcy velocity with the oil and gas fractional flow functions

fo (so, T , Xm ) =

ko/μo , ko/μo + kg /μg

fg = 1 − fo.

(2.10)

Dependence of the relative permeability functions on their respective saturations is chosen as

ko = k

(so − sor )2 (1 − sor )2

if

so ≥ sor ;

ko = 0

if

so ≤ sor ;

kg = ks2g , (2.11)

where sor is the residual oil saturation and k is the rock permeability and has unit of m2 . The oleic viscosity μo depends on the pure light and medium oil viscosities, μl and μm , and we use the expressions



μo =

−b

Xl

μ1l /b

+

Xm



,

μ1m/b

μ j = μ j0

Ej exp − RT

 ,

(2.12)

with b = 4 [24]. Here Ej is the activation energy for the viscosity of the j-component in the oleic phase and μj0 is a reference viscosity. For Xl = 1, one has μo = μl , while for Xm = 1 one has μo = μm . For simplicity, we assume that the pure medium component is very viscid (immobile), i.e, we consider the limit μm  μl , such that the oleic phase viscosity is modeled as

μo = μl Xl−b .

(2.13)

The gas (air) viscosity is given by [25]

μg =



7.5 T T + 120 291

3/2

Pa s.

(2.14)

Taking the sum of Eqs. (2.4)–(2.6) and using (2.8) and (2.9) with Yl + Yk + Yr = 1 yields the equality ugl + ugk + ugr = ug as well as the balance law for the total gas as

∂ ∂ ϕρ s + ρ f u = (νgm − 1 )Wrm + (νgl − 1 )Wrl + Wv . ∂t g g ∂ x g g

(2.15)

(2.16)

∂ ∂ ϕ X ρ s + ρ u = −νol Wrl − Wv , ∂ t l l o ∂ x l ol

(2.3)

∂ ∂ ϕY ρ s + ρ u = Wv , ∂ t l g g ∂ x g gl

Similarly, taking the sum of Eq. (2.2) (divided by the constant ρ m ) and (2.3) (divided by the constant ρ l ) and using (2.7) and (2.9), yields the balance law for the total liquid as

(2.4)

1 ∂ ∂ νom ν ϕs + f u=− W − ol W − W . ∂t o ∂ x o ρm rm ρl rl ρl v

∂ ∂ ϕY ρ s + ρ u = −Wrl − Wrm , ∂ t k g g ∂ x g gk

(2.5)

Assuming that the temperature of solid rock, liquid and gas are equal, we write the heat balance equation as

∂ ∂ ϕY ρ s + ρ u = νgl Wrl + νgmWrm , ∂ t r g g ∂ x g gr

(2.6)

∂  eff ∂ Cm + ϕCoso + ϕ cg ρg sg T + (C u + c ρ u )T ∂t ∂x o o g g g ∂ 2T = λ 2 + QrmWrm + Qrl Wrl − QvWv , ∂x

where uα j means the Darcy velocity of component j in the phase α . In summary, there are four components, viz., light (l), medium (m), oxygen (k), and the rest (r). The medium component can only exist in the oleic phase (o), the light component can exist in both the oleic (o) and the gaseous phase (g), whereas oxygen and the remaining gases can only exist in the gaseous phase. The expressions for velocities of medium and light components in the oleic phase can be derived as (neglecting mass diffusion, which is a small effect in this problem and will be discussed more ahead)

uom = Xm uo,

uol = Xl uo.

(2.7)

Similarly, we approximate the velocities of components in the gas phase by

ugl = Yl ug ,

ugk = Yk ug ,

ugr = Yr ug .

(2.8)

(2.17)

where T = T − Tini with the initial reservoir temperature Tini and eff Cm ≡ (1 − ϕ )Cm is the rock effective heat capacity. In Eq. (2.17), Qrm and Qrl are reactive enthalpies and Qv is the vaporization heat. The heat capacities Cm , Co , cg and the effective thermal conductivity λ are taken as constants, which is a reasonable approximation and facilitates the analysis. We neglect heat losses, which are usually very small in field applications (however, taking into account heat losses becomes essential for interpreting laboratory experiments). We use the ideal gas law to define

ρg = P/RT .

(2.18)

Pressure variations are assumed to be small compared to the prevailing pressure, so we take P ≈ const in (2.18), as well as in

54

M.A. Endo Kokubun et al. / Combustion and Flame 169 (2016) 51–62

the relationships (2.19) and (2.20) below. In other words, physical properties of the phases are evaluated at constant pressure. eq The partial pressure Yl P of the gaseous hydrocarbon in liquid– gas equilibrium is derived using the Clausius–Clapeyron relation with Raoult’s law as (see [25, p. 8.12])

 Q 1 v

Yleq P = xl γl Patm exp −

R

−

T



1 Tbn

,

(2.19)

where Tbn is the (normal) boiling point of light hydrocarbon measured at atmospheric pressure Patm and xl is the molar fraction of the light component. We choose the activity coefficient γ l such that xl γl = Xl , representing the volume fraction, which is considered to take into account part of the non-ideality of Raoult’s law. The purpose of this is that it simplifies the analysis; with some effort it is also possible to explicitly take into account the activity coefficient (see [25, Tables 8–3]). Relation (2.19) determines the eq equilibrium fraction of gaseous hydrocarbon Yl (recalling that only the light component of the oil vaporizes). The boiling temperature T = Tb of light hydrocarbon at pressure P is obtained from eq Eq. (2.19) by setting Yl = 1 and Xl = 1. We express the oxidation (MTO) reaction rates for light and medium components by the Arrhenius form [26,27]

Wrl = Wrm =

ϕ Arl Xl ρl so

 PY n k

Patm

ϕ Arm Xm ρm so



exp −

 PY n k

Patm

Tl



ac



,

T

exp −

Tmac T

 ,



Wv = kvl ϕ Yleq − Yl ρg soXl ,

(2.21)

which is proportional to the deviation of the light component fraceq tion Yl in the gas phase from its equilibrium value Yl with an empirical transfer parameter kvl . This formulation can be considered a consequence of non-equilibrium thermodynamics, see for instance [23,28]. We assume that kvl is large, describing the situation close to local thermodynamic equilibrium for the gaseous light fraction Yl , i.e., very fast vaporization or condensation. The vaporization rate Wv vanishes under the conditions

or

ρg∗ ρg∗ Qv , βm = , γv = , Cm Cm ρl ρm Qrl Q T Qv ϕv∗ γr = rm , θ0 = ini∗ , θh = , σ = . (3.3) ∗ Qrl T R T uin j

αo =

ϕCo

αg =

,

ϕ cg ρg∗

soXl = 0.

(2.22)

We note that the physical relations for relative permeability (2.11 ), viscosities (2.13), (2.14), and source terms (2.20), (2.21) are used as examples, since the developed theory does not rely on their specific form.

∂ ∂ (1 + αoso + αg Sg )θ + (αo fo + αg Fg )uθ ∂t ∂x ∂ 2θ σ = + (γ w + wrl − γv wv ), ∂ x2 Y in j r rm ∂ so ∂ u f o + = −νom βm wrm − νol βl wrl − βl wv , ∂t ∂x

(3.5)

∂ Sg ∂ uFg + = (νgm − 1 )wrm + (νgl − 1 )wrl + wv , ∂t ∂x

(3.6)

∂ Yl Sg ∂ uFgYl + = wv , ∂t ∂x

(3.7)

∂ Yk Sg ∂ uFgYk + = −wrm − wrl , ∂t ∂x

(3.8)

∂ Xm so ∂ Xm u fo + = −νom βm wrm , ∂t ∂x

(3.9)

where the temperature-corrected gas saturation and flux are defined as

S g ( s o , θ ) = ( 1 − s o ) ρg , 1

ρg =

1 + θ /θ0

(3.10)

The dimensionless forms of the reaction and vaporization rates (2.20) and (2.21) are given by



wrm = arm Xm soYkn exp −

θmac =

θ=

T − Tini , T ∗

u˜ =

u

ϕv∗

,

ρ˜g =

ρg , ρg∗

(3.1)

(3.11)



ρg∗ =

x∗

v∗

,

P , RTini

x∗ =

λ Cm v∗

,

v∗ =

T ∗ = Tb − Tini .

Tlac

T ∗

 θlac , arl = t ∗ Arl βl−1 (P/Patm )n , θ + θ0

,

(3.12)

wv = κvl (Yleq (θ , Xl ) − Yl )ρg Xl so,

κvl = kvl t ∗ .

(3.13)

For the fraction of light hydrocarbon in the gaseous phase in equilibrium with liquid, we have

where the reference quantities denoted by an asterisk are

t∗ =

 θmac −1 , arm = t ∗ Arm βm (P/Patm )n , θ + θ0

Tmac , T ∗

θlac =

x , x∗

Fg (so, θ , Xm ) = (1 − fo (so, θ , Xm ))ρg ,

.

In order to render the governing equations dimensionless, we introduce the ratios

x˜ =

(3.4)

k

3. Dimensionless equations

t , t∗

βl =

The dimensionless governing equations obtained from Eqs. (2.17), (2.16), (2.15), (2.4), (2.5) and (2.2) are given by (we drop the tildes for simplicity)

wrl = arl Xl soYkn exp −

t˜ =

,

(2.20)

where Arl , Arm are the frequency factors and Tlac , Tmac are the activation temperatures. We use volume fractions Xl and Xm instead of mole fractions again to take into account part of the non-ideality of the liquid. The vaporization rate is given by

s0 Xl > 0, Yl = Yleq

the rock temperature to Tb . Dimensionless parameters are introduced as

ρ

Qrl g∗ uin jYkin j , Cm T ∗

Yleq (θ , Xm ) = (1 − Xm )Y (θ ),



(3.2)

The dimensionless variable θ describes the temperature distribution with the reservoir condition θ = 0 and the boiling temperature θ = 1 for purely light oil. The reference quantities t∗ , x∗ and v∗ are obtained by considering a wave whose combustion heat raises

(3.14)

where we defined the function Y (θ ) as



Y (θ ) = exp

θh

θ0 + 1

−

θh

θ0 + θ



.

(3.15)

The following properties of the dimensionless parameters will be multiply used in our analysis

M.A. Endo Kokubun et al. / Combustion and Flame 169 (2016) 51–62

55

Fig. 1. Left: schematic of the physical problem under consideration; Right: wave sequence solution with the thermal, combustion, and saturation waves traveling with different speeds, vT < v < vS . Shown are the schematic profiles of temperature θ , oxygen fraction Yk and oil saturation so in the rock cylinder.

σ 1, βl 1, Y (0 ) 1,

βl Cm T ∗ =  1, σ ϕ Qrl ρl Ykin j

(3.16)

where σ is the ratio between the combustion wave and the air injection velocities and β l is the ratio between the gas and the light oil densities. The last approximation is given by realizing that β l /σ represents the ratio between the heat accumulated in the rock and heat released by the reaction of light oil. Since the heat of vaporization is small, in comparison with the heat released by the reaction, most of the heat is used to raise the rock temperature. The condition Y (0 ) 1 implies a small fraction of gaseous light comeq ini ) 1, and it is a ponent at the reservoir temperature, Yl (0, Xm consequence of large θ h in Eq. (3.15). The system (3.4)–(3.9) contains six equations for the variables θ , so , u, Yl , Yk and Xm dependent on t and x . The initial (reservoir) conditions are given by

θ = 0, so = sini o , ini Yl = (1 − Xm )Y (0 ), Yk = 0, t = 0, x ≥ 0 :

5. Combustion wave

Xm = Xmini , (3.17)

corresponding to the reservoir containing oil with a specified ini . The injection conditions are given by medium fraction Xm

x = 0, t ≥ 0 : Yk = Ykin j ,

αg Fg uθ −

u = σ −1 ,

∂θ = 0, ∂x

so = Yl = 0, (3.18)

corresponding to the injection of an oxidizer (air) at reservoir temperature and the dimensionless injection speed σ −1 given by its dimensional value uinj in Eq. (3.3). It is assumed that there is no gaseous light hydrocarbons in the injected gas, Yl = 0. 4. Wave solutions For large times, the solution of the problem contains three distinct waves: thermal, combustion and saturation waves, as shown in Fig. 1. The thermal wave is slower and it is located upstream of the combustion wave because of the high thermal capacity of the rock. In the thermal wave region there is no volatile liquid components, and no reaction occurs (then the oxidant fraction is constant at its injection value). From the analysis of the thermal wave (see e.g., [9]), one derives the wave speed and the gas flux as

vT =

αg u , uFg = ≈ σ −1 . σ 1 + θ /θ0

matching condition for the combustion wave. The saturation wave is the fastest wave and is located downstream of the combustion wave and resembles a classical Buckley–Leverett solution. In this wave, the oil saturation so is changing, while the other variables remain constant. The combustion (MTO) wave, where oxidation and phase transitions occur, is the focus of this work. The gaseous hydrocarbons, which are generated by vaporization of the light component, flow downstream of the combustion wave towards the region of lower temperatures, where they condense. Hence, the combustion wave is composed of the vaporization/condensation front, which is sustained by the heat released in the oxidation reaction. This wave is analyzed in detail in the next section. We will show that the processes in the combustion wave are dominated by the light oil component, even when the oil is predominately non-volatile.

(4.1)

The thermal wave speed is small, vT 1, because it is proportional to the ratio of the air to rock heat capacities, see Eq. (3.3). The second expression in Eq. (4.1) describes the change of the gas velocity due to thermal expansion and it will be used as a

Our analysis will show that almost no medium oil component is present in the reaction region of the combustion wave. This surprising result is a consequence of the vaporization/condensation process. Therefore, we facilitate the derivations by neglecting the medium oil reaction, i.e., setting wrm = 0, which will be confirmed by the final result. For further simplification, we consider νgl = 1 (no net gas production by the reaction) and ν ol 1, which means that the amount of light oil consumed by reaction is small. We also consider that γ v 1, meaning that the heat of vaporization is much smaller than the heat of combustion. Thus, as a first approximation we set

νol = 0, γv = 0.

(5.1)

Validity of these simplifications is also confirmed by numerical simulations described below. We assume that the combustion wave has a stationary profile traveling at a constant speed v. In order to analyze its structure, we change to a moving reference frame with ξ = x − vt. In this new coordinate, the governing Eqs. (3.4)–(3.9) for the stationary profile take the form

d d2 θ σ (−v + αoψo + αg ψg )θ = 2 + in j wrl , dξ dξ Yk

(5.2)

d ψo = −βl wv dξ

(5.3)

d ψg = wv , dξ

(5.4)

56

M.A. Endo Kokubun et al. / Combustion and Flame 169 (2016) 51–62

where the light component fraction Yl is given by the equilibrium value (3.14), while the oil saturation sdo and gas speed ud are unknown. Therefore, the unknowns are the limiting states θ u , suo , sdo , ud and the wave velocity v. 5.1. Limiting states This Section derives the relations for the limiting states of the combustion wave that follow from the overall balance conditions. Manipulating Eqs. (5.2)–(5.7) in order to eliminate the source/sink terms, we obtain the following equations for the conserved scalars

d dξ Fig. 2. Combustion wave profile in the moving reference frame: air injected on the left side and the wave moves to the right.



(−v + αoψo + αg ψg )θ −

σ ψgYk dθ + dξ Ykin j



= 0,

(5.12)

d (ψo + βl ψg ) = 0, dξ

(5.13)

d (ψoXm + βl ψg (1 − Yl ) ) = 0, dξ

(5.14) (5.15)

dψgYl = wv , dξ

(5.5)

dψoXm = 0. dξ

dψgYk = −wrl , dξ

(5.6)

dψoXm = 0, dξ

We added Eq. (5.15) into Eq. (5.14) because it will simplify the analysis. Integrating Eqs. (5.12)–(5.15) from upstream to downstream using (5.9)–(5.11), we obtain

(5.7)

where we used Eq. (5.1) and defined

ψo = u f o − vso ,

ψg = uFg − vSg .

(5.8)

The quantities ψ o and ψ g represent, respectively, the oil and gas fluxes in the moving reference frame. All the dependent variables are now functions of ξ only. We expect, and it will be confirmed by the final result, that the wave speed v ∼ 1 due to the choice of the reference value v∗ in Eq. (3.2). In Fig. 2 we show a schematic profile of the combustion wave. Upstream of the wave the temperature is θ u (unknown), all light oil is vaporized (which is equivalent to Xm = 1) and there is no gaseous light oil in the injected gas, Yl = 0. Downstream of the wave the oil is at its initial temperature θ = 0 and composition ini . Since we neglected the reaction of the medium comXm = Xm ponent, the existence of a (pure medium) oil saturation suo at the upstream side is assumed. As we mentioned, the result will show that suo is negligibly small. Therefore, we have the upstream boundary conditions

ξ → −∞ :

θ = θ u, u=

Yl = 0,

1 + θ u /θ0

σ

Yk = Ykin j ,

Xm = 1,

,

so = su0 , (5.9)

where the last equality follows from the match with the thermal wave conditions (4.1). The upstream values of the fluxes ψ o and ψ g given by Eqs. (5.8) and (3.10) are determined by realizing that fo = 0 and fg = 1 for pure medium oil, see Eqs. (2.10)–(2.13), such that

ξ → −∞ :

ψo = −v

suo ,

1 v(1 − suo ) ψg = − ≈ , σ 1 + θ u /θ0 σ 1

(5.10)

where in the last approximation we used (3.16) and the fact that v ∼ 1. We assume complete oxygen consumption in the combustion wave, leading to the following downstream conditions

ξ → +∞ :

θ = 0, so = sdo ,

Yl = (1 − Xmini )Y (0 ), u = ud ,

Yk = 0,

Xm = Xmini , (5.11)



−v − αovsuo +

αg  u θ + 1 = 0, σ

βl − vsuo = ψod + βl ψgd , σ βl − vsuo = ψod Xmini + βl ψgd (1 − (1 − Xmini )Y (0 )), σ −vsuo = ψod Xmini ,

(5.16) (5.17) (5.18) (5.19)

where the upstream states are at the left-hand side and the downstream states are at the right-hand side in the equations above. In Eq. (5.16) we used that dθ /dξ = 0 for both upstream and downstream states. By manipulating Eqs. (5.17) and (5.18) we can find expressions for ψod and ψgd as

  Y (0 ) β ψod = − l − vsuo , σ 1 − Y (0 ) ψgd =

1

βl



 1 βl − vsuo . σ 1 − Y (0 )

(5.20)

(5.21)

The upstream temperature θ u is obtained from Eq. (5.16) as

θu =

1

(v + αovsuo − αg /σ )

.

(5.22)

Substituting Eq. (5.20) in Eq. (5.19) gives the upstream saturation suo as

suo =

  βl Xmini Y (0 ) . σ v 1 − (1 − Xmini )Y (0 )

(5.23)

Using the conditions (3.16) and the consideration v ∼ 1, yields

suo ≈

βl ini X Y ( 0 ) 1, σv m

(5.24)

such that we can neglect the term vsuo in Eqs. (5.20)–(5.22). Then, the upstream temperature θ u is given by

θu =

1

(v − αg /σ )

,

(5.25)

M.A. Endo Kokubun et al. / Combustion and Flame 169 (2016) 51–62

57

Thus, we distinguish a thin vaporization region upstream and a much wider condensation region (CR) downstream. Most of the reaction occurs at highest temperatures, i.e., in the upstream part of the CR, as we describe below. Eqs. (5.16)–(5.19) were obtained by integrating Eqs. (5.12)–(5.15) from upstream to downstream constant states. In order to obtain the wave profile, we integrate Eqs. (5.12) and (5.15) from downstream up to some point ξ inside the wave, and Eqs. (5.13) and (5.14) from upstream to the same point ξ , which yields

(−v + αoψo + αg ψg )θ − Fig. 3. Schematic graph for the left- and right-hand sides of Eq. (5.30) as functions of the oil saturation so .

and the downstream mass flux of oil and gas, using Y (0 ) 1 from Eq. (3.16), are given by

ψod = −

βl Y (0 ) 1, σ

ψgd =

1

σ

 1.

(5.26)

Note that the term α g /σ in the denominator of Eq. (5.25) is equal to the thermal wave speed vT , see Eq. (4.1), which was shown to be small. The values of the downstream fluxes ψod and ψgd can also be determined by Eq. (5.8) in terms of the downstream oil flux fod and saturation sdo as

ψod = ud fod − vsdo , ψgd = ud (1 − fod ) − v(1 − sdo ).

(5.27)

The sum of these expressions yields the gas velocity downstream as

u =v+ψ +ψ d

d o

d g.

(5.28)

We substitute Eq. (5.26) in Eq. (5.28) and obtain

ud = v −

1 βl Y (0 ) 1 + ≈ . σ σ σ

(5.29)

Using the first equation in Eq. (5.27) with the neglected small term ψod 1 from Eq. (5.26) and with ud = σ −1 from Eq. (5.29), we obtain the downstream saturation sdo implicitly through

fo (sdo , 0, Xmini ) = σ vsdo ,

(5.30)

where we expressed the downstream fractional flow function f od in the left-hand side. A typical graph of the left- and right-hand sides of Eq. (5.30) as functions of sdo is shown in Fig. 3. The left-hand side is the s-shaped fractional flow function, while the right-hand side is a linear function of sdo with a small positive slope σ v and starting at so = 0. The two roots for Eq. (5.30) are given by sdo = 0 and sdo > 0, where the relevant root is the positive one [8]. Eqs. (5.24), (5.25), (5.29) and (5.30) determine the conditions at the limiting states suo , θ u , ud and sdo , respectively. In order to close the problem we must obtain an expression for the wave speed v. This expression is obtained from the analysis of the internal wave profile, which will be done in the next Section. 5.2. Combustion wave profile We will show that the wave speed v is determined by the resonance condition at the point where the vaporization of the light component is changed by its condensation. Since vaporization is assumed to be the fastest process, the vaporization region (VR) is thin and located in the upstream part of the wave, where the injected air gets into contact with the oleic phase, see Fig. 2. Downstream of the VR, the liquid and gaseous oil are close to thermodynamic equilibrium. As the temperature gradually decreases in the downstream direction, slow condensation occurs.

ψo + βl ψg =

dθ Y + σ ψg ink j = 0, dξ Y

(5.31)

k

βl , σ

(5.32)

ψoXm + βl ψg (1 − Yl ) =

βl , σ

(5.33)

ψoXm = ψod Xmini ,

(5.34)

where the right-hand sides were determined from Eqs. (5.16)– (5.19), and we neglected the small term vsuo according to Eqs. (3.16) and (5.24). Eqs. (5.32) and (5.33) yield

ψo = − ψg =

 βl  Yl , σ 1 − Yl − Xm  1−X  m

1

σ 1 − Yl − Xm

(5.35)

.

(5.36)

Eq. (5.34) with ψ o from Eq. (5.35) and ψod from (5.26) can be solved with respect to Xm as

Xm =

Y (0 )Xmini (1 − Yl ) Y (0 )Xmini + Yl

.

(5.37)

We can use (5.37) in Eqs. (5.35) and (5.36) in order to express ψ o and ψ g in terms of Yl as

ψo = − ψg =

βl (Y (0 )Xmini + Yl ) βl Yl ≈− , σ (1 − Yl ) σ (1 − Yl )

Y (0 )Xmini + 1 ≈ σ (1 − Yl )

1

σ (1 − Yl )

,

(5.38)

(5.39)

where the approximations were made because Y (0 ) 1. From the definitions of ψ o and ψ g , Eqs. (5.8) and (3.10), we express the gas velocity u as

u = v + ψo + ψg (1 + θ /θ0 ).

(5.40)

Eqs. (5.38) and (5.39) with conditions (3.16) allow to express the gas velocity (5.40) in terms of Yl and θ as

(1 + θ /θ0 ) − βl Yl (1 + θ /θ0 ) ≈ . σ (1 − Yl ) σ (1 − Yl )

u=v+

(5.41)

Substituting (5.38) and (5.39) in (5.31), after some manipulations and with the use of θ u given by Eq. (5.25), yields



dθ =− dξ

1

θu

+

 (αoβl − αg )Yl Y θ + in j k . σ (1 − Yl ) Y (1 − Yl )

(5.42)

k

The above expressions are valid at all parts of the wave profile. Now let us discuss the properties specific for each region, Fig. 2. In the VR, dθ /dξ given by Eq. (5.42) is not large. Since this region is very thin, we neglect the temperature variation in the VR, such that the temperature is constant and given by θ u in Eq. (5.25). For the same reason, one can neglect reaction in the VR, assuming the in j oxygen fraction Yk ≈ Yk to be approximately constant. In the VR, the gaseous light hydrocarbon Yl increases from 0 to its equilibeq rium value Yl (θ u , Xm ), while the fraction of medium oil Xm given by Eq. (5.37) decreases from 1 to some value.

58

M.A. Endo Kokubun et al. / Combustion and Flame 169 (2016) 51–62

gency condition



v=



∂ fo

u . ∂ so r

(5.45)

The wave speed v is determined from the so-called resonance condition, which is satisfied at the point where the VR and CR meet. For this purpose, we substitute ψ o from Eq. (5.8) and the gas velocity u given by (5.41) in Eq. (5.38), and obtain

Here the wave speed (left-hand side of Eq. (5.45)) is equal to the Buckley–Leverett characteristic speed (right-hand side of Eq. (5.45)), which we call resonance [9], denoted by “r” . In the following, we show that resonance corresponds to the point where the VR and CR meet and the gaseous oil fraction Yl attains a maximum. In the VR, the light component fraction Yl increases due to the vaporization, while the temperature remains approximately constant as we described earlier. The increase in Yl yields larger values of fo according to Eqs. (2.10)–(2.13) and (5.37) as shown in Fig. 4 (horizontal arrow), while the straight-line moves downwards according to Eq. (5.44) (vertical arrow in Fig. 4). Hence, the coordinate so of the lower intersection in Fig. 4 increases as we move along the wave profile in the downstream direction. This tendency cannot be extended beyond the tangency configuration (crossed dot in Fig. 4), as this would lead to no solution (no intersection). In the wide CR, the temperature changes from θ u to its downstream value θ = 0 according to Eq. (5.42). Since α o β l > α g , as the molar heat capacity of the liquid oil is higher than for the gas, see Eq. (3.3), and Yk ≈ 0 (most of oxygen is consumed at highest temperatures), we see from Eq. (5.42) that dθ /dξ < 0. Thus, the change of the temperature is monotonic, decreasing in the downeq stream direction. The equilibrium light oil fraction Yl = Yl given by Eqs. (3.14) and (3.15) lowers because of the decrease in the temperature θ and increase in Xm given by Eq. (5.37). This decrease in Yl and θ yields smaller values of the left-hand side of Eq. (5.44) according to Eqs. (2.10)–(2.13), while the straight-line moves upwards according to Eq. (5.44). This corresponds to the divergence away from the tangent configuration in Fig. 4. The upper solution so (bold dot) increases as we move along the wave profile downstream until it reaches the limiting constant state sdo , see Figs. 2 and 3. Therefore, our analysis proves that the point, where the VR and the CR meet, must correspond to a tangency configuration in Fig. 4 (crossed dot). Because the gaseous oil fraction Yl decreases on both sides of the resonance point, it attains the maximum given eq r ). Since ∂ X /∂ Y < 0 according to Eq. (5.37), we by Ylr = Yl (θ u , Xm m l conclude that the medium fraction in oleic phase attains the minr at the resonance point. imum Xm As we showed in Eq. (5.43) with Y (0 ) 1, the minimum value r 1. Hence, the medium oil fraction can be neglected in the Xm flow function fo in Eq. (5.44) at the resonance point. As a result, Eq. (5.44) at the resonance point yields

(1 + θ /θ0 ) fo (so, θ , Xm ) = vσ (1 − Yl )so − βYl .

(1 + θ u /θ0 ) fo (sro, θ u , 0 ) = vσ (1 − Y (θ u ))sro − β Y (θ u ).

Fig. 4. Left-hand side (dashed curve) and right-hand side (dashed straight line) of Eq. (5.44) as functions of oil saturation so . The arrows indicate the change of the graphs with increasing Yl along the wave profile. Open dot represents the upstream solution, black dot represents the downstream solution and crossed dot represents solution at the resonance point (tangency of the solid curve and solid line).

In the CR, the gaseous light component is close to thermodyeq namic equilibrium, i.e., Yl ≈ Yl (θ , Xm ) = (1 − Xm )Y (θ ). Using this expression in Eq. (5.37) and solving with respect to Xm yields two roots. The first (irrelevant) root is Xm = 1, and the second root is

Xm = Xmini

Y (0 ) (1 − Y (θ ) ). Y (θ )

(5.43)

We assume that the oxidation reaction occurs in the region of high temperatures with complete oxygen consumption. In fact, the more detailed analysis [29] shows that most of the reaction is confined in a small region near the highest temperature θ u if the air injection rate is not large. For large injection rates, the reaction may spread into a wider part of the CR. Since Y (0 ) Y (θ ) for high temperatures θ ≈ θ u , one has Xm 1 in Eq. (5.43), i.e., a very small amount of medium component is available in the reaction ini . This justifies the iniregion, irrespective of the initial amount Xm tial hypothesis that the medium non-volatile component plays a minor role in the reaction. 5.3. Resonance condition

(5.44)

Eq. (5.44) determines implicitly the saturation so in the wave profile. The left- and right-hand sides of Eq. (5.44) are shown schematically in Fig. 4 in terms of so with the other variables fixed. The left-hand side of Eq. (5.44) is proportional to the fractional flow function, determined by Eqs. (2.10) and (2.11), and it is sshaped. The right-hand side of Eq. (5.44) is a straight line with a small inclination, which crosses the vertical axis at the small negative value −β Yl . The solution of Eq. (5.44) is given by the intersection between the two dashed curves in Fig. 4. Upstream of the VR, the lower intersection point is relevant (open dot) because this region starts with so = suo 1, according to ( 5.24). Downstream of the CR, the upper intersection point is relevant (black dot) because in this region so = sdo > 0, where sdo is not small, see Fig. 3. By continuity, there must exit a tangency represented by the crossed dot in Fig. 4 at some internal point of the wave profile. Since the slopes of the left- and right-hand sides of Eq. (5.44) are equal at this internal point, we use Eq. (5.41) and obtain the tan-

(5.46)

The resonance condition (5.45), with u from Eq. (5.41), is written in a similar way as

v=

1 + θ u /θ0 σ (1 − Y (θ u ))





∂ fo

. ∂ so (sr ,θ u ,0)

(5.47)

o

Substituting θ u from Eq. (5.25), the solutions of the system of two nonlinear equations (5.46) and (5.47) for the saturation sro at the resonance point and the wave speed v can be found numerically. With the wave speed v determined, we proceed to obtain the upstream temperature θ u from Eq. (5.25) and the downstream saturation sdo from Eq. (5.30). This determines all unknowns quantities of the combustion wave limiting states. 5.4. Remarks on medium component oxidation The a priori assumption that the oxidation reaction of the medium component is negligible was justified by showing that

M.A. Endo Kokubun et al. / Combustion and Flame 169 (2016) 51–62

a

b

1

c

8

0.8

59

0.1

0.08

suo (x 10 )

6

0.06

Y max k

sdo

3

0.6

4

0.4

0.04

2

0.2

0.02

0 0

0.2

0.4 0.6 Xini m

0.8

1

0

0 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

ini Xm

0.6

0.8

1

Xini m

Fig. 5. (a) Dependence of the downstream oil saturation sdo in the combustion wave on the initial medium fraction Xmini for the reservoir pressures P = 10( ), 40( ) and 70( ) bar. The crosses show the results obtained from numerical simulations at P = 10 bar. (b) Upstream oil saturation suo (×103 ). (c) Oxygen fraction that reacts with the residual medium oil.

only a small fraction Xm remains in the region of high temperatures in the CR according to Eq. (5.43), and a small amount suo remains upstream according to Eq. (5.24). In this Section we make this argument more precise. The non-volatile and immobile medium oil of saturation suo may react with the oxygen at the upstream side of the combustion wave, where the temperature is elevated, see Fig. 2. In dimensionless variables, the rate of immobile medium oil deposition is given by |ψou | = vsuo , which by the stoichiometric relation (2.1) can react with the amount vsuo /(νom βm ) of oxygen, see Eqs. (3.8) and (3.9). With the upstream gas flux ψgu = σ −1 , as given by Eq. (5.10), and suo from Eq. (5.24), one finds the maximum oxygen fraction that can be consumed by the reaction of this residual medium oil as

Ykmax =

|ψou | σ vsuo β X ini = = l m Y ( 0 ). u νom βm ψg νom βm νom βm

(5.48)

According to Eq. (3.3), the ratio βl /(νom βm ) = ρm /(νom ρl ), where ρ m and ρ l are the pure oleic molar densities of the medium and light components, is not expected to be very large. Thus, the value of Ykmax is dominated by the small equilibrium fraction Y (0 ) for the light oil component at the initial reservoir temperature. In practical applications, especially at elevated pressures (see Fig. 5(c) in the next Section), the value of Ykmax is small compared to the in j Yk ,

injected fraction which, once again, confirms our initial assumption of a small effect of the medium component reaction. The actual form of how the residual medium oil is oxidized depends on the magnitude of the reaction rate wrm in Eq. (3.11). If this rate is not too small at the upstream temperature θ u , all medium oil is consumed in the upstream part of the combustion wave. This leads to a small temperature elevation and a small decrease of the oxygen fraction by the amount of Ykmax . In order to take this effect into account, one can use the corrected value in j in j Yk − Ykmax instead of Yk in the analysis of Section 5. It is possible, however, that the reaction rate wrm is very small, e.g., due to a small oil saturation or a higher activation energy. In this case, we may consider the possibility of the high-temperature

oxidation (HTO) wave formation after a certain transient period. The HTO process is well understood, see e.g., [3,18]. In the case of low fuel saturation suo 1, a considerable temperature increase is required for HTO. This is possible when the HTO wave speed vH is close to the thermal wave speed vT = αg /σ , resulting in the superadiabatic √ heat accumulation [30] . The thermal wave width grows as LT ∼ t [9,30], while the separation between the two waves increases linearly as L ∼ (vH − vT )t. Since we need L > LT , one recovers the time for the HTO wave formation by an order of magnitude as

tH ∼

1

( vH − vT )2

.

(5.49)

As vH − vT vT v, expression (5.49) typically yields very large values of tH , considerably exceeding the time in which the combustion wave traverses from the injection to the recovery well in practical applications. This argument shows that the HTO wave does not form at relevant times of our problem. Finally, let us comment on the possible effect of mass diffusion in the gaseous phase, which was neglected in the model. The size of gas diffusion layer can be estimated, in dimensional variables, as ldiff ∼ ϕ Dg /ug , which is much smaller than the size of ther-

mal diffusion layer lthermal ∼ λ/(Cm v ) because Dg ≈ λ/Cm and the wave speed v is much smaller than the gas speed ug /ϕ [18]. Moreover, the diffusion transport only affects the region upstream of the resonance point, while the combustion parameters are determined mainly by the processes in the downstream region (dominated by the equilibrium of gaseous and oleic phases), as we demonstrated in Section 5. Thus, the effect of mass diffusion may be expected to be small. Also, this agrees with the numerical simulations presented below, where diffusion was taken into account. eff

eff

6. Numerical results and simulations Here we present some results based on the theory developed in this work and compare them with numerical simulations. We

60

M.A. Endo Kokubun et al. / Combustion and Flame 169 (2016) 51–62 Table 1 Values of reservoir parameters used in numerical simulations. Arl = Arm cg eff Cm Co El kvl n Pres Qrl = Qrm

= = = = = = = = =

4060 1/s 29 J/mol K 2 MJ/m3 K 1.53 MJ/m3 K 8364 J/mol 100 1/s 1 106 Pa 400 kJ/mol O2

Qv sor sini o Tlac Tmac Tbn Tini uinj Ykin j

= = = = = = = = =

Table 2 Maximum temperatures Tu and speeds v of the combustion (MTO) wave at different pressures P. Pressure (bar)

Temperature (K)

Wave speed (m/s)

10 40 70

383 430 454

1.67 × 10−7 4.33 × 10−7 6.43 × 10−7

use the liquid, gas and rock parameters given in Table 1, taken from [22], which correspond to air injection at low rates into a rock of medium permeability filled with a two pseudo-component oil. Physical properties of the oil are chosen by considering heptane as a light component and we take hexadecane as a model medium component, with all the physical properties of hexadecane, except for its higher viscosity, thus avoiding the necessity of using complex mixtures for the model oil. Liquid combustion can be described in an analogous way to pore diffusion [31]. Then, the activation energies are chosen by taking half of their values for gaseous reaction [23,31,32]. Numerical simulations were carried out using the finite element method implemented by COMSOL software. For such, we introduce the model equations in their weak form. We use fifth-order Lagrange elements. The grid size is 0.025 m with an adaptive time step, which is appropriate for capturing the multiscale processes in the problem. In the simulations we added the terms describing effective molecular (Dg = 10−6 m2 /s) and capillary (Dcap = 5 × 10−8 m2 /s) diffusion into the governing equations in order to improve numerical convergence. The MTO temperature and wave speed are obtained as described in Section 5.3 using Eqs. (5.46), (5.47) and (5.25). Considering different reservoir pressures, i.e., P = 10, 40 and 70 bar, the dimensional results for maximum temperatures and combustion wave speeds are given in Table 2. Both the temperature and the speed increase with the pressure. Note that the combustion wave speed remains comparable to the air injection Darcy velocity uin j = 8 × 10−7 m/s. This fact is a consequence of the phasetransfer mechanism leading to strongly enhanced liquid–gas flow, a phenomenon that distinguishes the MTO process from HTO. The theoretical values in Table 2 are independent on the initial oil composition, as we explained in Section 5.3, because the medium component is expelled from the reaction region and does not appear in the resonance condition. This surprising prediction agrees reasonably well with the results of numerical simulations: in the case of a pressure of P = 10 bar, almost identical temperatures T u = 396, 398, 401, 406 K and speeds v = 1.08 × 10−7 , 1.06 × 10−7 , 1.01 × 10−7 , 9.4 × 10−8 m/s are obtained for different medium ini = 0.2, 0.4, 0.6, 0.8 in the initial oil, respectively. We fractions Xm note, however, that our model due to the assumption of absence of cracking of the heavy component is not applicable for a very heavy mixture. The discrepancies in the wave parameters obtained theoretically and numerically are due to approximations used in the derivation as well as due to numerical errors. This leads to some decrease on the wave speed and an increase on the combustion temperature, when compared to the results given by the theory.

31.8 kJ/mol 0.1 0.9 7066 K 15, 550 K 371 K 300 K 8 × 10−7 m/s 0.21

λ μlo ν gm ν om ν ol ν gl ρl ρm ϕ

= = = = = = = = =

3 W/m K 1.32 × 10−2 Pa s 1.36 0.065 0.090 1.36 6826 mol/m3 5130 mol/m3 0.3

However, the numerical simulations fully confirm the overall structure of the combustion wave. Though not affecting the temperature and speed, the initial oil composition strongly influences the oil saturation profile. In Fig. 5(a) we show the dependence of the downstream oil saturaini . This saturation remains tion sdo on the initial oil composition Xm ini , increasing when the fraction of nonlarge for all values of Xm volatile component gets larger. The results for three different pressures, shown by solid (10 bar), dashed (40 bar) and dotted (70 bar) lines, demonstrate a rather weak pressure dependence. The flow at the downstream side depends on the initial oil composition because of the change in the oil mobility, quantified by the change in the oleic viscosity in Eq. (2.13). As the initial fraction of medium oil in the oleic mixture increases, the mobility of the oil decreases. Thus, in order to keep the wave speed unchanged, the downstream ini leading to a considerable insaturation must increase with Xm crease of gas drive for small gas saturations sdg = 1 − sdo . Comparing the theoretical results in Fig. 5(a) with the values obtained by numerical simulations, shown by crosses, we see a reasonable agreement with the theory. Upstream of the combustion wave, only the medium (immobile) oil is left behind. Figure 5(b) shows the upstream oil saturation suo given by Eq. (5.24). As predicted, this value is small for all pressures and initial oil compositions. Still, this small residual oil may react with some fraction of oxygen in the injected air, since it is located at the hot upstream part of the combustion wave. Figure 5(c) shows the estimate (5.48) of this oxygen fraction Ykmax . We see that the oxygen consumption by the residual medium oil is negligible for larger pressures, even if the initial oleic mixture is predominantly non-volatile. Only for lower pressures this value may become considerable in comparison with the in j injected fraction Yk = 0.21.

The large downstream oil saturations sdo in Fig. 5(a) lead to high flow rates of oil downstream of the combustion wave, while the small upstream oil saturations suo in Fig. 5(b) imply almost complete recovery of oil from the reservoir. Also the relatively low temperature increase in the wave requires less oil to be burned compared to the HTO process for heavy oils. All these factors make the described combustion (MTO) mechanism attractive for applications. Now we describe the results of numerical simulations focusing on the case of a pressure P = 10 bar and initial oil with medium ini = 0.4. Figure 6 shows the dependent variables at a fraction Xm specific time t = 1.2 × 109 s when the stationary combustion process is developed. Dimensional temperature T, oxygen fraction Yk and oil saturation so are shown in Fig. 6(a), while the oleic medium fraction Xm and gaseous light fraction Yl are shown in Fig. 6(b) together with the products Xl so and Xm so describing the light and medium oil saturations. One can recognize the combustion wave in the middle (around x ≈ 110 m) and the slower thermal wave forming in the upstream (left) part of the domain; the faster saturation wave moved out of the computation region at earlier times. The combustion wave pushes the oleic mixture forward and is responsible for displacing the non-volatile oil.

M.A. Endo Kokubun et al. / Combustion and Flame 169 (2016) 51–62

61

Fig. 6. Combustion wave profile for Xmini = 0.4 and P = 10 bar. (a) Temperature T, oxygen fraction Yk and oil saturation so . (b) Medium oil saturation Xm so , light oil saturation Xl so , gaseous hydrocarbon Yl and oleic medium fraction Xm . The arrows indicate the resonance point (RP).

One can see a good agreement of numerical simulation results with the combustion wave structure described theoretically in Section 5, see also Fig. 2. The vaporization/condensation front and the reaction front travel together, inside the combustion wave, with constant speed v. The arrows in Fig. 6(a,b) indicate the location of the resonance point (RP); the inset in Fig. 6(b) shows a zoom of the figure near this point. In agreement with theoretical results, one can see that the oleic medium fraction Xm attains a minimum at the resonance point, while the fraction Yl of the gaseous light component generated by the vaporization has a sharp peak at the resonance point, decreasing on the downstream side because of condensation at lower temperatures. A small amount of medium component Xm so remains behind the resonance point, see the inset in Fig. 6(b). This residual medium oil does not burn because of the high value of its activation energy. We note that even if the described traveling wave solution is qualitatively the same for all initial oil compositions, the transient process of the combustion wave formation is very different for oils of high or low volatility. A bank of immobile medium oil appears ini is large, for at the initial stage of the combustion process if Xm ini = 0.8. Combustion of this oil leads to high temperexample, Xm atures reaching typical values for the HTO process. This fact, observed earlier in numerical simulations [22], was interpreted as a bifurcation from the MTO to the HTO combustion with an increase of the non-volatile component in the initial oil. Our simulations, performed for larger time intervals, showed that the MTO combustion wave is formed eventually after (probably long) transient behavior; the bifurcation is not observed in the current model, as stated previously. In view of the results obtained, more detailed models are necessary to see the transition to HTO, which must include the cracking process (pyrolysis) leading to coke deposition at high temperatures.

saturations are attained downstream of the combustion wave, leading to considerable increase in the wave speed. These aspects are demonstrated analytically by studying a two-component liquid model originating from the problem of light oil recovery by air injection, where one pseudo-component is volatile and has a high mobility (called light), while the other pseudo-component is non-volatile and immobile (called medium). The dominant role of the light component in the combustion process has several practical consequences: the total temperature increase remains small being limited by the light oil boiling point, no oil is left behind the wave, while the downstream oil flux may reach magnitudes comparable to the injected air flux. These properties make the described combustion mechanism attractive for oil recovery, in comparison with the high-temperature oxidation process applied for heavy oils, because of lower oil consumption in the reaction and applicability to abandoned reservoirs with the potential of recovering residual oil. Note that the suggested mechanism of oil recovery is not just due to decrease of oil viscosity, but it is related to a complex interplay between the flow and phase transition mechanism. Our results allow to make quantitative estimates for recovery rates based on the conditions downstream of the combustion wave. For example, the downstream oil saturation ranges between 0.5 and 0.9. Theoretically, we found a qualitatively new combustion mechanism controlled by the successive vaporization and condensation of the liquid phase sustained by the reaction. Since the multicomponent effect is dominant in this process, our analysis is crucial for understanding filtration combustion of light oils, which are intrinsically multicomponent. However, our results show that further study is necessary in order to fully understand the complex multicomponent phenomena in this process, by considering models with a large number of components having variable physical properties.

7. Conclusions Acknowledgments For filtration combustion of a volatile multicomponent liquid fuel, multiphase flow and phase transitions become as important for the development of the combustion wave, as the oxidation mechanism. In addition to this fact, we demonstrated that the multicomponent aspects of the flow are equally important, if the liquid components have essentially different physical properties, leading to qualitative changes in the combustion process. In particular, a component with low mobility and volatility is expelled from the reaction region in the downstream direction, opposite to what would be intuitively expected. As a result, large liquid

This work was supported by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) under Grant no. 302351/2015-9 and by Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ) under Grant nos E26/201.210/2014, E-26/110.658/2012, E-26/110.114/2013 and Pensa Rio E-26/210.874.2014. The first author acknowledge the financial support by Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) and Instituto de Matemática Pura e Aplicada (IMPA).

62

M.A. Endo Kokubun et al. / Combustion and Flame 169 (2016) 51–62

References [1] A. Aldushin, Filtration combustion, in: P. Zarchan (Ed.), Advances in Combustion Science: In Honor of Ya. B. Zel’dovich, AIAA, Reston, USA (1997), pp. 95–115. [2] A. Mailybaev, J. Bruining, D. Marchesin, Analytical formulas for in-situ combustion, SPE J. 16 (3) (2011) 513–523. [3] D. Schult, B. Matkowsky, V. Volpert, A. Fernandez-Pello, Forced forward smolder combustion, Combust. Flame 104 (1996) 1–26. [4] A. Aldushin, B. Matkowsky, D. Schult, Downward buoyant filtration combustion, Combust. Flame 107 (1996) 151–175. [5] I. Akkutlu, Y. Yortsos, The dynamics of in-situ combustion fronts in porous media, Combust. Flame 134 (2003) 229–247. [6] C. Wahle, B. Matkowsky, A. Aldushin, Effects of gas–solid nonequilibrium in filtration combustion, Combust. Sci. Technol. 175 (2003) 1389–1499. [7] A. Aldushin, A. Bayliss, B. Matkowsky, On the transition from smoldering to flaming, Combust. Flame 145 (2006) 579–606. [8] A. Mailybaev, D. Marchesin, J. Bruining, Resonance in low-temperature oxidation waves for porous media, SIAM J. Math. Anal. 43 (2011) 2230–2252. [9] A. Mailybaev, J. Bruining, D. Marchesin, Recovery of light oil by medium temperature oxidation, Transp. Porous Media 97 (2013) 1–27. [10] N. Gargar, A. Mailybaev, D. Marchesin, H. Bruining, Effects of water on light oil recovery by air injection, Fuel 137 (2014) 200–210. [11] N. Gargar, A. Mailybaev, D. Marchesin, H. Bruining, Recovery of light oil by air injection at medium temperature: experiments, J. Pet. Sci. Technol. 133 (2015) 29–39. [12] P. Sarathi, In-situ combustion handbook – principles and practices, BDM Petroleum Technologies, Bartlesville, OK, USA, 1999. [13] L. Castanier, W. Brigham, Upgrading of crude oil via in situ combustion, J. Pet. Sci. Technol. 39 (2003) 125–136. [14] S.M. Farouq Ali, Heavy oil-evermore mobile, J. Pet. Sci. Eng. 37 (2003) 5–9. [15] M. Onyekonwu, G. Falade, Recovery of light oil using in situ combustion thermal recovery methods, Energy 14 (1989) 153–159. [16] S. Ren, M. Greaves, R. Rathbone, Air injection LTO process: an IOR technique for light-oil reservoirs, SPE J. 7 (1) (2002) 90–99.

[17] L. Adetunji, R. Teigland, Light-oil air-injection performance: sensitivity to critical parameters, SPE Annual Technical Conference and Exhibition (2005). Paper SPE-96844. [18] A. Mailybaev, J. Bruining, D. Marchesin, Analysis of in situ combustion of oil with pyrolysis and vaporization, Combust. Flame 158 (6) (2011) 1097–1108. [19] M. Pascual, D. Crosta, P. Lacentre, D. Coombe, Air injection into a mature waterflooded light oil reservoir–laboratory and simulation results for Barrancas field, Argentina (SPE94092), 67th EAGE Conference & Exhibition (2005). [20] C. Clara, M. Durandeau, G. Quenault, T. Nguyen, Laboratory studies for light-oil air injection projects: potential application in handil field, SPE Reserv. Eval. Eng. 3 (3) (20 0 0) 239–248. [21] M. Greaves, S. Ren, R. Rathbone, T. Fishlock, R. Ireland, Improved residual light oil recovery by air injection (LTO process), J. Can. Pet. Technol. 39 (1) (20 0 0). [22] N. Gargar, A. Mailybaev, D. Marchesin, H. Bruining, Compositional effects in light/medium oil recovery by air injection: vaporization vs. combustion, J. Porous Media 17 (2014) 937–952. [23] O. Levenspiel, Chemical reaction engineering, Wiley, 1999. [24] E. Koval, A method for predicting the performance of unstable miscible displacement in heterogeneous media, SPE J. 3 (2) (1963) 145–154. [25] B. Poling, J. Prausnitz, O. John Paul, R. Reid, The properties of gases and liquids, McGraw-Hill, New York, 2001. [26] N. Freitag, B. Verkoczy, Low-temperature oxidation of oils in terms of SARA fractions: why simple reaction models don’t work, J. Can. Pet. Technol. 44 (3) (2005) 54–61. [27] F. Moritis, R. Bell, The low temperature liquid phase oxidation of hydrocarbons: a literature survey, J. Inst. Pet. 44 (1958) 260–272. [28] I. Prigogine, Introduction to thermodynamics of irreversible processes, Wiley, 1967. [29] F. Santos, A. Mailybaev, D. Marchesin, Oxidation wave structure and oxygen breakthrough for air injection into light oil reservoirs, Comput. Geosci. (2016). In Press [30] A. Aldushin, I. Rumanov, B. Matkowsky, Maximal energy accumulation in a superadiabatic filtration combustion wave, Combust. Flame 118 (1999) 76–90. [31] R. Bird, W. Stewart, E. Lightfoot, Transport phenomena, 2002. [32] P. Harriott, Chemical reactor design, 2003.