Application of Fractional Calculus to Oil Industry

Abstract

In this chapter, we present a discussion about the practical application of the fractal properties of the medium in the mathematical model through the use of fractional partial derivatives. We present one of the known models for the flow in saturated media and its generalization in fractional order derivatives. In the middle section, we present one of the main arguments that motivate the use of fractional derivatives in the porous media models, this is the Professor Nigmatullin’s work. The final part describes the process for obtaining the coupled system of three equations for the monophase flow with triple porosity and triple permeability, briefly mentioning the method used for the first solutions of the system.

1. Introduction

The objects of nature rarely have a classical geometric form; in the particular case of oil reservoirs, the ground where the wells are found has been considered with Euclidean geometry; this is not sufficient in many cases to give good approximations in the mathematical models. Since its forms are closer to the fractal geometry, the knowledge of this can be useful to develop models that allow us to better manage the wells. This work presents an approach in fractional derivatives for the triple porosity and triple permeability monophasic saturated model, based on the one proposed by Camacho et al. [1, 2] and generalized partially by Fuentes et al. [3]. The main contribution is to consider the link between fractional equations and fractal geometry through the revision of Alexander-Orbach’s conjecture [23], taken to the particular case.

2. Background of the approach of models of diffusion on fractal media

Fractional calculus was originated as a way to generalize classic calculus; however, it is more difficult to find a direct physical interpretation than in the classical version. When we consider an oil well as a fractal, it is important to choose which of its properties can be useful for elaborating a mathematical model [20, 21, 26, 27].

Alexander and Orbach [4] calculated the “spectral dimension (fracton)”; this parameter is associated to volume and fractal connectivity by being considered as an elastic fractal net of particles connected by harmonic strings. Thus, we consider the particle movement over this fractal and we find a relation of root mean square of an r aleatory walker dependent of time over the fractal, which is in accordance to the following relation:

$⟨r^2(t)⟩\approx t^{2/(2+\theta )},$E1

where r is in euclidean space. Alexander and Orbach defined $d_s=\frac{2d_f}{2+\theta }$ as the spectral dimension or fracton, where $d_w=2+\theta$ is the dimension of the walk, θ gives us the dependence of the diffusion constant over the distance and d_f is the effective dimension [24, 25].

O’Shaughnessy and Procaccia [5] used the concepts of Alexander and Orbach to formulate their fractal diffusion equation:

$\frac{\partial p(r,t)}{\partial t}=\frac{1}{r^d{f}^{-1}}\frac{\partial }{\partial r}\left[K_1{r}^{-\theta }{r}^{d_f-1}\frac{\partial p(r,t)}{\partial r}\right],r>\mathrm{0, }{K}_1\text{constant,}$E2

with solution.

$P(r,t)=\frac{2+\theta }{d_f\text{Γ}(d_f/(2+\theta ))}{\left[\frac{1}{{(2+\theta )}^2{K}_{1}t}\right]}^{d_f/(2+\theta )}\exp \left[-\frac{r^{2+\theta }}{{(2+\theta )}^2{K}_{1}t}\right],$E3

of which one finds a power law

$⟨r^2(t)⟩=\frac{\text{Γ}\left[\frac{d_f+2}{2+\theta }\right]}{\text{Γ}\left[\frac{d_f}{2+\theta }\right]}{\left[{(2+\theta )}^2{K}_{1}t\right]}^{2/(2+\theta )}=⟨r^2(1)⟩r^{2/(2+\theta )}\mathrm{.}$E4

Metzler et al. [6] started with the characterization of an anomalous diffusion process 1. Here, they consider $d_w=\theta +2$ as the anomalous diffusion exponent, they are referencing the work of Havlin and Ben-Avraham [7] to calculate diffusion with a media (1). They obtain an “approach exponential”:

$P(r,t)\approx A{t}^{-d_f/d_w}{e}^{-c{(r/R)}^u},$E5

valid in the asyntotic range $r/R\gg 1$, and t → ∞ with R and r defined by

$\{\begin{array}{l}R=\sqrt{⟨r^2(t)⟩}\hfill \\ u=\frac{d_w}{d_w-1}\hfill \end{array}.$E6

Thus, it is possible to obtain the solution of the fractional derivative diffusion equation:

$\frac{{\partial }^{2/d_w}}{\partial t^{2/d_w}}P(r,t)=\frac{1}{r^{d_s-1}}\frac{\partial }{\partial r}\left(r^{d_s-1}\frac{\partial }{\partial r}P(r,t)\right)$E7

where

$\frac{{\partial }^{2/d_w}}{\partial t^{2/d_w}}P(r,t)=\frac{1}{\text{Γ}\left(1-\frac{2}{d_w}\right)}\frac{\partial }{\partial t}{\displaystyle {\int }_0^t}d\tau \frac{P(r,\tau )}{{(t-\tau )}^{2/d_w}},0\le \frac{2}{d_w}<1.$E8

3. Brief history of fractional calculus

In mathematics, one way to obtain new concept is to generalize by extending one definition or context for values not previously considered. For example, it is possible to generalize the power concept of xⁿ, for natural n values such as the concept of x, n times, to negative integers n, as the product of $\frac{1}{x}$, n times, then to n rational values such as $\sqrt[q]{x^p}$, if $n=\frac{p}{q}$, with positive p and q. In each step, the generalization modifies the concept a little, but it keeps the previous one as a particular case. This process can continue all the way to a complex n. In the same way as generalizations in differential and integral calculus have been made, in this case the generalization goes toward the n order of the $\frac{d^{n}y}{d{x}^n}$ derivative [22].

Leibniz: In a letter dated September 30, 1695, L’Hôpital, he has been inquired about the meaning of $\frac{d^{n}y}{d{x}^n}$, if $n=\frac{1}{2}$, in response he wrote: “You can see by that, sir, that one can express by an infinite series a quantity such as $d^{1/2}xy$ or $d^{1:2}xy$. Although infinite series and geometry are distant relations, infinite series admits only the use of exponents that are positive and negative integers, and does not, as yet, know the use of fractional exponents.” Later in the same letter, Leibniz continues: “Thus, it follows that $d^{1/2}x$ will be equal to $x\sqrt{dx:x}$. This is an apparent paradox from which, one day, useful consequences will be drawn.”

In his correspondence with Johan Bernoulli, Leibniz mentioned to him general order derivatives. In 1697, he established that differential calculus can be used to achieve these generalizations and used the $d^{1/2}$ notation to denote order $\frac{1}{2}$ derivative [22].

Euler: In 1730, Euler proposed derivatives as rate between functions and variables that can be expressed algebraically; the solution with this approach when the orders are not integers is the use of interpolations.

The (fractional) non-integer order derivative motivated Euler to introduce the Gamma function. Euler knew that he needed to generalize (or, as he said, interpolate), the $1\cdot 2\cdots n=n!$ product for non-integer n. He proposed an integral:

${\displaystyle {\displaystyle \prod _{k=1}^n}}n={\displaystyle {\int }_0^1}{(-\log x)}^{n}dx,$E9

and used it to partially solve the Leibniz paradox. He also gave the basic fractional derivative (with modern notation $\text{Γ}(n+1)=n!$):

$\frac{d^{\alpha }{x}^{\beta }}{d{x}^{\alpha }}=\frac{\text{Γ}\left(\beta +1\right)}{\text{Γ}\left(\beta -\alpha +1\right)}{x}^{\beta -\alpha },$E10

which is valid for non-integer α and β [22].

Laplace and Lacroix: Laplace also defined his fractional derivative via an integral. In 1819, Lacroix, applying the (10) formula and the Legendre symbol for Gamma function, was able to calculate the derivative with y = x and $n=\frac{1}{2}$. He was also the first to use the term “fractional derivative.” He thus achieved

$\frac{d^{1/2}y}{d{x}^{1/2}}=\frac{2\sqrt{x}}{\sqrt{\pi }}\mathrm{.}$E11

Fourier: Joseph Fourier (1822), in his famous book “The Analytical Theory of Heat” making use of this expression of a function and an interpretation of the sines and cosines derivatives gave his definition of a fractional derivative:

$f(x)=\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }}f(\alpha )d\alpha {\displaystyle {\int }_{-\infty }^{\infty }}\cos p(x-\alpha )dp,$E12

