APPLICATION OF MONTE CARLO METHODS FOR SOLVING THE REGULAR AND DEGENERATE PROBLEM OF TWO-PHASE FILTRATION

Tastanov Meyrambek G.; Utemissova Anar A.; Mayer Fedor F.

2022

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

Математика и механика Tomsk State University Journal of Mathematics and Mechanics

№ 80

Original article UDK 519.245

doi: 10.17223/19988621/80/13

Application of Monte Carlo methods for solving the regular and degenerate problem of two-phase filtration

Meyrambek G. Tastanov1, Anar A. Utemissova2, Fedor F. Mayer3

12•3 Kostanay Regional University named after A. Baitursynov, Kostanay, Kazakhstan

1 tastao@mail.ru 2 anar_utemisova@mail. ru 3 maiyer@mail. ru

Abstract. One of popular mathematical models of filtration is the classical elastic regime model describing the nonstationary equilibrium filtration. It is also called the Muskat-Leverett model. Solving filtration problems by Monte Carlo methods makes it possible to find the solution of the problem at an individual point of the domain and to estimate derivatives of the solution. This paper is devoted to applying algorithms of the Monte Carlo method to problems of filtration. The Monte Carlo algorithms of random walk by spheres and on boundaries are used for solving the stationary problem of filtration of two immiscible inhomogeneous incompressible fluids in a porous medium and for estimating the solution and the derivatives of the solution of this problem.

Keywords: Monte Carlo method, continuity equation, Dirichlet problem, Markov chains, estimate of the solution and its derivatives

For citation: Tastanov, M.G., Utemissova, A.A., Mayer, F.F. (2022) Application of Monte Carlo methods for solving the regular and degenerate problem of two-phase filtration. Vestnik Tomskogo gosudarstvennogo universiteta. Matematika i mekhanika -Tomsk State University Journal of Mathematics and Mechanics. 80. pp. 147-156. doi: 10.17223/19988621/80/13

Научная статья

Применение методов Монте-Карло для решения регулярной и вырожденной задачи двухфазной фильтрации

Мейрамбек Габдулиевич Тастанов1, Анар Алтаевна Утемисова2, Федор Федорович Майер3

12•3 Костанайский региональный университет им. А. Байтурсынова, Костанай, Казахстан

1 tastao@mail.ru 2 anar_utemisova@mail. ru 3 maiyer@mail. ru

Аннотация. Одной из популярных математических моделей фильтрации является классическая модель упругого режима, описывающая фильтрацию с нестационар© M.G. Tastanov, A.A. Utemissova, F.F. Mayer, 2022

ным равновесием. Ее также называют моделью Маскета-Леверетта. Решение задач фильтрации методами Монте-Карло позволяет найти решение задачи в отдельной точке области и оценить производные решения. Данная статья посвящена применению алгоритмов Монте-Карло к задачам фильтрации. Алгоритмы случайного блуждания по сферам и по границам методом Монте-Карло используются для решения стационарной задачи фильтрации двух несмешивающихся неоднородных несжимаемых жидкостей в пористой среде и для оценки решения и производных от решения этой задачи.

Ключевые слова: метод Монте-Карло, уравнение неразрывности, задача Дирихле, цепи Маркова, оценка решения и его производных

Для цитирования: Тастанов М.Г., Утемисова А.А., Майер Ф.Ф. Применение методов Монте-Карло для решения регулярной и вырожденной задачи двухфазной фильтрации // Вестник Томского государственного университета. Математика и механика. 2022. № 80. С. 147-156. doi: 10.17223/19988621/80/13

Formulation of the problem

Let us consider the initially boundary problem for saturation and reduced pressure (5,p) in a given finite domain QeRn (n > 2) with the boundary SQ, Q = Q x [0, T],

G=dQ.x [0, 7]:

—■

m—= divCK0ciVs + Kyp + f0), (x,t )eO. (1)

dt

div(KVp + /) = 0, (x, t) e Q, (2)

s(x, t) = s0 (x, t), (x, t) e G, (3)

p(x, t) = p0 (x, t), (x, t) e G, (4)

s(i,0) = S°(i,0),i£Q, (5)

where the coefficients K0, a, K. f0 , K, and f , as well as the boundary and initial conditions, are given [1].

For the approximate solution of problem (1)-(5), two methods were proposed in [2]: Method 1.

3s — —

V,+1 - -m-+ div(K(x, s, )VsM ) + B(x, s, )VsM +

+D(x, s, )Vpi+1Vsi+1 = 0, (x, t) e Q,

s+ (x, t) = s0 (x, t), (x, t) e G, (7)

