On Solving the Regular Problem of Two-Phase Filtration Taking Into
Account The Energy and one Model Problem of Filtration in Potentials
by Monte Carlo Methods
M. G. TASTANOV, A. A. UTEMISSOVA*, F. F. MAIYER, D. S. KENZHEBEKOVA,
N. M. TEMIRBEKOV
Kostanay Regional University named after A.Baitursynuly,
Kostanay,
REPUBLIC OF KAZAKHSTAN
*Corresponding Author
Abstract: - It should be noted that in some practical tasks, it is impossible not to take into account the
temperature change. In this case, the energy equation should be added to the filtration equations. The
algorithms of «random walk by spheres» and «random walk along boundaries» of Monte Carlo methods are
used to solve regular degenerate filtration problems of two immiscible inhomogeneous incompressible liquids
in a porous medium. The derivatives of the solution are evaluated using Monte Carlo methods. A model
problem of filtration of a two-phase incompressible liquid with capillary forces is considered.
Key-Words: - Monte Carlo method, Dirichlet problem, Markov chains, Levy function, unbiased estimation of
the solution, integral operator.
1 Introduction
In the previous work, two approximate methods for
solving two-phase filtration problems - regular and
degenerate were proposed. In each of these
methods, a linear differential problem is solved at
an intermediate stage, which can be approximated
by known difference schemes. In both methods, a
linear elliptic problem is solved first concerning a
given saturation value, [1]. If the temperature
change is taken into account, then the energy
equations are added to the filtration equations, [2],
[3].
2 Formulation of a Regular Filtration
Problem Taking Into Account the
Energy
It is known if the conditions are met:
1
),(,0
,),(,0
StxnF
QtxFdivx
(1)
then the filtering task is called regular. It follows
from (1) that almost everywhere
in
0 0 0 1
0 min ( , ) ( , ) max( ( , ) 1 1.s x t s x t s x t
and for a stationary
task
01
0 ( ) 1 1, .s x x
On
regular solutions
0 0 1
( , ) . 0.K const
Also, in
the regular case, the function
( , ) .,
s
a x s C const
because
( ) .
c
p x const
is
in a homogeneous environment. So, we got to
determine
s
again the Dirichlet problem for the
Poisson equation. We
denote
6 2 5
C C C
и
61
(1/ ) .C divW W
Then
)3(),(
)2(,
0
1
onxss
inWs
The problem (2) – (3) is solved using the
«random walk by spheres» algorithm. In [2], (see
paragraph 11. Unsolved problems, p. 5, p. 303) it is
noted that in practical problems it is impossible to
neglect temperature changes. That is, to the
filtration equations
)5(,0),()(
)4)(,()( 1010
psVdivfpKdiv
psVdivfpKsaKdiv
t
s
m
Received: April 21, 2023. Revised: November 17, 2023. Accepted: December 13, 2023. Published: March 8, 2024.
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2024.23.18
M. G. Tastanov, A. A. Utemissova,
F. F. Maiyer, D. S. Kenzhebekova, N. M. Temirbekov
E-ISSN: 2224-2880
154
Volume 23, 2024
the energy equation should be added
)6(,)(
Udiv
t
c
where
is the temperature of the mixture,
и c
heat capacity and thermal conductivity
of the mixture,
average density of the
mixture,
U
average speed.
We consider in a homogeneous isotropic
medium
with a boundary
stationary
filtration problem taking into account the
temperature
)14(
)13(,)(
)12(,)(
)11(,0)(
)10(,)()()(
)9(,)()()(
)8(,)()(
)7(,)(-hp(x)
0
0
2
3
0
1
on
inUpdiv
onxss
inWsacdiv
inxFxVbxW
inxfxpcxV
onxpxp
inx
We construct an algorithm for solving the
stationary problem of two-phase filtration taking
into account temperature in a homogeneous and
isotropic medium
with a boundary
.
.For
doing so, first, by solving the problem using the
«random walk by spheres» algorithm, [4]
)16()()(
)15(,)()(
0
1
onxpxp
inxhxp
we evaluate
( ), ( ) и ( ),p x p x p x
(here
and further index
shows on
bias of estimates),
[5].
We determine:
)()()(
),()()(
3
3
xfdivxpcxVdiv
xfxpcxV
e
e
and
)()()(
),()()(
xFdivxVdivxWdiv
xFxVbxW
in the area of
.
We assume that the filtering task is regular.Then
to determine the saturation
s
we get the task
)18(,)(
)17(,
0
1
onxss
inWs
Where
16
( ) (1/ ) .W x c divW
Solution of the
problem (17), (18)
()sx
we also evaluate it by
«random walk by spheres». Now we consider the
stationary energy equation in the domain
)19(,0)( inUdiv
We suppose that on the boundary
areas
0
and
, и U
constant
values.
In this case, they get the task
)21(,
)20(,0
0
1
7
on
in
x
cn
ii
Where
7., ( 1,2,3)
i
с const u i
vecto
r components
.U
Modeling of the Markov chain,
on the trajectories of which estimates of the
solution of the problem (20), (21) are constructed,
is based on the von Neumann selection method, [6].
We consider a more general case. Let be given
an elliptic operator
n
jj
j
n
ji i
ij xc
x
xb
xjx
xa
11,
2)22(),(
~
)()(
where are the
coefficients
( ), ( )
ij j
a x b x
and
()cx
real
measurable functions defined in
. Suppose that
the matrix of higher coefficients
,1
() n
ij ij
A x a
is
symmetrical. Through
()Dx
let's denote a
collection of regions where the Green's function for
a given operator is known.For example, a collection
of balls of maximum radius centered
in
. ( ) { : ( )}.x D x y r R x
Let
12
( ) ( ) ... ( )
n
x x x
matrix
eigenvalues
( ),Ax
r
distance between
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2024.23.18
M. G. Tastanov, A. A. Utemissova,
F. F. Maiyer, D. S. Kenzhebekova, N. M. Temirbekov
E-ISSN: 2224-2880
155
Volume 23, 2024
points
x
and
, .y r x y
Then by virtue of the
ellipticity condition
1( ) ,x
and from the
conditions on the coefficients
()
ij
ax
it follows
that
1
( ) .
nx const
We determine the function
( , )xy
equality
),(max
,)))()(((),(
1
1,
2
1
yxR
yxyxyayx
y
n
ji jjiiij
For an arbitrary summable on
1
[0, ]R
we
determine the functions
()p
0
( ) ( )
R
q R p d
and
1
1/2
( ) ( ) ( 2) ( ) ,
n
x q R n A x
where
n
the surface area of a sphere of radius 1 in
.
n
R
Now for the
case
2
,1
( ) ,
n
ij
ij ij
ax xx
where
Laplace
operator, we describe the process of modeling a
Markov chain.To do this, we use the following
density
(majorant)
1
1
2( , ) 2 ( )( 2) ( ) n
p x y x n p
.In
general, as the majorant density, by choosing the
density
2( , ),p x y
it is easy to construct a Markov
chain using the Neumann method, since the
modeling
2( , )p x y
does not cause difficulties. For
modeling the main Markov chain with
density
1( , )p x y
we use the
inequality
12
( , ) ( , ),p x y p x y
where
[0,1]
evenly distributed on[0,1]random
variable. Modeling efficiency
1( , )p x y
with
majorant
2( , )p x y
equal to 1/2.Hence, the average
number of samples to obtain one implementation is
2, [7].
3 Solutions of One Model Filtration
Problem in Potentials
Problem statement.We consider the model problem
of filtration of a two-phase incompressible liquid
taking into account capillary forces. In the
cylinder
(0, )T
with a boundary
(0, )QT
we look for a solution to the
problem
)23(,0)()( 21 QinffMgradPdivNgradRdiv
Qinff
NgradPdivMgradRdiv
t
R
,0
)()(
21
(24)
),()0,(2 xxR
(25)
onRP
n,0)(
(26)
onRP ,
(27)
where
1 2 1 2
0, .M k k N k k
Here the desired functions are
PandR
1 2 1 2
2 , 2 .R U U P U U
We discretize (23)-(27) only by the time
variable and we will assume that the moment of
time
tn
we know
n
P
and
n
R
.Then to
determine
1n
P
on time layers from(23), (27) we
will have a Dirichlet problem for an elliptic
equation
infNgradPdiv nn ,)( 1
(28)
onRP nn ,
1
(29)
where
12
( ) .
n n n n
f div N gradR f f
Using (24), (25) and (26) to determine
1n
R
on
time layers, we obtain the Neumann problem for
the elliptic equation, [8],
,)(R 111n ingRgradMdiv nn
(30)
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2024.23.18
M. G. Tastanov, A. A. Utemissova,
F. F. Maiyer, D. S. Kenzhebekova, N. M. Temirbekov
E-ISSN: 2224-2880
156
Volume 23, 2024
.,
1
11
ong
n
P
n
Rn
nn
(31)
Where
1 1 1 1
12
( ( ) )
n n n n n
g div N gradP f f R
.
We solve the problem (28), (29) using either
«random walk by spheres» or «random walk along
boundaries» algorithms, [9], [10], [11], [12], [13],
[14]. The derivatives of the solution are evaluated
using Monte Carlo
methods
1n
gradP
and
1
( ).
n
div N gradP
We will assume that the area
is bounded
convex, and the boundary
is smooth enough
that Green's formula is valid. Also we
assume
0.M c const
Then
n 1 1
( gradP ) n
div M M P
and in these
assumptions using the Levy function for the
operator
( ), ( ) 0a x a x
it is possible to
determine an integral equation for solving the
problem (30), (31) with the norm of the integral
operator acting in
( ),c
a smaller unit.
Equation (30) is written as:
1 1 1
( ) , .
n n n
R a x R g x
(32)
where
1
( ) 1/ 0, / .
n
a x const g g a
Let
2
11
( ) , 0.a x c c const
Then the
function
21
1
( , ) exp( )( ( 2) ) ,
m
m
v x y c r m r
where
r x y
is a Levy function for the
operator
( ),ax
in doing so, this operator is
formally self-adjoined. The integral equation for
the problem (32), (31) will have the form, [9]
xFydyRyxpyxqxR nn ,)()(),(
ˆ
),(1()( 11
where
11
1
2
cos
( ) , ,
1 ( )
ˆ( , )
( )/ , ,
xy
m
m
m
k r y
Ix
p x y r
k r r y
1
2
0,
( , ) ( ) exp( ) ,
( )( 2)
y
q x y a y r c r y
k r m
12
1 ( )
( ) ( ( ) ( )).
m
Ix
F x F x F x
Here
()Ix
is the boundary indicator,
measure
determined on
- algebra of Borel
subsets
equality
( ) ( ) ( ), A A S
- Lebesgue
measure in
m
R , S
-surface area,
),exp(
2
)3(
)(
),exp()
2
1()(
1
1
2
1
2
1
1
1
rc
m
mcc
rk
rc
m
rc
rk
.)(
)exp(
)(
,)(
)exp(
)(
1
21
2
1
21
1
dyyh
r
rc
xF
dyyg
r
rc
xF
n
m
n
m
Since the area
is convex, then
( , ) 0.p x y
If
11
( ) 0, ( ) 0, ( ) 0,
nn
a x h x g x
then
11
n
R
is a solution to the Neumann problem.
It is clear that
0 ( , ) 1,q x y
because
2
1
( ) .a x c
If
( ) 0,ax
then the norm of the integral
operator in equation (10) acting in
()c
less than
one. Consequently, the Neumann-Ulam scheme is
applicable to the integral equation. Now it is
possible to construct unbiased estimates of the
solution of the problem (10), (9), and with finite
variances, [15], [16], [17], [18]. To solve the
problem (8), (9), we can apply a scheme
generalizing the Neumann-Ulam scheme – a
method for isolating the main part of the operator.
All processes of the «random walk by spheres»
type with allocation
-the neighborhood of the
boundary fits into the scheme of allocation of the
main part of the operator with a small norm of the
integral operator.
In the work [14], an algorithm of random walk
along the boundary for the external Neumann
problem for a self-adjoined elliptic equation of the
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2024.23.18
M. G. Tastanov, A. A. Utemissova,
F. F. Maiyer, D. S. Kenzhebekova, N. M. Temirbekov
E-ISSN: 2224-2880
157
Volume 23, 2024
second order of the general form is proposed.
Having solved the tasks (6), (7) and (8), (9)
for
( 1)n
-th step in time, we move on to the next
time layer, etc. When numerically solving problems
of filtration of a two-phase incompressible liquid,
some features of these problems should be taken
into account. For example, it is necessary to take
into account when constructing difference schemes
a feature associated with highly varying
discontinuous coefficients in sub domains. This
difficulty associated with the choice of a step is
easily eliminated when Monte Carlo algorithms
associated with the modeling of Markov chains are
used to evaluate the solution. In this case, the
transition to those points where the coefficients
suffer a gap is not carried out. That is, the
trajectories of the chain should «be able» to start at
those points, and at the transition
yx
get to
those points in the area where the coefficients do
not have break points.
4 Conclusion
The Masket–Leverett model describing the process
of liquid filtration in a porous medium is a system
of equations taking into account saturation and
pressure obtained using nonlinear laws, solvable
Monte Carlo methods and the method of
differences in probability. We were able to apply
the algorithms of «random walk by spheres» and
«random walk along boundaries» Monte Carlo
methods to solve regular, degenerate, stationary
and non-stationary filtration problems of two
immiscible inhomogeneous incompressible liquids
in a porous medium, [19], as well as using Monte
Carlo algorithms, we solved the filtration problem
taking into account temperature (i.e., an energy
equation is added to the filtration equations) and
one model filtration problem in potentials.
References:
[1] Shakenov K.K. (2007) Solution of problem
for one model of relaxational filtration by
probability–difference and Monte Carlo
methods. Polish Academy of Sciences.
Committee of Mining. Archives of Mining
Sciences. Vol. 52, Issue 2, Krakow, 2007. P.
247 – 255.
[2] Antontsev S.N., Kazhikhov A.V.,
Monakhov V.N. (1983). "Boundary value
problems of mechanics of inhomogeneous
fluids" (Kraevye zadachi mekhaniki
neodnorodnykh zhidkostey). Novosibirsk:
Nauka.
[3] Loitsyansky L.G. (2003) Mechanics of
liquid and gas.7th edition, revised. – M.:
Bustard, 2003, p.840.
[4] H.Amann, G.P. Galdi, K. Pllecas and V.A.
Solonnikov (1997), Navier-Stokes
Equations and Related Nonlinear Problems.
Proceedings of the Sixth International
Conference NSEC-6, Palanga, Lithuania,
May 22–29, p.440.
[5] Todor Boyanov, Stefka Dimova, Krassimir
Georgiev, Geno Nikolov. (2007) Numerical
Methods and Applications, Springer Nature,
p.732.
[6] New Directions in Matematical Fluid
Mechanics. The Alexander V. Kazhikhov
Mtmorial Volume. (2010) Springer Science
and Business Media, p.414.
[7] 7.Ermakov S.M., Nekrutkin V.V., Sipin
A.A. (1984), Random processes for solving
classical equations of mathematical physics.
Nauka, –M., 1989, p.208.
[8] Shoujun Xu, Hao Wu. The Gevrey
normalization for quasi-periodic systems
under Siegel type small divisor conditions
(2016), Journal of Differential Equations,
Vol. 284, p.1223-1243.
[9] Alexander Keller, Stefan Heinrich, Harald
Niederreiter Editors. Monte Carlo and
Quasi-Monte Carlo Methods 2006 2008th
Edition, Kindle Edition (2008). Springer,
p.708.
[10] Mikhailov G.A. (1987). Optimization of
Monte Carlo weight methods. Nauka, – M.,
1987, p.239.
[11] Mikhailov G.A. (1974). Some questions of
the theory of Monte Carlo methods. Nauka,
Novosibirsk, 1974, p.145.
[12] Dynkin E.B. (1963) Markov processes.
Fizmatgiz, –M., 1963, p.432.
[13] Dynkin E.B., Yushkevich A.A. (1967)
Theorems and problems about Markov
processes. Nauka, –M., 1967, p.232.
[14] Sobol I.M. (1973) Numerical Monte Carlo
methods. Nauka, –M., 1973, p.312.
[15] Meyer P.A. (1976) Probability and
potentials. Mir, –M., 1976, p.436.
[16] Ermakov S.M. (1975) The Monte Carlo
method and related issues. Second edition,
expanded. Nauka, –M., 1975, p.471.
[17] Yelepov B.S., Kronberg A.A., Mikhailov
G.A., Sabelfeld K.K. (1980) Solving
boundary value problems by the Monte
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2024.23.18
M. G. Tastanov, A. A. Utemissova,
F. F. Maiyer, D. S. Kenzhebekova, N. M. Temirbekov
E-ISSN: 2224-2880
158
Volume 23, 2024
Carlo method. Nauka, Novosibirsk, 1980,
p.176.
[18] Mikhailov G.A. Voitishek A.V. (2006)
Numerical statistical modeling. Monte Carlo
methods. –M.: Publishing Center
«Academy» 2006, p.368.
[19] M. Tastanov, A. Utemissova, F. Maiyer, R.
Ysmagul. (2022) On a Solution to the
Stationary Problem of Two-Phase Filtration
by theMonte-Carlo Method.WSEAS
Transactions on Mathematics Volume 21,
2022, pp.386-394 DOI:
10.37394/23206.2022.21.88
Contribution of Individual Authors to the
Creation of a Scientific Article (Ghostwriting
Policy)
The authors equally contributed in the present
research, at all stages from the formulation of the
problem to the final findings and solution.
Sources of Funding for Research Presented in a
Scientific Article or Scientific Article Itself
No funding was received for conducting this study.
Conflict of Interest
The authors have no conflicts of interest to declare
that are relevant to the content of this article.
Creative Commons Attribution License 4.0
(Attribution 4.0 International, CC BY 4.0)
This article is published under the terms of the
Creative Commons Attribution License 4.0
https://creativecommons.org/licenses/by/4.0/deed.e
n_US
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2024.23.18
M. G. Tastanov, A. A. Utemissova,
F. F. Maiyer, D. S. Kenzhebekova, N. M. Temirbekov
E-ISSN: 2224-2880
159
Volume 23, 2024