then

$\frac{d^n}{d{x}^n}\cos p(x-a)=p^n\cos \left[p(x-\alpha )+\frac{1}{2}n\pi \right],$E13

for an integer n. Formally replacing n with an arbitrary u, he obtained the generalization:

$\frac{d^u}{d{x}^u}f(x)=\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }}f(\alpha )d\alpha {\displaystyle {\int }_{-\infty }^{\infty }}{p}^u\cos \left[p(x-a)+\frac{1}{2}u\pi \right]dp\mathrm{.}$E14

Fourier thus establishes that the u number can be regarded as any quantity, positive or negative [8, 22].

Abel: In 1823, N. H. Abel published the solution of a problem presented by Hyugens in 1673: The tautochrone problem. Abel gave his solution in the form of an integral equation that is considered the first application of fractional calculus. The integral he worked with is

${\displaystyle {\int }_0^x}{(x-t)}^{-\frac{1}{2}}f(t)dt=k\mathrm{.}$E15

This integral is, except for the $1/\text{Γ}\left(\frac{1}{2}\right)$ factor, a fractional integral of 1/2 order, Abel wrote the left part as $\pi \left[\frac{d^{-\frac{1}{2}}}{d{x}^{-\frac{1}{2}}}\right]f(x)$, thus he worked with both sides of the equation as

$\sqrt{\pi }f(x)=\frac{d^{\frac{1}{2}}}{d{x}^{\frac{1}{2}}}k\mathrm{.}$E16

The first integral equation in history had been solved. Two facts may be observed: the regard for the sum of the orders, and that unlike in classical calculus, the derivative of a constant is not zero [8, 22].

Liouville: In 1832, Liouville made the first great study of fractional calculus. In his work, he considered $\left(\frac{d^{1/2}}{d{x}^{1/2}}\right)e^{2x}$. The first formula he obtained was the derivative of a function:

$f(x)={\displaystyle {\displaystyle \sum _{n=0}^{\infty }}}{c}_n{e}^{a_{n}x},\mathfrak{R}(a_n)>\mathrm{0,}$E17

from which he got

$D^{\nu }f(x)={\displaystyle {\displaystyle \sum _{n=0}^{\infty }}}{x}_n{a}_n^{\nu }{e}^{a_{n}x}\mathrm{.}$E18

that can be obtained using the extension

$D^{\nu }{e}^{ax}=a^{\nu }{e}^{ax},$E19

for an arbitrary number ν. A second definition was achieved by Liouville from the defined integral:

$I={\displaystyle {\int }_0^{\infty }}{u}^{a-1}{e}^{-xu}du,a>\mathrm{0,}x>\mathrm{0,}$E20

of which, after a change of variable and a suitable rewriting is obtained

$D^{-\nu }{x}^{-a}=\frac{{(-1)}^{\nu }\text{Γ}\left(a+\nu \right)}{\text{Γ}(a)}{x}^{-a-\nu },a>0.$E21

Liouville also tackled the tautochrone problem and proposed differential equations of arbitrary order.

In 1832, he wrote about a generalization of Leibnitz’s rule about the nth derivative of a product:

$D^{\nu }f(x)g(x)={\displaystyle {\displaystyle \sum _{n=0}^{\infty }}}\left(\begin{array}{c}\nu \\ n\end{array}\right)D^{n}f(x)D^{\nu -n}g(x),$E22

where Dⁿ is the ordinary n order differential operator, D^ν−n fractional operator, and $\left(\begin{array}{c}\nu \\ n\end{array}\right)$ the generalized binomial coefficient, expressed in terms of the Gamma function, $\frac{\text{Γ}(\nu +1)}{n!\text{Γ}(\nu -n+1)}$.

Liouville expanded the coefficients in Eq. (18) as

$\begin{array}{c}{a}_n^{\nu }={\displaystyle \underset{h\to 0}{{\displaystyle \lim }}}\frac{1}{h^{\nu }}{\left(1-e^{-h{a}_n}\right)}^{\nu },a_n>\mathrm{0,}\\ (-1{)}^{\nu }{a}_n^{\nu }={\displaystyle \underset{h\to 0}{{\displaystyle \lim }}}\frac{1}{h^{\nu }}{\left(1-e^{h{a}_n}\right)}^{\nu },a_n<0.\end{array}$E23

And inserted those equations in Eq. (18) to get

$\begin{array}{l}\frac{d^{\nu }}{d{x}^{\nu }}f(x)={\displaystyle \underset{h\to 0}{{\displaystyle \lim }}}\left\{\frac{1}{h^{\nu }}{\displaystyle {\displaystyle \sum _{n=0}^{\infty }}}\left[{(-1)}^m\left(\begin{array}{c}\nu \\ n\end{array}\right)f(x-mh)\right]\right\},\\ =(-1{)}^{\nu }{\displaystyle \underset{h\to 0}{{\displaystyle \lim }}}\left\{\frac{1}{h^{\nu }}{\displaystyle {\displaystyle \sum _{n=0}^{\infty }}}\left[{(-1)}^m\left(\begin{array}{c}\nu \\ n\end{array}\right)f(x+mh)\right]\right\}\mathrm{.}\end{array}$E24

These formulas would be retaken by Grünwald in 1867.

Riemann: Riemann developed his Fractional Calculus theory when he was preparing his Ph.D. thesis, but his oeuvre was published posthumously around 1892. He searched for a generalization of Taylor’s series, in which he defined the n-th differential coefficient of a f(x) function as the hⁿ coefficient in the f(x + h) expansion with integer powers of h. Thus, he generalizes this definition to non-integer powers and demands that

$f(x+h)={\displaystyle {\displaystyle \sum _{n=-\infty }^{n=\infty }}}c-n+\alpha ({\partial }_x^{n+\tau }f)(x)h^{n+\alpha },$E25

be valid for $n\in \mathbb{N},a\in \mathbb{R}$. The $c_{n+\alpha }$ factor is determined by the ${\partial }^{\beta }({\partial }^{\alpha }f)={\partial }^{\beta +\alpha }f$ condition, and he found that it was $\frac{1}{\text{Γ}(n+\alpha +1)}$. Riemann then derived Eq. (25) expression for negative α:

${\partial }^{\alpha }f=\frac{1}{\text{Γ}(-\alpha )}{\displaystyle {\int }_k^x}{(x-t)}^{-\alpha -1}f(t)dt+{\displaystyle {\displaystyle \sum _{n=1}^{\infty }}}{K}_n\frac{x^{-\alpha -n}}{\text{Γ}\left(-n-\alpha +1\right)},$E26

where k, K_n are finite constants. Then, he extended the result to non-negative α.

Sonin and Letnikov: The Russian mathematicians N.Ya. Sonin (1868) and A.V. Letnikov (1868–1872) [29] made contributions taking as basis the formula for the nth derivative of the Cauchy integral formula given by

$D^{n}f(z)=\frac{n!}{2\pi i}{\displaystyle {\int }_C}\frac{f(\xi )}{{(\xi -z)}^{n+1}}d\xi \mathrm{.}$E27

They worked using the contour integral method, with the contribution of Laurent (1884), they achieved the definition:

${}_c{D}_x^{-\alpha }f(x)=\frac{1}{\text{Γ}\left(\alpha \right)}{\displaystyle {\int }_c^x}{(x-t)}^{\alpha -1}f(t)dt,\mathfrak{R}(\alpha )>0.$E28

For an integration to an arbitrary order, when x > c has the Riemann definition, but without a complementary definition, when c = 0 we get the shape known as Riemann-Liouville fractional integral:

${}_0{D}_x^{-\alpha }f(x)=\frac{1}{\text{Γ}(\alpha )}{\displaystyle {\int }_0^x}{(x-t)}^{\alpha -1}f(t)dt,\mathfrak{R}(\alpha )>0.$E29

Assigning c values in Eq. (19), we get different integrals of fractional order, which will be fundamental to define fractional derivatives.

If c = −∞, we get

${}_{-\infty }{D}_x^{-\alpha }f(x)=\frac{1}{\text{Γ}(\alpha )}\underset{-\infty }{\overset{x}{{\displaystyle \int }}}{(x-t)}^{\alpha -1}f(t)dt,\Re (\alpha )>0.$E30

Using integration properties, more definitions will be given.

Grünwald: Another contribution is that of Grünwald (1867) and Letnikov (1868). This extension of the classical derivative to fractional order is important because it lets us apply it in numerical approximations. They started with the definition of derivative as a limit given by Cauchy (1823):

$\frac{df}{dx}={\displaystyle \underset{h\to \infty }{{\displaystyle \lim }}}\frac{\left[f(x)-f(x-h)\right]}{h}\mathrm{.}$E31

First generalizing for a nth integer derivative we get

$\begin{array}{l}{D}^{n}f(x)=\underset{h\to 0}{\lim }\frac{{\displaystyle \sum }_{j=0}^n[{(-1)}^j\left(\begin{array}{c}n\\ j\end{array}\right)f(x-jh)]}{h^n},\hfill \\ n\in \mathbb{N},\text{and}f\in C^n[a,b],a<x<b.\hfill \end{array}$E32

Grünwald generalizes Eq. (32) for an arbitrary q value, expressing it as

$D_a^{q}f(x)={\displaystyle \underset{N\to \infty }{{\displaystyle \lim }}}{h}_N^{-q}\left[{\displaystyle {\displaystyle \sum _{j=0}^N}}{(-1)}^j\left(\begin{array}{c}-n\\ j\end{array}\right)f(x-j{h}_N)\right],q\in \mathbb{R},$E33

where the binomial coefficient is

$\left(\begin{array}{c}q\\ j\end{array}\right)=\frac{q(q-1)(q-2)\cdots (q-j+1)}{j!},$E34

also showing that

${(-1)}^j\left(\begin{array}{c}q\\ j\end{array}\right)=\left(\begin{array}{c}j-q-1\\ j\end{array}\right)=\frac{\text{Γ}(j-q)}{\text{Γ}(-q)\text{Γ}(j+1)}\mathrm{.}$E35

and with those previous results, it is possible to establish this important property for $\alpha \in \mathbb{R}$ and $n\in \mathbb{N}$

$\frac{d^n}{d{x}^n}{D}_a^{\alpha }f(x)=D_a^{n+\alpha }f(x).$E36

In the twentieth and twenty-first centuries, more definitions will rise, but they will be given in terms within the Riemann-Liouville fractional integral and will be part of the Modern Fractional Calculus Theory, in all their fundamental definitions [22].

4. Fractional calculus

We will now present the assorted definitions and notations of fractional derivatives that will be used throughout this work. It is worth pointing out that this is necessary because such notation is currently standardized [18, 19].

4.1. Riemann-Liouville fractional derivative

The Riemann-Liouville derivative is the basis to define most fractional derivatives; it generalizes the Cauchy’s formula for derivatives of high order. For an f function defined in a [a, b] interval, a $\alpha \in \mathbb{C}$ value with $\mathfrak{R}(\alpha )>0$ defines the left and right Riemann-Liouville integrals by

$\left({}_a^{RL}{I}_x^{\alpha }f\right)(x)=\frac{1}{\text{Γ}(\alpha )}{\displaystyle {\int }_a^x}\frac{f(t)}{{(x-t)}^{1-\alpha }}dt,\left(x>a\right),$E37

$\left({}_x^{RL}{I}_b^{\alpha }f\right)(x)=\frac{1}{\text{Γ}(\alpha )}{\displaystyle {\int }_x^b}\frac{f(t)}{{(t-x)}^{1-\alpha }}dt,\left(x<b\right)\mathrm{.}$E38

Following Riemann’s notion of defining fractional derivatives as the integer order derivative of an fractional integral, we have the left and right derivative proposal as follows:

${(}_a^{RL}{D}_x^{\alpha }f)(x)={\left(\frac{d}{dx}\right)}^n\left({(}_a^{RL}{I}_x^{n-\alpha }f)(x)\right),x>a,$E39

${(}_x^{RL}{D}_b^{\alpha }f)(x)={\left(-\frac{d}{dx}\right)}^n\left({(}_b^{RL}{I}_x^{n-\alpha }f)(x)\right),x<b,$E40

with $n=-\lfloor \mathfrak{R}(\alpha )\rfloor$, i.e., $n=\lfloor \mathfrak{R}(\alpha )\rfloor +1$.

As shown in Refs. [8–10], these operators generalize the usual derivation. In other words, when $\alpha \in {\mathbb{N}}_0$, then

$\begin{array}{l}{(}_a^{RL}{D}_x^{0}f)(x)={(}_x^{RL}{D}_b^{0}f)(x)=f(x),\text{if}\alpha =\mathrm{0,}\end{array}$E41

$\begin{array}{l}{(}_a^{RL}{D}_x^{n}f)(x)=f^{(n)}(x),{(}_x^{RL}{D}_b^{n}f)(x)=(-1{)}^n{f}^{(n)}(x).\end{array}$E42

It is also possible to prove that the semigroup propriety about the order of integral operators (i.e., for $\alpha \in \mathbb{C}(\mathfrak{R}(\alpha >0)),\beta \in \mathbb{C}(\mathfrak{R}(\beta )>0)$) is achieved:

$\left({}_a^{RL}{I}_x^{\alpha }\left({}_a^{RL}{I}_x^{\beta }f\right)\right)(x)=\left({}_a^{RL}{I}_x^{\alpha +\beta }f\right)(x),$E43

$\left({}_x^{RL}{I}_b^{\alpha }\left({}_x^{RL}{I}_b^{\beta }f\right)\right)(x)=\left({}_x^{RL}{I}_b^{\alpha +\beta }f\right)(x).$E44

For the derivatives, we have

$\left({}_a^{RL}{D}_x^{\alpha }\left({}_a^{RL}{D}_x^{\beta }f\right)\right)(x)=\left({}_a^{RL}{D}_x^{\alpha +\beta }f\right)(x)-{\displaystyle {\displaystyle \sum _{j=1}^m}}\left({}_a^{RL}{D}_x^{\beta -j}f\right)(a)\frac{{(x-a)}^{-j-\alpha }}{\text{Γ}(1-j-\alpha )}\mathrm{.}$E45

For $f(x)\in L_p(1\le p\le \infty )$, the following relationships are valid:

$\left({}_a^{RL}{D}_x^{\beta }\left({}_a^{RL}{I}_x^{\alpha }f\right)\right)(x{)=}_a^{RL}{I}_x^{\alpha -\beta }f(x),$E46

$\left({}_x^{RL}{D}_b^{\beta }\left({}_x^{RL}{I}_b^{\alpha }f\right)\right)(x{)=}_x^{RL}{I}_b^{\alpha -\beta }f(x).$E47

