TRAVELING WAVES FOR IN-SITU COMBUSTION IN POROUS MEDIA ONDAS VIAJANTES PARA COMBUSTÃO IN-SITU EM MEIOS POROSOS

This work presents a mathematical model describing the in-situ combustion process, which can be used in enhanced oil recovery. The hyperbolic part of the system was solved as a Riemann Problem. Necessary conditions for the existence of a traveling wave solution were verified. Furthermore, theoretical results are verified by using numerical simulations.


INTRODUCTION 2 FRACTIONAL FLOW THEORY
One way to model the fluid injection and oil production is through the Buckley-Leverett theory (BUCKLEY;LEVERETT, 1942), which allows the analysis of unidimensional and twophase flow in a porous medium. According to this model, after the fluid injection, the oil from the reservoir is displaced; however, its saturation is never zero, since a portion of residual oil is always kept in pores. Likewise, after the injection, the fluid cannot be completely removed, so that its saturation will never be zero. Thus, the interval in which the relative saturation of the injected fluid varies is [S f c , 1 − S or ], where S f c is the residual fluid saturation and S or is the residual oil saturation. In this case, S o and S g are oil and gas saturation, respectively. It will be considered that S o and S g are normalized saturations, that is, they belong to the interval [0, 1] and S g = 1 − S o . The fractional flow of a phase is defined as the flow of this phase in relation to the total flow of fluids being produced. Thus, the functions that describe fractional flows of oil and gas can be defined as: where λ o and λ g are oil and gas mobilities. Note that mobility depends on the relative permeability, commonly described by the Brooks-Corey model (BROOKS;COREY, 1964).
In this case, f o and f g are functions of saturation. For purposes of simplifying the model, constant relative permeabilities and therefore constant functions will be considered, maintaining the relationship f o + f g = 1. Thus, the Darcy velocities of oil and gas phases are given by are also constants and thus u o + u g = u, where u is the Darcy velocity of the resulting fluid, composed of oil and gas.

MATHEMATICAL MODELING OF THE IN-SITU COMBUSTION
The model proposed in this paper is similar to the model proposed in (AKKUTLU;YORTOS, 2003). In this case, it will be considered the burning of oil instead of solid fuel, disregarding the effects of diffusion. It is also assumed that the flow is one-dimensional, with air (oxygen) injected to the left of a porous rock cylinder containing a non-reactive gas and a mobile fuel.
Some additional simplifications have been adopted. The molar density of gas ρ g is assumed constant and approximately equal to the oxygen molar density ρ o . In the heat balance Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 6, p. 278-01, 278-14, 2020 DOI: 10.21575/25254782rmetg2020vol5n61289 equation, the term that accompanies the partial derivative with respect to x is (c o ρ o S o u o + c g ρ g S g u g )(T − T res ), where T and T res represent local and initial temperatures of the reservoir, c o and c g -thermal capacities of oil and gas, C m -thermal capacity of the porous medium.
However, the thermal capacity of oil is much higher than the thermal capacity of gas, as well as the density. Then, the term c g ρ g S g u g is considered negligible compared to c o ρ o S o u o and, therefore, the equation is as in (3). It was also assumed that porosity (ϕ) and Darcy's speed (u) are constant. The model with time coordinate t and space coordinate x is composed of the equations of heat balance (3), molar balance for oxygen (4) and molar balance for mobile fuel (5), is given by: In the reaction, ω f moles of fuel react with ω o moles of oxygen generating ω g moles of gaseous products and solid waste products, which are neglected in this model. For simplicity, consider the case ω o = ω f = ω g = 1, such as in the reaction C + O 2 CO 2 . The reaction rate is proportional to W r , given by the Arrhenius' Law: where k p is a pre-exponential factor, E r is the activation energy, R is the gas constant.

Constant Parameters
Parameter values describing oil properties vary significantly from one reservoir to another due to the change in the oil composition. Here we adopt values related to mineral oil because it possesses more consistent parameter values. Thermal capacity (c o ) and average molar density (ρ o ) were taken from (OLIVEIRA; GOUVEIA, 2006;HARTWIG, 2002). As for the reservoir, it will be considered a porous rock block with length L = 20 cm (typical value for laboratory experiments) with a mean pressure gradient ∇P = 10 6 P a and permeability k = 10mD.
The Darcy velocity is estimated using: where k is the permeability of the porous rock, µ is the viscosity of the fluid, ∇P is the pressure gradient and L is the length of the porous medium. Considering that the average viscosity of the mixture is approximately equal to that of the gas (µ ≈ µ g ), for (8), we obtain u = 2.03 × 10 −3 m/s. Using (2) Table 1.

