Integrability of Buckley – Leverett’s Filtration Model*

IFAC Conference on Manufacturing Modelling, IFAC Conference Manufacturing Modelling, Management and on Control IFAC Conference Conference on Manufacturing Modelling, Modelling, IFAC Manufacturing Management and on Control Available online at www.sciencedirect.com June 28-30, 2016. Troyes, France Management and Control Management and Control June 28-30, 2016. Troyes, France June 28-30, 28-30, 2016. 2016. Troyes, Troyes, France France June

ScienceDirect

IFAC-PapersOnLine 49-12 (2016) 1251–1254

Integrability of Buckley – Integrability of Buckley – Integrability of Buckley – Filtration Model Filtration Model Filtration Model

Leverett’s Leverett’s Leverett’s ⋆⋆ ⋆

Atlas V. Akhmetzianov ∗∗ Alexei G. Kushner ∗∗ ∗∗ Atlas ∗ Alexei G. ∗∗ ∗∗∗ Kushner ∗ Atlas V. V. Akhmetzianov Akhmetzianov Alexei G. Kushner ∗∗ Valentin V. Lychagin Atlas V. Akhmetzianov Alexei G. ∗∗∗ Kushner Valentin V. Lychagin ∗∗∗ ∗∗∗ Valentin V. Lychagin Valentin V. Lychagin ∗ ∗ Institute of Control Sciences of the Russian Academy of Sciences, Institute of Control Sciences of the Russian Academy of Sciences, ∗ ∗ Institute of Sciences of Russian Academy Russian Federation, Moscow (e-mail: [email protected]). Institute of Control Control Sciences of the the Russian Academy of of Sciences, Sciences, Russian Federation, Moscow (e-mail: [email protected]). ∗∗ Russian Federation, Moscow (e-mail: [email protected]). Institute of Control Sciences of the Russian Academy of Sciences, Russian Federation, Moscow (e-mail: [email protected]). ∗∗ Institute of Control Sciences of the Russian Academy of Sciences, ∗∗ ∗∗ Institute of Control Control Sciences Sciences ofMoscow the Russian Russian Academy of(e-mail: Sciences, Russian Federation, Moscow; of State Academy Universityof Institute of the Sciences, Russian Federation, Moscow; Moscow State University (e-mail: Russian Federation, Moscow; Moscow State University (e-mail: [email protected]) Russian Federation, Moscow; Moscow State University (e-mail: [email protected]) ∗∗∗ [email protected]) Sciences of the Russian Academy of Sciences, [email protected]) ∗∗∗ Institute of Control Institute of Control Sciences of Academy of ∗∗∗ ∗∗∗ Institute of Control ControlMoscow; SciencesTromso of the the Russian Russian Academy of Sciences, Sciences, Russian Federation, University, Tromso, Norway Institute of Sciences of the Russian Academy of Sciences, Russian Federation, Moscow; Tromso University, Tromso, Russian Federation, Federation, Moscow; Tromso University, University, Tromso, Tromso, Norway Norway (e-mail: [email protected]) Russian Moscow; Tromso Norway (e-mail: [email protected]) (e-mail: (e-mail: [email protected]) [email protected]) Abstract: One-dimensional Buckley – Leverett’s model of frontal displacement of petroleum Abstract: One-dimensional Buckley –– Leverett’s model oftofrontal displacement of petroleum Abstract: One-dimensional Buckley model displacement of by water is considered. We apply Backlund transformation the Buckley – Leverett equation. Abstract: One-dimensional Buckley – Leverett’s Leverett’s model of oftofrontal frontal displacement of petroleum petroleum by water is considered. We apply Backlund transformation the Buckley – Leverett equation. by water is considered. We apply Backlund transformation to the Buckley – Leverett equation. This allows us to find Cauchy problem’s solutions for this equation in quadratures. These results by water is considered. We apply Backlund transformation to the Buckley – Leverett equation. This allows us to find Cauchy problem’s solutions for this equation in quadratures. These Thisbe allows us to find Cauchy Cauchy problem’sofsolutions solutions for this this equation equation in in quadratures. quadratures. These These results results can usedus into control of development oil fields.for This allows find problem’s results can be used in control of development of oil fields. can be used in control of development of oil fields. can be used in control of development of oil fields. © 2016, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Frontal displacement, Backlund transformation, Monge-Ampere equation, Keywords: Frontal Frontal displacement, Backlund transformation, Monge-Ampere equation, Keywords: displacement, multivalued solutions, differentialBacklund forms, jettransformation, space, LegendreMonge-Ampere transformation.equation, Keywords: Frontal displacement, multivalued solutions, differentialBacklund forms, jet jettransformation, space, Legendre LegendreMonge-Ampere transformation.equation, multivalued solutions, differential forms, space, transformation. multivalued solutions, differential forms, jet space, Legendre transformation. 1. INTRODUCTION 1. INTRODUCTION INTRODUCTION 1. 1. INTRODUCTION Modelling of oil recovery process by water-flooding can be Modellingby of the oil recovery recovery process by water-flooding water-flooding can be be Modelling of oil process by can described differential equation of Buckley – Leverett Modelling of the oil recovery process by water-flooding can be described by differential equation of Buckley – Leverett described described by by the the differential differential equation of Buckley Buckley –– Leverett Leverett mst + U F ′′equation (s)sx = 0.of (1) ms + U F (1) ′′ (s)sx = 0. t + U F (s)s = 0. (1) ms t x mst + U F (s)sx = 0. (1) Here s = s(t, x) is the relative volume of water in pore Here s = s(t, x) is the relative volume of water in pore Here = x) of in space ss(water m is a volume coefficient of a capillary Here = s(t, s(t,saturation), x) is is the the relative relative volume of water water in pore pore space (water saturation), m is a coefficient of a capillary space (water saturation), m is a coefficient of a capillary porosity, U is the sum of filtration velocities of water and space (water saturation), m is a coefficient of a capillary porosity, U is the sum of filtration velocities of water and porosity, U is the sum of filtration velocities of water and oil, F (s) is the Buckley – Leverett function. The graph of porosity, U the is the sum of– filtration velocities The of water and oil, F (s) is Buckley Leverett function. graph oil, F Ffunction (s) is is the the Buckley – Leverett Leverett function. function. The The graph graph of of this seeBuckley on Fig.1. oil, (s) – of this function see on Fig.1. this function see on Fig.1. this function see on Fig.1. By Darcy’s law By By Darcy’s Darcy’s law law By Darcy’s law U = −κh(s)q, (2) U (2) U= = −κh(s)q, −κh(s)q, (2) U = −κh(s)q, (2) where κ is the oil formation permeability coefficient, q is where κ is the formation coefficient, q is where κ oil formation permeability coefficient, the gradient of oil pressure: q = ppermeability x (t, x), where κ is is the the oil formation permeability coefficient, qq is is the gradient of pressure: q = p (t, x), x the gradient gradient of of pressure: pressure: qqf = =p pxx (t, (t, fx), x),(s) the w (s) o + ffo (s) , h(s) = ffw (s) w o (s) h(s) = fµ fw (s) (s) µ(s) o w + w ,, Fig. 1. The graph of the Buckley – Leverett function. + h(s) = , h(s) = µw + µ 1. The graph of the Buckley – Leverett function. µw w Fig. µw and oilµ penetrability respec- Fig. fw (s) and fo (s) are waterµ w w Fig. 1. 1. The The graph graph of of the the Buckley Buckley –– Leverett Leverett function. function. (s) and f (s) are water and oil penetrability respecftively, w o does not depend on coordinate x, i.e. the total derivative (s) and f (s) are water and oil penetrability respecffw µ and µ dynamic viscosities of water and oil o w o (s) and f (s) are water and oil penetrability respecdoes not depend on coordinate x, i.e. the total derivative w o tively, µ viscosities of water and oil on does depend on w and o dynamic x not of the function U is zero: x, not depend on coordinate coordinate x, i.e. i.e. the the total total derivative derivative tively, µ and µ dynamic Entov, viscosities of and respectively (seeµ Ryzhik (2010)). o tively, µw µBarenblatt, viscosities of water water and oil oil does on x of the function U is zero: w and o dynamic Entov, respectively (see Barenblatt, Ryzhik (2010)). on dUzero: respectively (see (see Barenblatt, Barenblatt, Entov, Entov, Ryzhik Ryzhik (2010)). (2010)). on x x of of the the function function U U is is zero: respectively dU The Buckley – Leverett function has the form = 0. (3) dU The has the form 0. (3) dU dx = The Buckley Buckley – – Leverett Leverett function function has the form = 0. (3) dx The Buckley – Leverett function has the form (3) fw (s) dx = 0. ffw (s) dx , F (s) = (s) w Thus, we consider the following system of two partial F (s) = fw (s)fw +(s) µfo (s) ,, F Thus, following system of + µf , F (s) (s) = = ffw (s) o (s) Thus, we we consider consider the following system(1): of two two partial partial differential equationsthe instead of equation (s) + µf (s) o Thus, we consider the following system of two partial fw (s) + µf (s) where w o differential equations instead of equation (1):  µ where differential equations instead of equation (1): w differential equations instead of′ equation (1): where  . µ := µ where  t − kh(s)qF ′ (s)sx = 0,  ms µw w  µ := µ  ms kh(s)qF ′′ (s)s w o .. t− x = 0, µ :=  ms − (4) tt − kh(s)qF µ := µ  ms o . kh(s)qF (s)s (s)sxx = = 0, 0, µ (4) o  ′ µ o (4) h (s)qs + h(s)q = 0, In Buckley – Leverett’s model water and oil are considered x x  ′ (4)  h (s)qs + h(s)q = 0, ′ In Buckley – Leverett’s model water and oil are considered x x  ′ h (s)qs + h(s)q = 0, In Buckley – Leverett’s model water and oil are considered x x incompressible. Therefore the water sum of filtration velocities h (s)qsx + h(s)qx = 0, In Buckley – Leverett’s model and oil are considered incompressible. Therefore the sum of filtration velocities incompressible. Therefore the sum filtration velocities initial and boundary conditions: incompressible. Therefore theRussian sum of ofScience filtration velocities with ⋆ This work was supported by the Foundation with initial and boundary ⋆ This work was supported by the Russian Science Foundation (the with and boundary conditions: conditions: (the with initial initial and ⋆ q(0, x) = q0 (x). (5) grants 15-19-00275). s(0, x) boundary = s0 (x), conditions: This No work was supported supported by by the the Russian Russian Science Science Foundation Foundation (the (the ⋆ This work was (x), q(0, x) = qq0 (x). (5) grants No 15-19-00275). s(0, x) = s 0 (x), q(0, x) = (x). (5) grants No 15-19-00275). s(0, x) = s 0 0 q(0, x) = q0 (x). (5) grants No 15-19-00275). s(0, x) = s0 (x), Copyright 2016 IFAC 1251Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © 2016, IFAC (International Federation of Automatic Control) Copyright 2016 responsibility IFAC 1251Control. Peer review© of International Federation of Automatic Copyright © 2016 1251 Copyright ©under 2016 IFAC IFAC 1251 10.1016/j.ifacol.2016.07.685

IFAC MIM 2016 1252 June 28-30, 2016. Troyes, France

Atlas V. Akhmetzianov et al. / IFAC-PapersOnLine 49-12 (2016) 1251–1254

The corresponding effective differential 2-form on the space of 1-jets J 1 (R2 ) (see Kushner, Lychagin, Rubtsov (2007)) is ω = 2u1 du2 ∧ dx1 + dx1 ∧ du1 − dx2 ∧ du2 . Here x1 = t, x2 = x, u, u1 , u2 are canonical coordinates on the jet-space J 1 (R2 ) (see Krasil’shchik, Lychagin, Vinogradov (1986)). Two Laplace forms of this equation are zero (see Kushner (2009)). Therefore equation (9) can be integrated in quadratures. To construct its solution, we apply Legendre transformation (see Kushner, Lychagin, Rubtsov (2007)) ϕ : (x1 , x2 , u, u1 , u2 ) �−→ (−u1 , −u2 , u−x1 u1 −x2 u2 , x1 , x2 ), which takes the form ω to the form ϕ∗ (ω) = −2x1 du1 ∧ dx2 + dx1 ∧ du1 − dx2 ∧ du2 . The last form corresponds to linear hyperbolic equation ux1 x2 − x1 ux1 x1 = 0. Fig. 2. The graph of the function Φ.

General solution of this equation is

and

u(x1 , x2 ) = e−x2 F1 (x1 ex2 ) + F2 (x2 ), 0

s(t, 0) = s (t),

0

q(t, 0) = q (t).

(6)

We suppose that the pumped well is located at the point x = 0. Boundary conditions (6) define the mode water injection into the oil reservoir. Note 1. Typically, the pressure gradient and the saturation of water are kept constant at pumped well. This means that the functions s0 and q 0 in (6) could be constant. Here we describe a method of constructing explicit solutions of system (4). This article can be considered as a continuation of our work on the application of differential geometric methods in control theory of oil fields development (see Akhmetzianov, Kushner, Lychagin (2013-a), Akhmetzianov, Kushner, Lychagin (2013-b), Akhmetzianov, Kushner, Lychagin (2015)). 2. BACKLUND TRANSFORMATION AND INTEGRABILITY Due to equations (3) and (1), we get: sx stx − st sxx F ′′ (s) . = ′ 2 st sx F (s) The last equation can be written in the following form: � � ∂ ∂ st (7) = Φ(s), ∂x sx ∂t where The graph of the function Φ see on Fig.2.

(8)

to this equation we get the following hyperbolic MongeAmpere equation: utx − ut uxx = 0.

(11)

where F1 and F2 are arbitrary functions. Applying inverse transformation ϕ−1 to solution (11) we get multivalued solutions (see Kushner, Lychagin, Rubtsov (2007)) of equation (9) in parametric form:  t = F1′ (−u1 e−u2 ),         x = −eu2 F1 (−u1 e−u2 ) − u1 F1′ (−u1 e−u2 )+ F2′ (−u2 ), (12)      ′  u2 )F1′ (−u1 e−u2 )+   u = uu22F2 (−u2 ) + u1 (1 −−u e (1 − u2 )F1 (−u1 e 2 ) + F2 (−u2 ). Let’s introduce two parameters a := −u1 e−u2 , b := −u2 , and two functions k(a) := F1 (a), r(b) := F2 (b). Then the last formula can be rewritten in compact form:  t = k ′ (a),         x = e−b (ak ′ (a) − k(a)) + r′ (b),     u = −br′ (b) − (1 + b)e−b (ak ′ (a) − k(a)) + r(b), (13)       u1 = −ae−b ,       u2 = −b.

From the geometrical point of view formula (13) gives a surface L in 5-dimensional space J 1 (R2 ). It is not hard to check that a restriction of the differential form ω on L is zero. Therefore the surface L is a multivalued solution of equation (9).

Φ(s) := ln F ′ (s).

Applying Backlund transformation st ut = , ux = Φ(s) sx

(10)

(9)

The projection

π : L → R2 (t, x) of this surface to the plane of independent variables can have singular points that form a curve Σ on the surface (caustic).

1252

IFAC MIM 2016 June 28-30, 2016. Troyes, France

Atlas V. Akhmetzianov et al. / IFAC-PapersOnLine 49-12 (2016) 1251–1254

1253

The set π0 (L \ Σ) consists of graphs of classical solutions. Here π0 : L → R3 (t, x, u). is a projection. Shock waves (discontinuous solutions) are formed near caustics. Any solution of equation (1) can be found by second part of formula (8). To do this one should solve functional equation ux = Φ(s) with respect to the function s. 3. INITIAL AND BOUNDARY CONDITIONS Initial and boundary conditions (5) and (6) generate conditions for the function u. Let us find them. It follows from (5) that ux |t=0 = φ(x),

(14)

Solving this equation we find b as a function of z: b = B(z) and boundary condition (15) takes the form

where φ(x) := Φ(s0 (x)). Howewer, from equation (6) and (4) we get:   s| = (s0 (t))′ ,   t x=0 0

   sx |x=0 =

′

−ae−B(ak (a)−k(a)+ξ) = γ(k ′ (a)).

′

m(s (t))

κh(s0 (t))F ′ (s0 (t))q 0 (t)

Fig. 3. The curve b + φ(−ze−b ) = 0.

This is a first order ordinary differential equation with respect to function k.

.

Therefore the boundary condition for the function u is (15) ut |x=0 = γ(t), where

κ h(s0 (t))F ′ (s0 (t))q 0 (t), m and one should find the functions k = k(a) and r = r(b) such that above conditions (14) and (15) are satisfied. γ(t) :=

Suppose that the initial moment of time corresponds to value a0 of the parameter a, i.e. k ′ (a0 ) = 0. Then condition (14) generates the equation b + φ(r′ (b) − ξe−b )) = 0, where ξ := k(a0 ). This equation can be solved with respect to r′ on intervals of monotony of the function φ: r′ (b) = ξe−b + φ−1 (−b). (16) The image of line x = 0 on the (a, b)-plane is the curve e−b (ak ′ (a) − k(a)) + r′ (b) = 0,