If α = β_, we have the identity operator and the operators turn out to be inverted. On the other hand, if the order of the operators is inverted, it will have

$\left({}_a^{RL}{I}_x^{\beta }\left({}_a^{RL}{D}_x^{\alpha }f\right)\right)(x)=f(x)-{\displaystyle {\displaystyle \sum _{j=1}^n}}\frac{f_{n-\alpha }^{(n-j)}(a)}{\text{Γ}(\alpha -j+1)}{(x-\alpha )}^{\alpha -j},$E48

where $\mathfrak{R}(\alpha )>\mathrm{0,}n=\lfloor \mathfrak{R}(\alpha )\rfloor +1$ and $f_{n-\alpha }(x)=\left({}_a^{RL}{I}_x^{n-\alpha }f\right)(x)$ in analogy for the right derivative.

All these properties can be used in the phenomena modeling and its solution; such models have shown to improve usual approaches. However, when using equations with Riemann-Liouville type fractional derivatives, the initial conditions cannot be interpreted physically; a clear example is that the derivative Riemann-Liouville of a constant is not zero, contrary to the impression that the derivatives gives a notion about the change that the function experiences when advancing in the time or to modify its position. This was the motivation for another definition that is better coupled with physical interpretations; this is the derivative of Caputo type.

4.2. Caputo fractional derivative

Michele Caputo [11] published a book in which he introduced a new derivative, which had been independently discovered by Gerasimov (1948). This derivative is quite important, because it allows for understanding initial conditions, and is used to model fractional time. In some texts, it is known as the Gerasimov-Caputo derivative.

Let [a, b] be a finite interval of the real line $\mathbb{R}$, for $\alpha \in \mathbb{C}(\mathfrak{R}(\alpha )\ge 0)$. The left and right Caputo derivatives are defined as

$\left({}_a^C{D}_x^{\alpha }y\right)(x)=\frac{1}{\text{Γ}(n-\alpha )}{\displaystyle {\int }_a^x}\frac{y^{(n)}dt}{{(x-t)}^{\alpha -n+1}}=\left({}_a^{RL}{I}_x^{n-\alpha }{D}^{n}y\right)(x),$E49

$\left({}_x^C{D}_b^{\alpha }y\right)(x)=\frac{{(-1)}^n}{\text{Γ}(n-\alpha )}{\displaystyle {\int }_x^b}\frac{y^{(n)}dt}{{(t-x)}^{\alpha -n+1}}=(-{1)}^n\left({}_x^{RL}{I}_b^{n-\alpha }{D}^{n}y\right)(x),$E50

where $D=\frac{d}{dx}$ and $n=--\mathfrak{R}(\alpha )$, i.e., $n=\mathfrak{R}(\alpha )+1$ for $\alpha \notin {\mathbb{N}}_0$ and $n=\alpha$ for $\alpha \in {\mathbb{N}}_0$. And if $0<\mathfrak{R}(\alpha )<1$

$\begin{array}{l}\left({}_a^C{D}_x^{\alpha }y\right)(x)=\frac{1}{\text{Γ}(1-\alpha )}{\displaystyle {\int }_a^x}\frac{y^{\prime }dt}{{(x-t)}^{\alpha }}=\left({}_a^{RL}{I}_x^{1-\alpha }Dy\right)(x),\end{array}$E51

$\begin{array}{l}\left({}_x^C{D}_b^{\alpha }y\right)(x)=-\frac{1}{\text{Γ}(n-\alpha )}{\displaystyle {\int }_x^b}\frac{y^{\prime }dt}{{(t-x)}^{\alpha }}=-\left({}_x^{RL}{I}_b^{1-\alpha }Dy\right)(x).\end{array}$E52

The connection between Caputo and Riemann derivatives is given by the relations

$\left({}_a^C{D}_x^{\alpha }y\right)(x)=\left({}_a^{RL}{D}_x^{\alpha }\left[y(t)-{\displaystyle {\displaystyle \sum _{k=0}^{n-1}}}\frac{y^{(k)}(a)}{k!}{(t-a)}^k\right]\right)(x),$E53

$\left({}_x^C{D}_b^{\alpha }y\right)(x)=\left({}_x^{RL}{D}_b^{\alpha }\left[y(t)-{\displaystyle {\displaystyle \sum _{k=0}^{n-1}}}\frac{y^{(k)}(b)}{k!}{(b-t)}^k\right]\right)(x).$E54

In particular, if $0<\mathfrak{R}(\alpha )<1$, Eqs. (53) and (54) relation take the following shapes:

$\left({}_a^C{D}_x^{\alpha }y\right)(x)=\left({}_a^{RL}{D}_x^{\alpha }[y(t)-y(a)]\right)(x),$E55

$\left({}_x^C{D}_b^{\alpha }y\right)(x)=\left({}_x^{RL}{D}_b^{\alpha }[y(t)-y(b)]\right)(x).$E56

For α = n, then the Caputo derivatives match classical derivatives except for the sign of the right derivative.

However, for $k=\mathrm{0,1,}\dots ,n-1$, we have

$\left({}_a^C{D}_x^{\alpha }{(t-a)}^k\right)(x)=\mathrm{0,}\left({}_x^C{D}_b^{\alpha }{(b-t)}^k\right)(x)=\mathrm{0,}$E57

in particular,

$\left({}_a^C{D}_x^{\alpha }1\right)(x)=\mathrm{0,}\left({}_x^C{D}_b^{\alpha }1\right)(x)=0.$E58

On the other hand, if $\mathfrak{R}(\alpha )>0$ and λ > 0, then

$\left({}_a^C{D}_x^{\alpha }{e}^{\lambda t}\right)(x)\ne {\lambda }^{\alpha }{e}^{\lambda t},\text{for}\alpha \in \mathbb{R}\mathrm{.}$E59

The Caputo derivatives behave like inverted operators for the left Riemann-Liouville fractional integrals ${}_a^{RL}{I}_x^{\alpha }$ and ${}_x^{RL}{I}_b^{\alpha }$, if $\mathfrak{R}(\alpha )>0$ and $y(x)\in C[a,b]$

$\begin{array}{cc}\left({}_a^C{D}_x^{\alpha }\left({}_a^{RL}{I}_x^{\alpha }y\right)\right)\left(x\right)=y\left(x\right),& \left({}_x^C{D}_b^{\alpha }\left({}_x^{RL}{I}_b^{\alpha }y\right)\right)(x)=y(x).\end{array}$E60

On the other hand, if $\mathfrak{R}(\alpha )>0$ and $n=-\lfloor -\mathfrak{R}(\alpha )\rfloor$, then for good conditions for y(x)

$\begin{array}{l}\left({}_a^{RL}{I}_x^{\alpha }\left({}_a^C{D}_x^{\alpha }y\right)\right)(x)=y(x)-{\displaystyle {\displaystyle \sum _{k=0}^{n-1}}}\frac{y^{(k)}(a)}{k!}(x-a{)}^k,\end{array}$E61

$\begin{array}{l}\left({}_x^{RL}{I}_b^{\alpha }\left({}_x^C{D}_b^{\alpha }y\right)\right)(x)=y(x)-{\displaystyle {\displaystyle \sum _{k=0}^{n-1}}}\frac{{(-1)}^k{y}^{(k)}(b)}{k!}(b-x{)}^k\mathrm{.}\end{array}$E62

In particular if, $0<\mathfrak{R}(\alpha )\le 1$, then

$\left({}_a^{RL}{I}_x^{\alpha }\left({}_a^C{D}_x^{\alpha }y\right)\right)(x)=y(x)-y(a),$E63

