Numerical research of the Barenblatt - zheltov - Kochina stochastic model Текст научной статьи по специальности «Математика»

MSC 60H30

DOI: 10.14529/mmp 160211

NUMERICAL RESEARCH OF THE BARENBLATT -ZHELTOV - KOCHINA STOCHASTIC MODEL

S.I. Kadchenko, Nosov Magnitogorsk State Technical University, Magnitogorsk; South Ural State University, Chelyabinsk, Russian Federation, [email protected], E.A. Soldatova, South Ural State University, Chelyabinsk, Russian Federation, [email protected],

S.A. Zagrebina, South Ural State University, Chelyabinsk, Russian Federation, [email protected]

At present, investigations of Sobolev-type models are actively developing. In the solution of applied problems the results allowing to get their numerical solutions are very significant. In the article the algorithm for numerical solving of the initial boundary value problem is developed. The problem describes the pressure distribution of the homogeneous fluid in the horizontal layer in the circle. The layer is opened by a vertical well of a small radius. In our research we suppose that random disturbing loads have an influence on the fluid. The problem was solved under two assumptions. Firstly, we suppose that an unstable fluid flow is axially symmetric, and secondly, that in initial moment the pressure in the layer is constant. After the process of the discretization we modify the original model to the Cauchy problem for the system of ordinary differential equations. For the numerical solution we use algorithms based on explicit one-step formulas of the Runge - Kutta type with the seventh-order accuracy and with the selection of the integration step. We also use the scheme of the eighth-order accuracy to evaluate the calculation accuracy on each steps of time. According to the results of this control, we choose the time-step. A lot of numerical experiments have shown high numerical efficiency of the algorithm that we use to solve the investigated initial-boundary problem.

Keywords: stochastic Sobolev type equation; numerical solution; Barenblatt - Zheltova - Kochina model; Cauchy problem.

Introduction. Let G C Rd be a bounded domain with a boundary dG of class C. Consider in the cylinder G x R+ the Barenblatt - Zheltova - Kochina equation [1]

(A - A)ut = aAu + f, (1)

with the Dirichlet condition

u(x, t) = 0, (x, t) e dG x R+. (2)

Model (1), (2) simulates the pressure dynamics of the fluid filtered in fractured porous

f

[3] and processes of the solid-to-fluid thermal conductivity in the environment with two temperatures [4]. In equation (1) a and A are real parameters, that characterize the

environment; the parameter a e R+, and the parameter A can also take the negative

f

find the random process u = u(x,t) that satisfies (1), (2) with initial Cauchy condition

u(x, 0) = u0(x). (3)

It must be noted that the problem (1), (2) in corresponding Hilbert spaces can be reduced to the linear stochastic Sobolev type equation

Ldu = Mudt + NdW.

You can find more historical aspects of the stochastic equations theory (including the stochastic Sobolev type equations research), for example, in [6].

1. Homogeneous Fluid Filtration. Consider the algorithm of construction of numerical solution to (1), (2). As an example we consider unstable filtration of the homogeneous fluid in the horizontal layer. The layer is opened by a vertical well of a small radius. At the initial moment pressure of fluid in the layer is constant and is equal to p0. Assume that a random load begins to influence the fluid. Let the domain G C R2 be presented by disk with the radius r0 and center on the axis of the well. Introduce a polar system of coordinates and assume that the fluid flow is axially symmetric. In this case we can find the fluid pressure in a fissure if we solve the initial-boundary value problem [1]

(A — Ar)pt = aArp + A sin(^t), , ,

\p(0,t)\ < p(ro,t) = 0, p(r, 0) = po, ^

d2 1 d

where Ar = —— +— —, A E R, A = 1/^, ^ is a specific characteristics of the fractured

dr2 r dr

ground, a = x/n> X is the piezoconductivity coefficient of the fractured ground.

Then we need to modify (4) to the system of ordinary differential equations. So we make a discretization of the domain G. We work on a segment [0, r0] with the equally

spaced grid rj = (j — 1)hr, j = 1, J, hr = r0/(J — 1) J is the number of discretization

p

j

Op j +jp

' M)j = hpH apjiUi + o(hnu-p), (5)

,dr-P,