Non-dimensional Model
To perform the model's analysis, it is useful to rewrite it in non-dimensional form, preventing operations without physical meaning and assisting numerical simulations. Next some dependent and independent dimensionless variables (denoted with bars) are introduced Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 6, p. 278-01, 278-14, 2020 DOI: 10.21575/25254782rmetg2020vol5n61289 as fractions of the dimensional and reference quantities (later denoted by stars): Eq.

STUDY OF THE HYPERBOLIC PART
In this section, the hyperbolic part of the equations describing the model is studied, i.e., with the source term (right side of the equations) equal to zero. This analysis is useful because outside the combustion zone, the value of the source terms are insignificant. As in the section " Fractional Flow Theory", f o and f g are considered constants. For simplicity we remove the bars and rewrite the dimensionless equations: where:   Notice that this equation is in the form: The eigenvalues of the matrix f (U ) are λ 1 = a 1 S o , λ 2 = a 2 and λ 3 = a 3 . The eigenvectors associated with eigenvalues λ 1 , λ 2 and λ 3 are: With this information, we can investigate each of the p-fields: We conclude that all fields are Linearly Degenerate, indicating that the constant states will be connected by contact discontinuities, whose velocities are the eigenvalues λ i .

Hugoniot Locus
A discontinuity (shock) (U l , U r ) can be a solution of the Riemann Problem if it satisfies the Rankine-Hugoniot Condition f (U l ) − f (U r ) = s(U l − U r ) for some velocity s, see (LEVEQUE, 1992) for more details.
In order to obtain such solution for the system, Hugoniot Locus will be studied for the initial states U l and U r . By definition, the Hugoniot Locus of a pointÛ is given by: where s is the jump velocity established by the Rankine-Hugoniot Condition.
From the Eq. (21), we obtain: The other cases lead to same results as those found and also to the identity equations. Thus we do not consider them separately. The three shock curves are: Note that the propagation speeds of the contact discontinuities coincide with the eigenvalues. Therefore, the Hugoniot Locus ofÛ denoted by H(Û ), is given by: Substituting the values of the constants given in (17) in the speeds of the system (28)-(30), we have that the velocities of gas, oil and thermal wave are, respectively: s g = 814.45, s o = 14.21, s t = 2.34Ŝ o .
Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 6, p. 278-01, 278-14, 2020 DOI: 10.21575/25254782rmetg2020vol5n61289 8 In this section, we look for the model's solution in the form of a traveling wave. It is equivalent to search a solution that travels at a constant speed c (here we consider c > 0, i.e., the same direction as the fluid) with the same profile. It is also equivalent to search for stationary solution in traveling coordinates (ξ, t), where ξ = x − ct. Mathematically, the traveling wave solution is an orbit connecting left equilibrium to right equilibrium as α-limit and ω-limit, see (CHAPIRO; MARCHESIN; SCHECTER, 2014; VOLPERT; VOLPERT; VOLPERT, 1994) for more details. One of the necessary conditions for the existence of such orbit is the suitable dimension of these equilibria. Obtaining these dimensions is the main goal of this section.
Considering the constant flux functions in the travelling wave variables, the system formed by equations (10)-(12) is rewritten as: To continue the calculations, we assume hypotheses c = a 1 S o , c = a 2 and c = a 3 , resulting in the following system of Ordinary Differential Equations (ODEs) (for simplicity, we denote derivative with respect to ξ by primes): Isolating S o in Eq. (35) and replacing the result into (34), we obtain: As both left (θ L , S L y , S L o ) = (0, 1, 0) and right (θ R , S R y , S R o ) = (0, 0, 1) equilibria are stationary solutions, applying the limit (ξ → ±∞) and using the boundary conditions we can find the traveling wave speed: Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 6, p. 278-01, 278-14, 2020 DOI: 10.21575/25254782rmetg2020vol5n61289 Replacing the constant values given in (17) into Eq. (37), the speed of traveling wave is c = 39.72. On the other hand, from Eq. (36), we obtain: As S y is the oxygen saturation, the previous relation results in Y = 1. This physically inaccurate approximation is due to the simplicity of the system. The study of a more realistic model will remain for future work. When replacing S y , found in (33) to (35), yields: Therefore, the Jacobian matrix J of the system (39) is given by: where At the equilibrium (θ L , S L o ) = (0, 0) and (θ R , S R o ) = (0, 1), the Jacobian matrix (40) becomes: where the eigenvalues are λ L 1 = 0 and λ L 2 = b 3 Φ/M > 0 and λ R 1 = 0 and λ R 2 = −(b 3 Φ)/M < 0.
Note that both equilibria possess eigenvalues equal to zero, classified as nodes. Rigorous demonstration of the existence of an orbit connecting these equilibria stays outside of the scope of this work.

CONCLUSIONS
This paper proposes a simplified model of the in-situ combustion and discusses several simplifications that can be adopted. The solution of the Riemann problem was obtained as a sequence of three contact waves and one traveling wave. Numerical simulations have shown profiles compatible with the theoretical analysis.