$\left({}_x^{RL}{I}_b^{\alpha }\left({}_x^C{D}_b^{\alpha }y\right)\right)(x)=y(x)-y(b).$E64

In his early articles and several after that, Caputo used a Laplace transformed of the Caputo fractional derivative, which is given by

$\left(\mathcal{L}{\{}_0^C{D}_x^{\alpha }y\}\right)(s)=s^{\alpha }(\mathcal{L}y)(s)-{\displaystyle {\displaystyle \sum _{k=0}^{n-1}}}{s}^{\alpha -k-1}(D^{k}y)(0).$E65

When $0<\alpha \le 1$, then

$\left(\mathcal{L}\left\{{}_0^C{D}_x^{\alpha }y\right\}\right)(s)=s^{\alpha }(\mathcal{L}y)(s)-s^{\alpha -1}y(0).$E66

These derivatives can be defined over the whole real axis resulting in the expressions:

$\left({}_{}^C{D}_x^{\alpha }y\right)(x)=\frac{1}{\text{Γ}(n-\alpha )}{\displaystyle {\int }_{-\infty }^x}\frac{y^{(n)}(t)dt}{{(x-t)}^{\alpha -n+1}},$E67

$\left({}_x^C{D}^{\alpha }y\right)(x)=\frac{{(-1)}^n}{\text{Γ}(n-\alpha )}{\displaystyle {\int }_x^{\infty }}\frac{y^{(n)}(t)dt}{{(t-x)}^{\alpha -n+1}},$E68

with $x\in \mathbb{R}$.

5. Fractal geometry and fractional calculus

The phenomenon of anomalous diffusion is mathematically modeled by a fractional partial differential equation. The parameters of this equation are uniquely determined by the fractal dimension of the underlying object.

There are some results that show the relationship between fractals and fractional operators [24]; two of the most important that motivated the particular study of the equations to determine the pressure deficit in oil wells are highlighted below.

5.1. Cantor’s Bars and fractional integral

In 1992, Nigmatullin [12] presents one of the most distinguished contributions to the search of the concrete relationship between the fractal dimension of a porous medium and the order of the fractional derivative to model a phenomena through such a medium; in this, he achieves the evolution of a physical system of a Cantor set type.

In his research, Nigmatullin proposes a relationship between the fractal dimension of a Cantor type set and the order of a fractional integral of the Riemann-Liouville type. The systems he considers are named phenomena with “memory.” The use of fractional derivatives given by assuming a transference function J(t) in relationship to a rectifiable function f(t) through the convolution operator with a distribution K*(t) establishes that

$J(t)=K^{\star }(t)*f(t)={\displaystyle {\int }_0^t}{K}^{\star }(t-\tau )f(\tau )d\tau \mathrm{.}$E69

Where the distribution to apply (see Refs. [13, 14]) is a so-called “Cantor’s Bars” $K_{T,\nu }(t)$, supported in the [0, T] interval, with a fractal dimension $\nu =\ln (2)/\ln (1/\xi )$, with $\xi \in [\mathrm{0,1}/2]$ being the compression factor, normalized in L¹.

Through the result of distribution values, he establishes the relation:

$\begin{array}{l}J(t)=⟨K_{T,\nu }(t)⟩*f(t)=B(\nu )T^{-\nu }{}_0^{RL}{D}_x^{-\nu }\left[f(t)\right]=\frac{B(\nu )T^{-\nu }}{\text{Γ}(\nu )}\left(t^{\nu -1}*f(t)\right),\end{array}$E70

$⟨K_{T,\nu }(t)⟩={\displaystyle {\int }_{-1/2}^{1/2}}{K}_{T,\nu }(t{\xi }^{-x})(2\xi {)}^{-x}dx=\frac{B(\nu )}{\text{Γ}(\nu )}{\left(\frac{t}{T}\right)}^{\nu -1},$E71

$B(\nu )={\displaystyle {\int }_{-1/2}^{1/2}}{q}_{\nu }(z+x\ln \xi )dx\mathrm{.}$E72

Thus, assuming a porous medium with a ν fractal dimension, we establish a fractional derivative of −ν order.

The initial results were strongly questioned by different authors, including Roman Rutman (see Refs. [15, 16]), who asserts that the relation is so artificial. However, recent works suggest that Nigmatullin’s statements are not far from reality, but it is necessary to reduce the set of functions and that of fractals for which the necessary convergence is fulfilled.

6. Fractional calculus for modeling oil pressure

In this section, the Equation Continuity which follows from the law of conservation of mass is established. Darcy’s law is used to relate fluid motion to pressure and gravitational gradients. The combination of the Continuity Equation and Darcy’s Law leads to a heat-conducting differential equation in mathematical physics describing the transfer of the fluid. We obtain a system formed by three partial differential equations, one for each fluid. This multiphase system must be solved considering the relevant boundary and initial conditions [30].

In the particular case of naturally fractured reservoirs (see Refs. [1, 2]), usually it is possible to discern three porosity types: matrix, fracture, and vugs; with this conception, it is accepted that the three porosities have associated a solid phase, and with this both Continuity Equation and Darcy’s law can be expressed for each fluid in each geometrical media. If we only consider oil (monophasic) in a isotropic and saturated media, we can obtain a three equations system; for this, we begin with standard continuity equation and standard Darcy’s law, respectively (see Ref. [17]):

$\frac{\partial (\rho \theta )}{\partial t}+\nabla \cdot p(\rho q)=\rho \Upsilon ,q=-\frac{1}{\mu }k(p)(\nabla p-\rho g\nabla D),$E73

where θ is the volumetric content of fluid; $q=(q_1,q_2,q_3)$ is the Darcy flux, with its spatial components (x, y, z), t is the time; ρ is the density of the fluid; μ is the dynamic viscosity of the fluid; g gravitational acceleration; $\Upsilon$ is a source term and represents a volume provided per volume unit of porous media in the time unity; p is the pressure; D is the depth as a function of spatial coordinates, usually identified to the vertical coordinate z; k is the permeability tensor of the partially saturated porous media and it depends on the pressure. The relations θ(p) and k(p) are the fluid-dynamics characteristics of the media.

General fluid transfer equation results combining the formulas in Eq. (73):

$\frac{\partial (\rho \theta )}{\partial t}=\nabla \cdot p[\frac{\rho }{\mu }k(p)(\nabla p-\rho g\nabla D)]+\rho \Upsilon \mathrm{.}$E74

This differential equation contains two dependent variables, namely the humidity content and fluid pressure, but they are related. For this reason, the saturation S(p) is defined so that

$\theta (p)=\varphi (p)S(p)$E75

where ϕ is the total porosity of the medium, and the specific capacity defined by

$C(p)=\frac{d(\rho \varphi S)}{dp}=\varphi S\frac{d\rho }{dp}+\rho S\frac{d\varphi }{dp}+\rho \varphi \frac{dS}{dp},$E76

in consequence

$\frac{\partial (\rho \theta )}{\partial t}=C(p)\frac{\partial p}{\partial t},$E77

6.1. Triadic media

The porous media is considered to be formed by three porous media: the matrix, fractured media, and vuggy media. The total volume of the porous media (V_T) is equal to the sum of the total volume of the matrix (V_m), of the total volume of the fractured medium (V_F) and of the total volume of the vuggy media (V_G). In other words

$V_T=V_M+V_F+V_G,$E78

each of the porous media contains solids and voids so that

$\begin{array}{c}{V}_M=V_{MS}+V_{MV}\end{array}$E79

$V_F=V_{FS}+V_{FV}$E80

$V_G=V_{GS}+V_{GV}.$E81

The porous medium as everything contains solids and voids, with the following relations:

$V_T=V_{TS}+V_{TV},V_{TS}=V_{MS}+V_{FS}+V_{GS},V_{TV}=V_{MV}+V_{FV}+V_{GV}\mathrm{.}$E82

The volume fraction occupied by the matrix is defined as (ν_M), the volume fraction occupied by the fractured media as ν_F, and the fraction that occupies the vuggy media as (ν_G) relative to the total volume of the porous medium given by

${\nu }_M=\frac{V_M}{V_T},{\nu }_f=\frac{V_F}{V_T},{\nu }_G=\frac{V_G}{V_T},{\nu }_M+{\nu }_F+{\nu }_G=1.$E83

The porosity of the porous media (ϕ), in matrix (ϕ_M), fracture media (ϕ_F) and vuggy media (ϕ_G) are defined by

$\varphi =\frac{V_{TV}}{V_T},{\varphi }_M=\frac{V_{MV}}{V_M},{\varphi }_F=\frac{V_{FV}}{V_F},{\varphi }_G=\frac{V_{GV}}{V_G}\mathrm{.}$E84

From the above equations, we deduce the relation between the porosities:

$\varphi ={\nu }_M{\varphi }_M+{\nu }_F{\varphi }_F+{\nu }_G{\varphi }_G$E85

When the empty space contains fluid partially, the total volumetric content of the fluid (θ) as the total fluid volume (V_TW) with respect to the total volume of the porous medium is $(V_T):\theta =V_{TW}/V_T$. In an analogous way, the volumetric content of fluid in the matrix is defined ${\theta }_M=V_{MW}/V_M$, in the fractured media ${\theta }_F=V_{FW}/V_F$ and vuggy media ${\theta }_G=V_{GW}/V_G$. It follows that

$\theta ={\nu }_M{\theta }_M+{\nu }_F{\theta }_F+{\nu }_G{\theta }_G$E86

which is reduced to Eq. (85) when the porous medium is fully saturated with fluid. It is satisfied: $0<\theta <\varphi \mathrm{,0}<{\theta }_M<{\varphi }_M\mathrm{,0}<{\theta }_F<{\varphi }_F\mathrm{,0}<{\theta }_G<{\varphi }_G$. The relation between the total volumetric flow of the fluid per unit area in the porous medium (q), the volumetric flow per unit area in the matrix (q_M), the volumetric flow per unit area in the fractured medium q_F, and the volumetric flow per unit area in the vuggy media (q_G) is analogous to Eq. (86), namely

$q={\nu }_M{q}_M+{\nu }_F{q}_F+{\nu }_G{q}_G$E87

The continuity equations for the matrix, the fractured medium, and the vuggy media considering Eq. (86) acquire the form

$\frac{\partial (\rho {\theta }_i)}{\partial t}+\nabla \cdot p(\rho q_i)=\rho {\mathrm{\Upsilon }}_i,i=M,F,G\mathrm{.}$E88

Darcy’s law for the matrix, the fractured medium, and the vuggy media, takes the form

$q_i=-\frac{1}{\mu }{k}_i(p_i)(\nabla p_i-\rho g\nabla D),i=M,F,G$E89

The equation of continuity of the porous medium, Eq. (73), is deduced from the sum of Eq. (88) previously multiplied by ${\nu }_m,{\nu }_F,{\nu }_G$, respectively, if the source terms are related by

$\mathrm{\Upsilon }={\nu }_M{\mathrm{\Upsilon }}_M+{\nu }_F{\mathrm{\Upsilon }}_F+{\nu }_G{\mathrm{\Upsilon }}_G$E90

from Eqs. (87) and (89), the following relationships are deduced:

$k(p)={\nu }_M{k}_M(p_M)+{\nu }_F{k}_F(p_F)+{\nu }_G{k}_G(p_G)$E91

$\text{Φ}(p)={\nu }_M{\text{Φ}}_M(p_M)+{\nu }_F{\text{Φ}}_F(p_M)+{\nu }_G{\text{Φ}}_G(p_G)$E92

where Φ represents the potential of Kirchoff which is generically defined as

$\text{Φ}(p)={\displaystyle {\int }_{-\infty }^p}k(u)du$E93

If there is no fluid gain or loss in the porous medium, then ϒ = 0 and in consequence:

${\nu }_M{\mathrm{\Upsilon }}_M={\mathrm{\Upsilon }}_{MF}+{\mathrm{\Upsilon }}_{MG}$E94

${\nu }_F{\mathrm{\Upsilon }}_F=-{\mathrm{\Upsilon }}_{MF}+{\mathrm{\Upsilon }}_{FG}$E95

${\upsilon }_G{\mathrm{\Upsilon }}_G=-{\mathrm{\Upsilon }}_{MG}-{\mathrm{\Upsilon }}_{FG}$E96

where ϒ_MF is the input of fluid that receives the matrix from the fractured medium, ϒ_MG is the fluid input that receives the matrix of the vuggy media, and ϒ_FG is the contribution of fluid that receives the fractured medium from the vuggy media.

The system of differential equations is defined as follows:

$\frac{\partial (\rho {\theta }_M)}{\partial t}=\nabla \cdot p[\frac{\rho }{\mu }{k}_M(p_M)(\nabla p_M-\rho g\nabla D)]+\frac{\rho }{{\nu }_M}({\mathrm{\Upsilon }}_{MF}+{\mathrm{\Upsilon }}_{MG})$E97

$\frac{\partial (\rho {\theta }_F)}{\partial t}=\nabla \cdot p[\frac{\rho }{\mu }{k}_F(p_F)(\nabla p_F-\rho g\nabla D)]-\frac{\rho }{{\nu }_F}({\mathrm{\Upsilon }}_{MF}+{\mathrm{\Upsilon }}_{FG})$E98

$\frac{\partial (\rho {\theta }_G)}{\partial t}=\nabla \cdot p[\frac{\rho }{\mu }{k}_G(p_G)(\nabla p_G-\rho g\nabla D)]-\frac{\rho }{{\nu }_G}({\mathrm{\Upsilon }}_{MG}+{\mathrm{\Upsilon }}_{FG}).$E99

The contributions of fluid in each porous medium are modeled with the following relations:

${\mathrm{\Upsilon }}_{MF}=a_{MF}(p_F-p_M),$E100

${\mathrm{\Upsilon }}_{MG}=a_{MG}(p_G-p_M),$E101

${\mathrm{\Upsilon }}_{FG}=a_{FG}(p_G-p_F),$E102

where a_MF, a_MG, and a_FG are transfer coefficients at each interface, which may depend on the pressures on the adjacent media.

6.2. Monophasic flow saturated in triadic media

In the case of the monophasic flow saturated in triadic means, the continuity equations in each porous medium can be written as follows:

$\frac{\partial (\rho {\varphi }_i)}{\partial t}+\nabla \cdot p(\rho q_i)=\rho {\mathrm{\Upsilon }}_i,i=M,F,G\mathrm{.}$E103

Darcy’s law for each porous media takes the form

$q_i=-\frac{1}{\mu }{k}_i(\nabla p_i-\rho g\nabla D),i=M,F,G\mathrm{.}$E104

The substitution of Darcy’s law in the continuity equation leads to the following equations:

$\frac{\partial (\rho {\varphi }_M)}{\partial t}=\nabla \cdot p[\frac{\rho }{\mu }{k}_M(\nabla p_M-\rho g\nabla D)]+\frac{\rho }{{\nu }_M}({\mathrm{\Upsilon }}_{MF}+{\mathrm{\Upsilon }}_{MG}),$E105

$\frac{\partial (\rho {\varphi }_F)}{\partial t}=\nabla \cdot p[\frac{\rho }{\mu }{k}_F(\nabla p_F-\rho g\nabla D)]-\frac{\rho }{{\nu }_F}({\mathrm{\Upsilon }}_{MF}-{\mathrm{\Upsilon }}_{FG}),$E106