(17)

or e−b (ak ′ (a) − k(a) + ξ) + φ−1 (−b) = 0. Finally, we get an equation on the parameter b: b + φ(−ze−b ) = 0,

(19)

Let us consider the case when γ is a constant function (see Note 1). Then on intervals of monotony of the function B equation (19) can be solved with respect to the derivative k′ : ak ′ (a) − k(a) + ξ = g(a), (20) where g(a) := B −1 (ln(−a/γ)) and k(a) = a

�a

a0

ξ g(τ ) − ξ dτ + a. τ2 a0

Let us find the value of a0 . If we put a = a0 in formula (20) we get g(a0 ) = 0, i.e. B −1 (ln(−a0 /γ)) = 0. Then B(0) = ln(−a0 /γ) and we get a0 = −γeB(0) . Note that the value of ξ can be arbitrary. We take ξ = 0. Then �a g(τ ) dτ. (21) k(a) = a τ2 a0

(18)

z := ak (a) − k(a) + ξ.

4. THE METHOD TO CONSTRUCT SOLUTIONS OF BUCKLEY – LEVERETT’S EQUATION

The curve b + φ(−ze−b ) = 0 on the (z, b)-plane see on Fig.3.

Here we describe the method to construct solutions of Buckley – Leverett’s equation.