s.+1(x,0) = s°(x,0),xeQ, (8)

L2pM = div(^(x, 5,. yvPl+1 + f{x, 5,.)) = 0, (x, О e Q, (9)

pM(x,t) = p0(x,t\(x,t)eG (10)

where mK = Ко a,mB = К V m + f'0-bf , mD = kb', b(s) = K\K~l, k = fa\ + ka, fa are the phase permeabilities for a homogeneous isotropic soil (i = 1, 2). Method 2.

After dividing the time interval [0, T] into N parts (т = T/N) for each time layer te [It, (l + 1)т], l = 0, ..., N - 1, the initially boundary value problem is solved for the functions s/+1 (x, t), p/+1 (x), (s (x, 't) = s' (x), s0 (x, 0) = s0 (x) = s0 (x, 0)) :

ds_h

V/+1 =m—^- + äiv(K(x,s!)VsM) + B(x,s!)VsM +

+D(x, s' )VpMVsM = 0, (x, t) e Q,

si+i(x, t) = s0(x, t), x e G, (12)

sM(x,h) = s!, ieQ, (13)

LaPm = div(K{x,s' ))VpM + f{x,s')) = 0, (x,i) e ft (14)

Pi+i(x) = p0(x, 'x), x e5Q, (15)

Let us describe the general scheme of using Monte Carlo algorithms for methods 1 and 2 [3].

For method 1, the Dirichlet problem is first solved for the linear elliptic equation (9),

(10) for p+\(x, t) at a given saturation value Si and fixed t = to, in particular, to = 0. Then, equation (6) is split only with respect to the time variable, i.e., for the iteration index i + 1 the time interval [to, T] is divided into M parts (t = (T-to)/M) and for each time layer tj = Tj + to, j = 0, ..., M - 1, by use of the implicit difference scheme for (6), the Dirichlet problem for an elliptical equation for the variable sj+1 (x) is obtained. The corresponding boundary and initial conditions (7) and (8) are written in the form

sj^1 (x) = sj+1 (x), x e 9Q, j = 0,M -1,

s,°+1 (x, 0) = s0 (x, 0), x efi. Now, omitting the subscript, we obtain

vi+1 - sJ — -

-m-+ divCK(x,si)W+1) + fi(x,si)W>+1 +

x (16)

+D(x,sj)Vpj+1 -Vsj+1 = 0, j = 0,1,...,M-1,

sj 1 (x) = sj1 (x), x e 9Q, j = 0,1, ..M -1, (17)

s0 (x, 0) = s0(x), x e Q. (18)

For method 2, similarly to method 1, the Dirichlet problem for the linear elliptical equation (14), (15) is first solved for pm(x) at a given saturation s'; in particular, at l = o from (13) we obtain s:(x, o) = so(x), x eQ. Now, using the purely implicit scheme, approximating only with respect to the time variable, for the initially boundary problem

(11)-(13) for the variable sm(x) = s'+1(x), (l = j), we obtain the Dirichlet problem for the elliptical equation, i.e., problem (16)-(18).

Thus, to determine p and s, one has to solve the Dirichlet problem for the elliptical type [4].

Omitting indices of the time layer j, we obtain for the determination of pressure p(x) (p (x) =pi+1(x)) the problem

diviAXv. s))Vp + f (x, s)) = 0, x e Q, p(x) = p0 (x), x e 5Q

or

K(x, s))Ap(x) + Y C,. (x, s) ■ ^^ + g(x, s) = 0, x e Q, (19)

/=1 CtXi

p(x) = p0 (x), x e 5Q, (20)

8

where C,(x, s) = —K(x. s). g(x. s) = div / (x. s): when determining the saturation

dxt

s(x) (s(x) = s^OO), we obtain

div(^(x, s))Vs)x)) + (B(x, s) + D(x, s)Vp(x))Vs(x) -ms(x) = -ms(x), x e i! s(x) = i0(x), x e SQ , s(x) = s°(x), x e Q

or

" ds(x) " 8s (x)

K(x, s) A?(x) + X E, (x, s) —^ + X B, (x, s) +

5x ,=, 8x,

^dpix)ds(x) + D(x, s) • У ———m(x)s(x) = —m(x)s(x), x e Q ,=1 6x 8xi

(21)

i(x) = 50(x), ig^Q (22)

7tl( X)

where .v(x) is a known function by virtue of the initial data, m(x) = ——-,

t

m(x) 0 _

m(x) = , Ej(a) = — K(x, 5), Bj(x, 5) are components of the vector B(x,i).

t 0xj