$\frac{\partial (\rho {\varphi }_G)}{\partial t}=\nabla \cdot p[\frac{\rho }{\mu }{k}_G(\nabla p_G-\rho g\nabla D)]-\frac{\rho }{{\nu }_G}({\mathrm{\Upsilon }}_{MF}+{\mathrm{\Upsilon }}_{FG}).$E107

When the fluid is considered at constant density and viscosity and the means of constant permeability, with D = z, we have

${\varphi }_M{c}_M\frac{\partial p_M}{\partial t}=\frac{k_M}{\mu }\text{Δ}{p}_M+\frac{1}{{\nu }_M}({\mathrm{\Upsilon }}_{MF}+{\mathrm{\Upsilon }}_{MG}),c_M=\frac{1}{{\varphi }_M}\frac{\partial {\varphi }_M}{\partial p_M},$E108

${\varphi }_F{c}_F\frac{\partial p_F}{\partial t}=\frac{k_F}{\mu }\text{Δ}{p}_F-\frac{1}{{\nu }_F}({\mathrm{\Upsilon }}_{MF}-{\mathrm{\Upsilon }}_{FG}),c_F=\frac{1}{{\varphi }_F}\frac{\partial {\varphi }_F}{\partial p_F},$E109

${\varphi }_G{c}_G\frac{\partial p_G}{\partial t}=\frac{k_G}{\mu }\text{Δ}{p}_G-\frac{1}{{\nu }_G}({\mathrm{\Upsilon }}_{MF}-{\mathrm{\Upsilon }}_{FG}),c_G=\frac{1}{{\varphi }_G}\frac{\partial {\varphi }_G}{\partial p_G}$E110

6.3. Triple porosity and triple permeability model

The porosity of each medium has been defined as the volume of the space occupied by the medium. However, the porosity can be defined as the volume of empty space in each medium with respect to the volume of the total space occupied by the porous medium as a whole. These new porosities will be denoted with subscripts in lowercase letters and clearly have

$\begin{array}{l}{\varphi }_m={\nu }_M{\varphi }_M,{\varphi }_f={\nu }_F{\varphi }_F,{\varphi }_g={\nu }_G{\varphi }_G\end{array}$E111

$\varphi ={\varphi }_m+{\varphi }_f+{\varphi }_g$E112

In an analogous way, the corresponding Darcy´s flow can be defined in each medium:

$q_m={\nu }_M{q}_M,q_f={\nu }_F{q}_F,q_g={\nu }_G{q}_G$E113

$q=q_m+q_f+q_g$E114

Eq. (113) implies that the permeability of the Darcy’s law in each medium is defined as

$k_m={\nu }_M{k}_M,k_f={\nu }_F{k}_F,k_g={\nu }_G{k}_G$E115

The nest system by Eqs. (108)–(110), by congruently changing the subscripts in uppercase by lowercase in the pressures, in terms of compressibility, is written as follows:

${\varphi }_m{c}_m\frac{\partial p_m}{\partial t}=\frac{k_m}{\mu }\text{Δ}{p}_m+({\mathrm{\Upsilon }}_{mf}+{\mathrm{\Upsilon }}_{mg}),c_m=\frac{1}{{\varphi }_m}\frac{\partial {\varphi }_m}{\partial p_m},$E116

${\varphi }_f{c}_f\frac{\partial p_f}{\partial t}=\frac{k_f}{\mu }\text{Δ}{p}_f-({\mathrm{\Upsilon }}_{mf}-{\mathrm{\Upsilon }}_{fg}),c_f=\frac{1}{{\varphi }_f}\frac{\partial {\varphi }_f}{\partial p_f},$E117

${\varphi }_g{c}_g\frac{\partial p_g}{\partial t}=\frac{k_g}{\mu }\text{Δ}{p}_g-({\mathrm{\Upsilon }}_{mg}+{\mathrm{\Upsilon }}_{mg}),c_g=\frac{1}{{\varphi }_g}\frac{\partial {\varphi }_g}{\partial p_g}$E118

with $p_m\equiv p_M,p_f\equiv p_M,p_g\equiv p_G,{\mathrm{\Upsilon }}_{mf}\equiv {\mathrm{\Upsilon }}_{MF},{\mathrm{\Upsilon }}_{mg}\equiv {\mathrm{\Upsilon }}_{MG},{\mathrm{\Upsilon }}_{fg}\equiv {\mathrm{\Upsilon }}_{FG},c_m\equiv c_M,c_f\equiv c_F,c_g\equiv c_G$.

The substitution of Eqs. (116)–(118) in Eqs. (100)–(102) leads to the system of differential equations that finalize the pressure in the matrix, fractured media, and vuggy media:

${\varphi }_m{c}_m\frac{\partial p_m}{\partial t}=\frac{k_m}{\mu }\text{Δ}{p}_m+a_{mf}(p_f-p_m)+a_{mg}(p_g-p_m),$E119

${\varphi }_f{c}_f\frac{\partial p_f}{\partial t}=\frac{k_f}{\mu }\text{Δ}{p}_f-a_{mf}(p_f-p_m)+a_{fg}(p_g-p_f),$E120

${\varphi }_g{c}_g\frac{\partial p_g}{\partial t}=\frac{k_g}{\mu }\text{Δ}{p}_g-a_{mg}(p_g-p_m)-a_{fg}(p_g-p_f),$E121

in which this system constitutes a triple porosity and triple permeability model. In polar coordinates, the system reduces to

${\varphi }_m{c}_m\frac{\partial p_m}{\partial t}=\frac{k_m}{\mu }\frac{1}{r}\frac{\partial }{\partial r}(r\frac{\partial p_m}{\partial r})+a_{mf}(p_f-p_m)+a_{mg}(p_g-p_m),$E122

${\varphi }_f{c}_f\frac{\partial p_f}{\partial t}=\frac{k_f}{\mu }\frac{1}{r}\frac{\partial }{\partial r}(r\frac{\partial p_f}{\partial r})-a_{mf}(p_f-p_m)+a_{fg}(p_g-p_f),$E123

${\varphi }_g{c}_g\frac{\partial p_g}{\partial t}=\frac{k_g}{\mu }\frac{1}{r}\frac{\partial }{\partial r}(r\frac{\partial p_g}{\partial r})-a_{mg}(p_g-p_m)-a_{fg}(p_g-p_f)$E124

6.4. Dimensionless variables

Now we will give a process of dimensionlessness to better manage the variables. This is a technique commonly used to make the parameters or variables in an equation having no units, bring to a range the possible values of a variable or constant in order that its value is known, and in this way, more manipulable.

The system of Eqs. (122)–(124) takes the following form after making the changes mentioned in the previous paragraph:

$(1-{\omega }_f-{\omega }_v)\frac{\partial p_{Dm}}{\partial t_D}=(1-{\kappa }_f-{\kappa }_v)\frac{1}{r_D}\frac{\partial }{\partial r_D}(r_D\frac{\partial p_{Dm}}{\partial r_D})+{\lambda }_{mf}(p_{Df}-p_{Dm})+{\lambda }_{mv}(p_{Dv}-p_{Dm})$E125

${\omega }_f\frac{\partial p_{Df}}{\partial t_D}={\kappa }_f\frac{1}{r_D}\frac{\partial }{\partial r_D}(r_D\frac{\partial p_{Df}}{\partial r_D})-{\lambda }_{mf}(p_{Df}-p_{Dm})+{\lambda }_{fv}(p_{Dv}-p_{Df})$E126