where ′

1253

IFAC MIM 2016 1254 June 28-30, 2016. Troyes, France

Atlas V. Akhmetzianov et al. / IFAC-PapersOnLine 49-12 (2016) 1251–1254

Since the function U does not depend on x, we have (see (2)): l U dx . ∆p = − κ h(s) 0

Then

κ∆p U = − l . dx 0 h(s)

This formula allows us to control the filtration velocity U by a difference of pressure on wells. All the figures in the article were built by I. Boronin, A. Gorinov, and A. Shevlyakove. Authors express gratitude to them. REFERENCES A.V. Akhmetzianov, A.G. Kushner , V.V. Lychagin. Nonlinear Models of Oil Frontal Displacement and Shock Waves. 7th IFAC Conference on Manufacturing Modelling, Management, and Control, 2013. Manufacturing Modelling, Management, and Control, pp. 1170-1175. A.V. Akhmetzianov, A.G. Kushner , V.V. Lychagin. Geometric Theory of Special Modes in the DistributedParameter Control Systems. I. Automation and Remote Control, 11/2013; 74(11): 1786-1801. A.V. Akhmetzianov, A.G. Kushner , V.V. Lychagin. Integrable Models of Oil Displacement. IFAC-PapersOnLine 48-3 (2015) pp. 1264–1267. G.I. Barenblatt, V.M. Entov, V.M Ryzhik. Theory of Fluid Flows Through Natural Rocks. Kluwer, 396 pp., 2010. I.S. Krasil’shchik, V.V. Lychagin, A.M. Vinogradov. Geometry of jet spaces and nonlinear partial differential equations. Advanced Studies in Contemporary Mathematics. – 1. – New York: Gordon and Breach Science Publishers, xx+441 pp., 1986. A.G. Kushner. On contact equivalence of Monge-Ampere equations to linear equations with constant coefficients. Acta Appl. Math. 109(1) P. 197–210. Online First: DOI 10.1007/s10440-009-9447-z (2009). A.G. Kushner, V.V. Lychagin, V.N. Rubtsov. Contact geometry and nonlinear differential equations. Encyclopedia of Mathematics and Its Applications, 101. Cambridge: Cambridge University Press, xxii+496 pp., 2007.

Fig. 4. Evolution of the water saturation. (1) Solve the equation b + φ(−ze−b ) = 0. Let b = B(z) be a solution. (2) Find value a0 = −γeB(0) . (3) Calculate a function g(a) = B −1 (ln(−a/γ)). (4) Calculate a function k(a) by formula (21). (5) Calculate a function r′ (b) = φ−1 (−b). (6) Solve the equation t = k ′ (a) with respect to a. Let a = α(t) be a solution. (7) Solve the equation x = e−b (α(t)k ′ (α(t)) − k(α(t))) + r′ (b) with respect to b. Let b = β(t, x) be a solution. (8) Solve the equation Φ(s) + β(t, x) = 0 with respect to s. Let s = s(t, x) be a solution. (9) Solve the equation h′ (s)sx q + h(s)qx = 0 with condition (5) with respect to q. 5. CONTROL OF FILTRATION VELOCITY Suppose that there is an extracting well at the point x = l and let p0 and pl be pressures on pumped and extracting wells respectively. Difference of pressure on wells is ∆p = pl − p0 =

L

qdx.

0

1254