Let us construct a random process and algorithm for solving problem (19), (20). Consider the Dirichlet problem for a function 9 continuous at the boundary 5Q, a measurable functionf and an elliptical operator L

Lu(x) = —f (x), x eQ, (23)

u(x) = ^(x), x e 0Q. (24)

Let us construct random processes for numerical finding of the solution u. We suppose that the domain Q and operator L are such that problem (23), (24) has a unique solution, continuous in Q and regular in Q, for any sufficiently smoothf and 9 [5]. It is known that the integral representation

u(x) = J k(x, y)u(y)dy + J A(y, x)f (y)dy, (25)

V(x) V ( x)

where k(x, y) = N A(y, x) > 0, A(y, x) - is a Levy function, Ny is an operator adjoint to the operator L, V(x) is an ellipsoid,

V(x) = VR (x) = Iy : ct(y, x) = (A-1 (x)(y — x), (y — x))2 < R(x) |

is valid for the solution u(x) of the boundary problem (23), (24).

Here, R(x) is a maximum radius ball with a center at the point x lying in Q and A is the matrix of higher coefficients of the operator L; the matrix is symmetric [6].

Representation (25) is called the mean value theorem. Note that if the coefficient C < 0 at u(x) in equation (23) of the operator L, then the kernel k(x, y) is substochastic,

i.e., J k(x,y)dy < 1 [7]. Representation (22) allows one to construct unbiased esti-

V ( x )

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

mates for the solution of problem (23), (24). Any regular solution of problem (23), (24) satisfies equation (25) and boundary condition (24). In connection with this, we define operator K acting on functions from c(q) by the formula

[ k(x,y)u(y)dy, x eQ, (Ku)(x) = j v(x) (26)

u(x), x e SQ.

Consider the following problem: for ^e C(SQ) and F e C(Q), find u e C (q) such that

|u(x) =(Ku)(x) + F(x), x e Q,

[ u(x) = x), x e SQ.

If F(x) = J A(y, x)f (y)dy, then the solution of problem (27) is a solution for (23), (24).

v (x)

In [8], the uniqueness theorem was proved for the class of kernels k (x, y) = Ny A( y, x) of the operator K acting by formula (26) and having the properties

1) mes(V(x) A V(x2)) ^ 0 as x1 ^ x2,

2) diam(V(x)) ^ 0 for x e SQ,

3) kernelk(x,y) is substochastic and weakly polar, k(x,y) = W(x,y)/\x-y|n1, n > 3, where the function W(x,y) can be continuously extended from QxV(x) to QxQ,

4) |K||r = 1; then the integral operator acting in Lm (Q) by the formula (Ku)(x) = J k(x,y)u(y)dy and satisfying properties 1), 2), and 3) maps functions

V ( x )

bounded in Q into continuous ones. It is easy to establish that for an operator satisfying properties 1), 2), 3) the following conditions hold:

a) for any u e C(Q),

J k(x, y)u(y)dy ^ u(x0) as x ^ x0 e SQ,

V (x)

b) J k(x, y)dy ^ 1 as x ^ x0 e SQ.

V (x)

From the above, it follows that the operator defined by formula (26) is bounded in C(Q), i.e., problem (27) can be solved in C(Q) and problem (27) has not more than one solution in C(Q).

Construction of unbiased estimates of solution (27)

Let V(x) and kernel k(x, y) satisfy properties 1), 2), 3) and additional requirement b). By virtue of the last property, the spectral radius of the operator K is equal to unity; therefore, one cannot use the standard estimates for solving integral equations of the second kind by Monte Carlo methods. Here, the estimation scheme is based on the martingale theory. In this case, it is easy to analyze the variance of the estimates.

A terminating Markovian chain is determined with a transition density P(x, y) = k(x, y), y e V(x). The probability of termination q(x) = 1- J k(x, y)dy

V (x)

tends to zero at x e 5Q; therefore, the trajectory of the chain can have an infinite

length. For the set of trajectories B having an infinite length, the following Lemma is valid:

Lemma. Let problem (27) be solvable for F = 0,9 = 1. If there exists a function 0 < H e C(Q) such that problem (27) for F = H, 0 is solvable, then for a chain