${\omega }_v\frac{\partial p_{Dv}}{\partial t_D}={\kappa }_v\frac{1}{r_D}\frac{\partial }{\partial r_D}(r_D\frac{\partial p_{Dv}}{\partial r_D})-{\lambda }_{mv}(p_{Dv}-p_{Dm})-{\lambda }_{fv}(p_{Dv}-p_{Df})$E127

where

${\omega }_f=\frac{{\varphi }_f{c}_f}{{\varphi }_m{c}_m+{\varphi }_f{c}_f+{\varphi }_v{c}_v},{\omega }_g=\frac{{\varphi }_v{c}_v}{{\varphi }_m{c}_m+{\varphi }_f{c}_f+{\varphi }_v{c}_v},r_D=\frac{r}{r_w}$E128

${\kappa }_f=\frac{k_f}{k_m+k_f+k_v},{\kappa }_g=\frac{k_v}{k_m+k_f+k_v}$E129

${\lambda }_{mf}=\frac{a_{mf}\mu r_w^2}{k_m+k_f+k_v},{\lambda }_{mv}=\frac{a_{mv}\mu r_w^2}{k_m+k_f+k_v},{\lambda }_{fv}=\frac{a_{fv}\mu r_w^2}{k_m+k_f+k_v}$E130

$p_{Dj}=\frac{2\pi h(k_m+k_f+k_v)(p_i-p_j)}{Q_0{B}_0\mu },t_D=\frac{t(k_m+k_f+k_v)}{\mu r_w^2({\varphi }_m{c}_m+{\varphi }_f{c}_f+{\varphi }_v{c}_v)}$E131

Eqs. (128)–(131) represent dimensionless variables so they have no units. The boundary conditions to which the previous model is subjected are

${\displaystyle \underset{r_D\to 1}{{\displaystyle \lim }}}{r}_D(1-{\kappa }_f-{\kappa }_v)\frac{\partial p_{Dm}}{\partial r_D}+r_D{\kappa }_v\frac{\partial p_{Dv}}{\partial r_D}=-1$E132

$p_w(t)=p_{Dm}(r_D,t){|}_{r_D=1}=p_{Df}(r_D,t){|}_{r_D=1}=p_{Dv}(r_D,t){|}_{r_D=1}=-1$E133

Substituting derivatives $\frac{\partial p_{Di}}{\partial t_D}$ by Caputo fractional derivatives$\frac{{\partial }^{{\alpha }_i}{p}_{Dj}}{t_D^{{\alpha }_i}}$ with $0<{\alpha }_i\leqq 1$, and $\frac{\partial }{\partial r_D}(r_D\frac{\partial p_{Dj}}{\partial r_D})$ by Riemann-Liouville complementary derivative, i. e., with infinite limit of integration, (also called Weyl derivative) $\frac{{\partial }^{{\gamma }_i}}{\partial r_D^{{\gamma }_i}}(r_D\frac{{\partial }^{{\beta }_i}{p}_{Dj}}{\partial r_D^{{\beta }_i}})$, with $1<{\gamma }_i+{\beta }_i\le 2$, $i=\mathrm{1,2,3}$; $j=v,f,m$; ${\alpha }_i,{\beta }_i,{\gamma }_i$ rational numbers.

The choice of the derivatives, Caputo and Riemann-Liouville (Weyl), obeys the nature of the problem and the ease with which they can be manipulated.

The monophase flow model with triple porosity and triple permeability is expressed as follows: For the matrix

$(1-{\omega }_f-{\omega }_v)\frac{{\partial }^{{\alpha }_1}{p}_m}{\partial t^{{\alpha }_1}}=(1-{\kappa }_f-{\kappa }_v)\frac{1}{r}\frac{{\partial }^{{\gamma }_1}}{\partial r^{{\gamma }_1}}\left(r\frac{{\partial }^{{\beta }_1}{p}_m}{\partial r^{{\beta }_1}}\right)+{\lambda }_{mf}(p_f-p_m)+{\lambda }_{mv}(p_v-p_m),$E134

for fracture media

${\omega }_f\frac{{\partial }^{{\alpha }_2}{p}_f}{\partial t^{{\alpha }_2}}={\kappa }_f\frac{1}{r}\frac{{\partial }^{{\gamma }_2}}{\partial r^{{\gamma }_2}}\left(r\frac{{\partial }^{{\beta }_2}{p}_f}{\partial r^{{\beta }_2}}\right)-{\lambda }_{mf}(p_f-p_m)+{\lambda }_{fv}(p_v-p_f),$E135

for vuggs

${\omega }_v\frac{{\partial }^{{\alpha }_3}{p}_v}{\partial t^{{\alpha }_3}}={\kappa }_v\frac{1}{r}\frac{{\partial }^{{\gamma }_3}}{\partial r^{{\gamma }_3}}\left(r\frac{{\partial }^{{\beta }_3}{p}_v}{\partial r^{{\beta }_3}}\right)-{\lambda }_{mv}(p_v-p_m)-{\lambda }_{fv}(p_v-p_f).$E136

We reduce this system by applying semigroup properties with respect to the order of the Weyl derivative, assuming: $0<{\alpha }_i\le 1$ and $1<{\alpha }_i+{\beta }_i\le 2$. Let $\omega =1-{\omega }_f-{\omega }_v;\kappa =1-{\kappa }_f-{\kappa }_v$, putting p_m = p; p_f = f; p_v = u; r = x; ${\eta }_i={\gamma }_i+{\beta }_i$. Then, the previous system can be expressed as

$\frac{{\partial }^{{\alpha }_1}p}{\partial t^{{\alpha }_1}}=\frac{\kappa }{\omega }\left(-{\gamma }_1\frac{1}{x}\frac{{\partial }^{{\eta }_1-1}p}{\partial x^{{\eta }_1}-1}+\frac{{\partial }^{{\eta }_1}p}{\partial x^{{\eta }_1}}\right)+\frac{{\lambda }_{mf}}{\omega }(f-p)+\frac{{\lambda }_{mv}}{\omega }(u-p),$E137

$\frac{{\partial }^{{\alpha }_2}f}{\partial t^{{\alpha }_2}}=\frac{\kappa }{{\omega }_f}\left(-{\gamma }_2\frac{1}{x}\frac{{\partial }^{{\eta }_2-1}f}{\partial x^{{\eta }_2}-1}+\frac{{\partial }^{{\eta }_2}f}{\partial x^{{\eta }_2}}\right)-\frac{{\lambda }_{mf}}{{\omega }_f}(f-p)+\frac{{\lambda }_{fv}}{{\omega }_f}(u-f),$E138

$\frac{{\partial }^{{\alpha }_3}u}{\partial t^{{\alpha }_3}}=\frac{\kappa }{{\omega }_v}\left(-{\gamma }_3\frac{1}{x}\frac{{\partial }^{{\eta }_3-1}u}{\partial x^{{\eta }_3}-1}+\frac{{\partial }^{{\eta }_3}u}{\partial x^{{\eta }_3}}\right)-\frac{{\lambda }_{mv}}{{\omega }_v}(u-p)+\frac{{\lambda }_{fv}}{{\omega }_v}(u-f).$E139

The above approach can be solved by numerical methods as finite differences along with a predictor-corrector, such as Daftardar-Gejji works, for example in [19] and compared with previous ones, such as that presented by Camacho et al. [18, 28], the approximations are significantly improved. However, there is still work to be completed; the optimal solution method has not been found and the best way to determine the appropriate order, so far numerical methods, has been used to estimate the order that best approximates measurements.

The application of the fractional calculation can be very useful for the modeling of anomalous diffusion phenomena in which the fractal structure better reflects the real conditions of the medium, as it is the case of the reservoirs in which because of its very nature it is difficult to find a structure Euclidian.