i=j-jl

where ji(jp) is the amount of points to the left (right) from the j-th point in which the derivative is calculated, i is numbers of points that we use in the approximation of the derivative nu = ji + jp + 1 is the amount of points that we use in the derivative

approximation, nu < J. After we use the expansion of a function p(r, t) in a Taylor series rj

apji- We have already made a special program in Maple Math Software to calculate apji for p E N-order derivative and for j,i = 1, J. This program gives an opportunity to use (5) with difference approximation errors O^h^-1^ during the process of numerical simulations.

After we have done the discretization of initial-boundary conditions (4) we obtain a Cauchy problem for the system of ordinary differential equations for the function Uj (t). The system we will write in the matrix form

Mdm=PU (i)+

U (0) = Uo.

Where M =

пл

0

Л - &2,2

b2,3

2J

bj~ 1,1 bj-1,2 bj_ i

Л - br

j 1,j

Un

UOi Un2

Uoj _

v uoj y

p

V

b2,1

bj -1,1 О

О

О

b2,2

bj -1,2 О

( U1(t) ^

U2 (t)

U (t) =

О

b2,3 bj-1,3

О

/

0

b2,J

bj -1,j

1

V

Uj-1(t) Uj (t) y

О

A sin(wt)

W

\

A sin(wt) О

/

■Uj is values of pressure in the j-th point of digitalization. If we don't use the point in

a2ji a1ji

where the first

derrivativa approximation, then j = 0 otherwise j = +

index j indicates the number of the point in which we make the derivative approximation and the second index i indicates the number of point, that takes part in the derivative approximation.

It's easy to show that the determinant of the matrix P is not equal to 0. In rare cases when we set A it may be that the determinant of the matrix M equals to 0. In this case we need to change the amount of points that take part in derivatives approximation. It will

P

of the determinant of matrix M. Therefore, we assume that det(M) = 0.

M

for the first-order differential equations in the matrix form

dU (t) dt

= F(t,U(t)), О <t < t0 U (0) = Uo,

(7)

where F (t, U (t)) = M -1 ( PU (t) + W) •

For the numerical solution of (7) we use algorithms that are based on explicit one-step formulas of the Eunge - Kutta type with the seventh order accuracy [7]:

13

Un+1 — Un + ^J Pn,iki

i=1

13 13

ki = hF[tn + aht, фп +J2 fii,j kj, An + Y, fii,j kj) ■ v j=1 j=1 J

(8)

Where ht is a integration step in Eunge - Kutta method,

2 115 15

«1 = 0, a2 = 27, «3 = 9, «4 = ^, a5 = 12' a = 2' aj = 6'

1 2 1

«8 = 7, «9 = -, «10 = -, «11 = 1, «12 = 0, «13 = 1, 0 3 3

1 8 '

2111

@2,1 = 27, в3,1 = 30, в3,2 = 12, в4,1 = 24, в4,2 = 0' в4,3

5 25 25 1

e5,1 = 12, вб,2 = 0' вб,3 = —10, вб,4 = 10, вб,1 = 20, вб,2

1

0

0

3

0

0

1

1

1 1 25 125 65

A.,4 = 4 , ^6,5 = 5 , ^7,1 = - ^ , ^7,2 = ^ = 0 , ^ = —, ^ = - 27 ' 125 31 61 2

@7,6 = -54- ' ^8,1 = 300 ' ^8,2 = ^8,3 = ^8,4 = 0 ' ^8,5 = 225 ' B8,6 = — 9 '

13 _ 23 _ 704 _ 107

B8,7 = 900 ' BM =2 ' B9,2 = B9,3 =0 ' B9,4 = 108' B9,5 = 15 ' B9,6 = —9"'

67 91 23

P9,7 = 90 ' P9,8 = 3 ' ^10,1 = — 108 ' P1°>2 = P10,3 = 0 ' ^10,4 = '

B = 976 B =311 B = 19 B =17 B = 1

P 10,5 = —135 ' B1°,6 = ' B10,7 = — 60 ' B1°,8 = -Q ' B10,9 = —12 '

_ 2383 _ 341 4496 301

1 = 4100' ,2 = B11,3 = 0 ' ,4 = — 164 ' ,5 = 1025 ' B11,6 = — ' = 2133 B = 45 B = 45 B = 18

7 = 4100 ' ,8 = 82 ' B11,9 = 164 ' ,10 = — 41 '

3 6 3

B12,1 = 205 ' B12,2 = B12,3 = B12,4 = B11,5 = 0 ' B12,6 = — ^ ' ^12,7 = — 2^ '

a _ 3 _ 3 _6 1777

B12,8 = —7T' B12,9 — 77 ' B12,10 — 77' B12,11 = 0 ' B

12,8 = — 41 ' B12,9 = 41 ' B12,10 = 4! ' B12,11 = 0 ' B13,1 4100

341 4496 289

B13,2 = P13,3 = 0 ' B13,4 = — ' B13,5 = lnor ' B13,6

164' ' ' 1025' ' ' 82 ' 2193 51 33 12

B13,7 = — 4100 ' B13,8 = 82 ' B13,9 = 164 ' B13,10 = 41 ' B13,11 = 0 ' B13,12 = 1 ,

41 34 9

P7,1 = 840 ' P7,2 = P7,3 = P7,4 = P7,5 = 0 , P7,6 = ' P7,7 = P7,8 = 35 '

9 41

P7,9 = P7,10 = 280 ' P7,n = 840 ' P7,12 = P7,13 = 0

When we use explicit methods to solve the stiff problem we would like to note that the integration step is limited not only by calculation accuracy but also by the stability of the numerical method [8]. We can find the local error 5n of the method with the seventh order accuracy by using the next formula [7]

13

=YJ('P8,i — P7,i)ki' (9)

i=1

where coefficients P8i are equals

34 9

P8,1 = P8,2 = P7,3 = P7,4 = P7,5 = 0 ' P7,6 = 777' P8,7 = P8,8 = 77 '

105 35

9 41

P8,9 = P8,10 = 280 P8,n = P8,12 = 0 ' P8,13 = 840'

If you would like to evaluate the control of the calculation accuracy you can use the next inequality

\\Sn\\< s.

Here || • || is a norm in RJ, J is the dimension of the matrix column F(t, U(t)), £ is the necessary calculation accuracy. In our case 6n = O^ht^- Therefore integration step hac according to accuracy is calculated with the formula [8]

In the case q < 1 we need to repeat the calculation that we received at the step ht and

hac

time.

2. Results of the Calculation Experiment. The algorithm for solving of problem (4) that was described in previous paragraph, was implemented in a special program in Maple Math Software. This program gives us an opportunity to conduct numerical experiments that find the values of the pressure in breaches under a specific condition. This condition describes random loads, that have difference influence on the seam consumption. We use the random number generator to choose the constant A E R for function that determines random loads for the seam consumption. Estimates of experiments show that in equation (4) the parameter n f°r different earth materials take the value between very wide limits n = 10"4 m2 ^ 106 m2, and the coefficient x = 0,1 m2/s ^ 1 m2/s.

Numerous calculations have shown that the developed algorithm of the solution allows to find the numerical solution of the problem with high computational efficiency.

References

1. Barenblatt G.I., Zheltov Yu.P., Kochina I.N. Basic Concepts in the Theory of Seepage of Homogeneous Fluids in Fissurized Rocks. Journal of Applied Mathematics and Mechanics, 1960, vol. 24, no. 5, pp. 1286-1303._

Fig 1. Dependence of pressure p = p(t*,r) in a fissure at moments of time t* = 1, 27s, t* = 2,61s t* = 4s

p

p(t, r) t r

2. Zagrebina S.A., Soldatova E.A. Linear Sobolev Type Equations with Relatively p-Bounded Operators and Additive White Noise. Bulletin of the Irkutsk State University. Series: Mathematics, 2013, vol. 6, no. 1, pp. 20-34. (in Russian)

3. Hallaire M. On a Theory of Moisture-Transfer. Inst. Rech. Agronom., 1964, no. 3, pp. 60-72.

4. Chen P.J., Gurtin M.E. On a Theory of Heat Conduction Involving Two Temperatures. Z. Angew. Math. Phys., 1968, vol. 19, pp. 614-627.

5. Sviridyuk G.A. Couchy Problem for the Linear Singular Sobolev Type Equation. Differential Equations, 1987, vol. 23, no. 12, pp. 2169-2171.

6. Sviridyuk G.A., Manakova N.A. The Dynamical Models of Sobolev Type with Showalter - Sidorov Condition and Additive "Noise". Bulletin of the South Ural State University. Series: Mathematical Modelling, Programming and Computer Software, 2014, vol. 7, no. 1, pp. 90-103. (in Russian) DOI: 10.14529/mmpl40108.

7. Novikov E.A., Shornikov Yu.V. Control of the Stability for the Fulberg Method with Seventh-Order Accuracy. Computational Technologies, 2006, vol. 11, no. 4, pp. 65-72.

8. Novikov E.A. Algorithms of Variable Structure Development for the Solving of the Stiff Problems. Krasnoyarsk, 2014. 123 p. (in Russian)

Received January 20, 2016

УДК 517.9 DOI: 10.14529/mmpl60211

ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ СТОХАСТИЧЕСКОЙ МОДЕЛИ БАРЕНБЛАТТА - ЖЕЛТОВА - КОЧИНОЙ

С.И. Кадченко, Е.А. Солдатова, С.А. Загребина

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

Ключевые слова: стохастическое уравнение соболевского типа; численное решение; модель Баренблатта - Желтова - Кочиной; задача Коши.

Литература

1. Баренблатт, Г.И. Об основных представлениях теории фильтрации в трещиноватых средах / Г.И. Баренблатт, Ю.П. Желтов, И.Н. Кочина // Прикладная математика и механика. - 1960. - Т. 24, № 5. - С. 58-73.

2. Загребина, С.А. Линейные уравнения соболевского типа с относительно p ограниченными операторами и аддитивным белым шумом / С.А. Загребина, Е.А. Сол-датова // Известия Иркутского гос. ун-та. Серия: Математика. - 2013. - Т. 6, № 1. -С. 20-34.

3. Hallaire, М. On a Theory of Moisture-Transfer / M. Hallair // Inst. Rech. Agronom. - 1964. -№ 3. - pp. 60-72.

4. Chen, P.J. On a Theory of Heat Conduction Involving Two Temperatures / P.J. Chen, M.E. Gurtin // Z. Angew. Math. Phys. - 1968. - V. 19. - P. 614-627.

5. Свиридюк, P.A. Задача Коши для линейного сингулярного уравнения типа Соболева / P.A. Свиридюк // Дифференциальные уравнения. - 1987. - Т. 23, № 12. - С. 2169-2171.

6. Sviridyuk, G.A. The Dynamical Models of Sobolev Type with Showalter - Sidorov Condition and Additive «Noise» / Sviridyuk G.A., Manakova N.A.// Вестник ЮУрГУ. Серия: Математическое моделирование и программирование. - 2014. - Т. 7, № 1. - С. 90-103.

7. Новиков, Е.А. Контроль устойчивости метода Фельберга седьмого порядка точности / Е.А. Новиков, Ю.В. Шорников // Вычислительные технологии. - 2006. - Т. 11, № 4. -С. 65-72.

8. Новиков, Е.А. Разработка алгоритмов переменной структуры для решения жестких за-Дс1Ч« ДИ С« *«« Ту ТI « физ.-мат. наук / Е.А. Новиков. - Красноярск, 2014. - 123 с.

Сергей Иванович Кадченко, доктор физико-математических наук, кафедра «Прикладная математика и информатика:», Магнитогорский ГОСуДЯрСТВбННЫИ Т6ХНИЧ6-ский университет им. Г.И. Носова (г. Магнитогорск, Российская Федерация); кафедра «Уравнения математической физики», Южно-Уральский государственный университет (г. Челябинск, Российская Федерация), [email protected].

Екатерина Александровна Солдатова, ассистент, кафедра «Уравнения математической физики», Южно-Уральский государственный университет (г. Челябинск, Российская Федерация), [email protected].

Софья Александровна Загребина, доктор физико-математических наук, кафедра

«»

ныи университет (г. Челябинск, Российская Федерация), [email protected].

Поступила в редакцию 20 января 2016 г.