{xm with transition density P(x, y) almost all trajectories of infinite length approach the boundary oil:

P {dist(xn, 0Q) ^ > 0 |i?} = 1, and P [b) > 0. (28)

Proof. From Lemma 2.3.3 (see Lemma 2.3.3, [4] it follows that for the solution u(x) of problem (27) with F = H, 9 = 1 the maximum principle is valid: u(x) reaches

the smallest value on the boundary of the domain. Hence, u(x) > 0. Let be the indicator of an event (the moment of chain termination > i), {An be a sequence of

c-algebras generated by the chain up to the time moment n. Consider a sequence of random variables

n-1

• = Z H (x ^ +Xnu(xn ).

i=0

The sequence {n }J_o is a positive martingale with respect to {A }™_0. Indeed,

n-1

M |A*-i} = ZH(Xi)Xi +Mxn-t (Xnu(xn)} =

i=0

n-1

= Z H(Xi )Xi +Xn-1 • J k(Xn-1 > Xn )u (Xn )dxn = •-1.

i=0 T ( xn-1)

Then, by the martingale convergence theorem [9], there exists a random value such that Mxr\m<+<x and r|); —» r|„ as w-> 00 with probability 1. Consequently, almost everywhere on the set B II (x„ j —» 0 as n —> 00. Let

dQs ={xeQ':dist(x,aQ)<6},Qs =Q\dQs. It is evident that H (x) > const = c(e) > 0 011 CI. . If /], is a subset of trajectories from B such that H(xn) ->0 as n 00, thenPB (¿KfiJ = 0. Hence, if P^fiJ > 0, then

Px ldist(x„,aQ)->o|fi) = Px ldist(x„,aQ)->o|fii I.

\ ri—>cO / \ ri—>cO /

not

Let X = (x0,x.,...,x ,) e B., but dist(x ,3Q)->0. Then there also exists £0 and

«—»co

an increasing sequence {nk }™=0 such that dist(x, 9Q) > e0. Hence, H(x ) > c(e) > 0, so XeBv It follows that {x: dist(xn, cO) >()} id B . Therefore, if P (Bl) > 0, then P {dist(xn, cQ) ^ > 0\bx } = 1.

Now let us prove that Px (Bj) > 0. By condition of the lemma, there exists a solution to problem (27) with F = H and 9 = 1. Let us denote it as v(x). It can be shown that inf v(x) > 0, x e Q'. A martingale r\u = 7_„v(.v„) is uniformly integrable; therefore, by the martingale convergence theorem, v(x) = MX limr^ = 0, ifPX(B) = /' (/},) = 0. And

o

this contradiction proves the inequality Px (B1) > 0. The lemma is proved.

If we take H(x) = J A(x, y)dy as H(x), then we can construct estimates for the so-

V ( x)

lution u(x) of problem (27) on trajectories of the chain {xm with a transition density

\k(x, y), y e V(x),

x, y) = 'm

[0, y tV (x).

The sequence of estimates is determined by the equality

m-1

■qm F(x )%,■ +%,■u(x), where is the event indicator {the moment of the chain

1=0

termination > /'}. Obviously, Mx^ = u(x), i.e. estimates ^ are unbiased. The sequence {'nm}^_0 forms a martingale with respect to {Am - a sequence of c-algebras. Am is generated by the chain up to the time instant m. The last statement is proved in the same way as the Lemma. From this we have

Corollary. For a Markovian chain {xn determined by the transition density P(x, y) = k (x, y) = Ny A( y, x), (28) is fulfilled.

Let be the moment of chain termination, t2 be the moment when the chain enters the 5-neighborhood of the boundary xs = min(^, x2). A sequence {^m}™=0 of unbiased estimates for the solution u(x) of problem (23), (24) is called admissible if there exists a sequence of c-algebras {Ym }m=0 such that Am cYm and Am cYmt1, and E,m has the form = + xmu(xm), where \m Ym are measurable. For an admissible se-

1 w

quence of estimates , we define a random variable by the equality

x,5), (29)

where x denoted a border point closest to x . The definition is correct since ts < +w by virtue of the above Lemma. We finally obtain

Theorem. If an admissible sequence of estimates }W_0 forms a square integra-

'm.»m=0

ble martingale with respect to a family of c-algebras, {Ym }W=0, then the random variable is a e(5)-biased estimate for u(x), its variance is a bounded function of the parameter 5 (e(5) is the modulus of continuity of the function).

Proof. Let x denote the indicator of the event {x = т2} . By the theorem about the free choice transformation [6], Mx^ = u(x), therefore, \u(x)-Mx¡5| = -

< MX\u(\) -ф(хТ5)| < e(5), i.e. e(5) and the bias of the estimate is proven. Now let us prove that the variance is bounded. We introduce an abbreviation 4 for £ , then we obtain

= Mx (¡5 -Mx¡5 )2 = Mx (¡Ц - £ + £ -u(x) + u(x) -Mx^5)2 < < 4 [Mx (¡5 - г)2 + Mx (г - u(x))2 + (u(x) - Mx^5 )2 ] < 8s2 (5) + D

By virtue of square integrability of the martingale {¡n }™=1 , sup Mxi,l < there-

n

fore, the D£s < +да. The theorem is proved.

Consider now the sequence of estimates {\}" .

Lemma. If problem (27) is solvable when F(x) is replaced by |F(x)| and F 2( x), ф = 0, then the martingale {л„}Г=о is square integrable.

Proof. Let us put s = XF(x); and let v^ denote a solution of problem (27). Under conditions of the Lemma, F(x) ^ 0 as x ^ x0 e SQ; therefore, F(x) is bounded in Q. Let F = ||F|| = max|F(x)|; then, for the function H(x) = |F(x)|/ F, equation (27) is

solvable at ф = 0 разрешимо. Its solution vG0 can be represented as a series vG,0 = XKnG. Since 0 < G <1, the series XKnG2 also converges and yields a solu-

т=0 т=0

tion vG2 o of problem (27). Therefore, there exists v^,2 o. Then,

\ = [XS +Xnu(x„)| < 2 ^s,| + 2X„u2(x„), MxXnu\xn) < ||u||C(Q1),

f n-i n-1 n-2 n-1

Mx I Xs, I <MxXs,.2 + 2MxXs, X sk = I + 2I2,Ii < Vf2,0.

V ,=0 ) ,=0 ,=0 k=,+1

To estimate the second summand, we use the Markov property:

n-2 f f n-1 ] Л n-2

I2 <XMx W*[XKWf <MxS,V|f|,0(x,).

V

k=,+1

,=0

Therefore, I2 <

V|F| ,0

2 2 II ||2 II 112

. Thus, sup M \ < 2 u + v + V , . Lemma is

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

C(Q') n '0 11 1 '' 11

proved.

Now, it remains to determine for problem (23), (24) a sequence of unbiased estimates which are obtained from according to the estimate F(x) by a single random

node with the density A(y, x)/ h(x), h(x) = J A(y, x)dy. In this case, the sequence

V ( x)

m-1

of unbiased estimates h(x)f (y,)x, +xmu(xm) forms a martingale with respect

"Dm

to {Ym}™=0. Here, y0,y,...,ym_15... are independent random vectors with the densities A(y,x)/), Y is a o-algebra generated by x0,x,..., xm,y,y,..., ^-r The Lemma below follows from the previous one.

Lemma. The martingale {^m}™=0 is square integrable.

Finally, applying the theorem to the martingale 0, we obtain

Theorem. Let e(5) be the modulus of continuity of the solution u(x) of problem (23), (24); then, the estimate ^, determined according to {^m}™_0 by formula (29), is e(5)-biased for u(x). D^s is a bounded function of the parameter 5.

Conclusions

To implement the algorithm in practice, it is necessary to learn to model the chain {xm}™_0 on trajectories of which estimates of the solution are constructed, and algorithms of simulation of the Markovian chain A(y, x)/ h(x) are based on the von Neumann selection method. Modeling distributions requires special investigations, especially in cases where it is necessary to model the distributions regularly and repeatedly. Such investigations are usually carried out for each specific equation if it is solved several times.

References

1. Antontsev S.N., Kazhikhov A.V., Monakhov V.N. (1983) Kraevye zadachi mekhaniki neod-

norodnykh zhidkostey [Boundary value problems of mechanics of inhomogeneous fluids]. Novosibirsk: Nauka.

2. Tastanov M.G., Shakenov K.K. (1997) Ob odnoy statsionarnoy zadache dvukhfaznoy fil'tratsii

i primenenie algoritmov metodov Monte-Karlo [On one stationary problem of two-phase filtration and application of algorithms of Monte Carlo methods. Mathematical modeling of scientific, technological, and environmental problems in the oil and gas industry]. Almaty: Al-Farabi KazNU.

3. Ermakov S.M. (1975) Metod Monte-Karlo i smezhnye voprosy [Monte Carlo method and

APPLICATION OF MONTE CARLO METHODS FOR SOLVING THE REGULAR AND DEGENERATE PROBLEM OF TWO-PHASE FILTRATION Текст научной статьи по специальности «Математика»

Аннотация научной статьи по математике, автор научной работы — Tastanov Meyrambek G., Utemissova Anar A., Mayer Fedor F.

Похожие темы научных работ по математике , автор научной работы — Tastanov Meyrambek G., Utemissova Anar A., Mayer Fedor F.

Текст научной работы на тему «APPLICATION OF MONTE CARLO METHODS FOR SOLVING THE REGULAR AND DEGENERATE PROBLEM OF TWO-PHASE FILTRATION»