Построение трёхмерной модели геофильтрации флюида в многослойных пористых средах Текст научной статьи по специальности «Математика»

УДК 534.29:551.594.25

А.В. Шишеня, А.А. Афонин, АЛ. Сухинов ПОСТРОЕНИЕ ТРЁХМЕРНОЙ МОДЕЛИ ГЕОФИЛЬТРАЦИИ ФЛЮИДА В МНОГОСЛОЙНЫХ ПОРИСТЫХ СРЕДАХ

Целью данной работы является построение и исследование трёхмерной модели фильтрации флюида в многослойной среде. В результате работы были получены поля пьезометрической высоты и скорости фильтрации флюида в модельной области с двумя ти.

Уравнение Дарси; уравнение неразрывности; задачи фильтрации; трёхмерная визуа-.

A.V. Shishenya, A.A. Afonin, A.I. Sukhinov THREE-DIMENSIONAL MODEL OF FLUID FILTRATION IN A MULTI-LAYER POROUS MEDIA STRUCTURE

The main purpose of the work is to create and to investigate three-dimensional model of fluid filtration in a multi-layer environment. Fields of hydraulic head and filtration velocity in the model region had been received as a result.

Darcy equation; continuity equation; filtration; 3D visualization.

Подземные пространства, в которых протекают гидродинамические процессы, представляют собой, как правило, пористые или трещиновато-пористые тела, характеризующиеся сложной системой пор, каналов, трещин, по которым может происходить течение флюида - жидкости или газа. Течение через такие пористые тела, при котором сила трения флюида о скелет играет определяющую роль, назы-.

Численное решение задач фильтрации играет важную роль при строительстве дамб, скважин, ирригационных сооружений, а также при проектировании, постройке и эксплуатации гидротехнических и мелиоративных сооружений, в горном , . приводят к нелинейным уравнениям математической физики с весьма общими свойствами [1].

Целью данной работы является построение комплекса программ для нахождения трёхмерных полей пьезометрической высоты и скорости фильтрации флюида в многослойном грунте, а также средств, позволяющих визуализировать и изучать результаты расчетов.

В качестве модели рассмотрим процесс фильтрации флюида в грунте, состоящем из двух горизонтальных слоев с различными характеристиками. Будем , -ком флюида, например, рекой, а с другой - с поверхностью высачивания, напри, . . координатную ось z против вектора ускорения свободного падения, а оси x и у так, чтобы построенная система координат оказалась правой прямоугольной. В качестве расчетной области D возьмем прямоугольный параллелепипед (рис. 1).

i рунт 1 -го типа 1 раница грунтов

Рис. 1. Модельная область D с указанием типов граничных условий и положения

грунтов

Опишем модель геофильтрации следующими уравнениями [2]:

♦ уравнение неразрывности:

д

— {p™{p)s(p)) = -div(pv) + q, (1)

♦ зако н Дарси:

v = -K • grad(h), (2)

где K - в общем случае тензор 2-го ранга.

Для изотропной среды тензор K = k • E , где k - скалярная функция, E -

. K .

p - поровое давление (давление флюида), р - плотность носителя (флюида), m(p) - (

),

S (p) - функция водонасыщения, q - мощность источников, v - скорость фильтрации,

h = — + z (3)

Pg

( ),

K - тензор гидравлической проводимости.

Примем следующие упрощения для коэффициентов и функций, участвующих в уравнениях. Будем считать, что жидкость несжимаемая, а источники от:

р = const, q = 0.

Также допуская, что пористость слабо зависит от порового давления, положим, что пористость есть функция только координат:

т = т(х, у, г),

а функцию водонасыщения приближённо зададим полиномом [3]:

S (p ) = 3 a(x, y, z )■ p3 + 2 b(x, y, z )p2.

(4)

K

ki (, y, z) О О

О k 2 (x, y, z) О

О О k з (x, y, z)

V ” '“3V^^ “)J

(1), (2) ,

.

v = -K • grad(h)

д

Pm dt ( (p )) = -P■ div(x).

(1')

(2')

Выразим из (3) поровое давление и подставим в левую часть уравнения (2'). Тогда с учётом выражения (4) имеем:

р=р х -г X

Я х х - г У = И'г ЕР- [а^ y, [х (х - х )2 + Ьх y, г) • РЕ (х - г)).

Подставляя уравнение (1') в уравнение (2') и учитывая принятые упрощения, приходим к следующей системе уравнений:

К • gp • m( ^z ) • (a(( y, z У(pg(h - z )f+b( y, z ) • pg(h - z))=

/ / \' /

= ( ( y, z )h'x)x + ( ( y, z )h'y )y + (k3 ( y, z )tiz)z

v = -K • grad(h).

Поставим граничные условия в соответствии с описываемой моделью. 1. - . 2- :

dh ( )

— = Vi (y, z). dn

2. - ( ). 1- :

p = Ф2(x, У, z),

h = ^(w) + z.

Pg

(5)

(6)

(7)

3. Верхняя, нижняя, передняя, задняя стенки - непроницаемая граница

дк

0,

4. Граница двух пластов.

дп ^ = 0.

Рі = Р2 ,

Х1 = Х2, т.е. функция пьезометрической высоты непрерывна.

Выполним аппроксимацию уравнений модели.

Производная по временной координате.

Обозначим

I{х ^ г) = ЕР •m(х,У, г) (a(x,У, гУ Х(х - г))2 + b(x,У, г)' РЕ(х - г)),

тогда

кп+1 - кп

К- / (к У, 2 ) =----------------!г

Т

г,],к ■

Производные по пространственным координатам в центральных узлах в общем случае аппроксимируют следующим образом:

к + ,,, - к

к к - к

{кХ )х = а1 кг+— к— ~ к— кг-1-к

1г+1

к

а

{к2кУ ) У =

кг,, +1,к - кг,,,к кг,,,к - кг,,-1,к

« 2 , +1 ---------------2--------------а2 ,---------------^-------------

У

к ,/У кг,,,к+1 -кг,,,к

к3кх У 2 = а3

где

*3 к+1

а1г+1 =

ку

кг,,,к - кг,,,к-1

(9)

к2

*3 к

к2

1 Хг+1 1 [ к 1

ёх

\-1

кх X к1{к У, 2 У

2 ,+1

^ х Хг 1 V ^ ^ /

1_ - ёу

ку У к2 {x, У, Х)

а1к+1 =

1 2г+1

- (•

к

у

Л-1

к I кз(x, y, г)

Учитывая, что коэффициенты диффузии к1 (х, у, г), к2 (х, у, г), к3 (х, у, г) - 1-

грунтов, получаем:

а1г+1 = к1г+1, а1г = к1г-1,,,к ,

2 2

а2,+1 = ^2, + |,,,к , а2, = к2,-2,,,к ■

п

-1

а3к+1 Х3Х+1 ,,к , а3к Х3Х+1,,,к •

2 2

Учитывая полученные равенства в формулах (9), имеем:

(к1х; )х = к1 +1 .кх+1^к - к Лк - к1 1 ■кХы,к^ х~и ,к V 1 х) ь+?лк х2 1г - ?лк х.

(і ТІ кг,,+1, к - кг,—, к , кг, ,, к - кг, ,-1, к ,а,Л

\к2ку У у 2 2г,—+1,к ' 2 - к2г,—- 1,к ,2 , (9 )

2 ку, 7 2 к^у

хх 2 к, .к 1к—+■- ^ - х,. к 1 —-—1.

V 3 г; 3г,—,к+- к2 Зг,—,к-2 к2

П-

Рис. 2. Положение типов границ

Займемся аппроксимацией уравнений на границах. Обозначим через ух+, Ух-, Уу+ , У у-, Уг+, 7г- ™пы ГР аниц, ортогональные соответствующим коорди-

, ( . 2).

На границе с граничным условием 2-го рода (7) аппроксимация производных по пространственным координатам имеет вид

Ух+ : {ХХ У X 2 Х- | ({1ХХ У X ёх = Х-Г(кХ У+1—і - {кХ Уг-!,,Л = Хг[{кХ У+1,,,к - к1г — ,кРг,_

к - к к

, г+1,, ,к П1,],к 1г,, ,к

= Х +1 , ——-------------------------------------------—-—Ф

1г+—, !,к і 2 7. п,

ух- : {к1к'х У х 2 | {к1кХ У X Л = Х-Г{к1кХ ).+1,,,х - {к1к'х у.-1,,,Х 1 = тЧк1г,—х?г,,Х - {к1кХ Ц,

,—,х

к к - к

1г, —,Х » г, /,Х г-

й,—,Х - к1г-1,—,х----------------

Х 2 ^

к2

х

хх

х

у 1

2

у у 1 у

Уу+ : ^ 2 х'у )у * -1 | (к 2 х'у )уф = (к 2 ху )+1,к (к 2 ху ),7-Ьд] = х^( (к2 ху ^ - Члк^Лк

к ь Х1,—+1,к Х1^,к к21,—,к ?

' к2г,у+Ь,к , 2 / ^,./,к ,

2 Х х

у ь

2

Гу- :х2ху )' у * КТ ! х2х'у)'у *у = Х1 х2х'у 1+1,к - (к2х'у 1-1,к = ТТ к2- х2ху 1-1Л =

Гу--' 2 у/ хуу у/ ху

— 7

к2I,7,к . Х1,—,к - Х1,—-1,к

= - ?,—,к - к2^-Ь,к---------------------72-----■

1 > г,— ,к 2 l,J—,к 1 2

ху 2 х^

7г+ : (к3хгЬ */ (х3х0 гйЪ = -х~((зх'Х-,к+ь -(к3хг)1,-,к-Ь 1 = -х-Г(кзхг)1,-,к+ь -к31,—,к?1,—,к

, х—1,к х1,-,к к31,-,к

‘ к?; ; ь , 1

3—+2 х2 хг Гг,—,к:

7г-:(х3кь*х~ |(х3к)г&=(х3к)1,-,к+1 -(х3к)1,-,к-1)=хК^-.к?--(х3к

ко.. > Х- ■ 1 Х -71

31,-,к , I,-,к I, -,к-1

- к 3 — -1 —---------------------------

Хп+1 х1,—,к

Выпишем полученное сеточное уравнение в центральных узлах:

(.п к11+1,—,к k1l_Ь,—,к к21,—+\к к21 ,J_Ь,k к31,—,к +1 к31,—,к— ■>г,—,к + 2 + 2 + 2 + 2 + 2 +_ 2

ч к х2х х2х х2у х2у к2 к2

к11+1,—,к кИ——,к к 21,—+1,к

2 Х”+1 ^ Хп+1 2 хп+1 -

хх2 1+1,—,к хх2 l_1,J,k ху 1,—+1,к (10)

к 2 l,J_Ь,k к31,—,к+1

2 ХП+1 2 ХП+1 2 ХП+1 =

Х2 г—1,к х2 г,—,к+1 х2 ^

п

гп

1г,—,к х Х Х— ■

'ч

Исследуем разностную схему на устойчивость.

Канонический вид сеточных уравнений в общем случае даётся формулой

ДрМр)- Е^ч)у(ч)=-Г(р)VPе ®.

9^'(р )

ж,«)=кх2-+к<ч),

Хх\у\г

у

г

г

/

р {р )=ГШх {ру

А{р) = IХрУ+ X вХр.д)

ч^Ш'кр)

Коэффициенты теплопроводности являются положительными функциями как и функция /(х, у, г), поэтому

в{p, дУ> °

А{р)> 0,

я(р)=А{рУ- Xв{р,дУ=^ПкрУ > 0,

' (p)

-

, .

Исследуем разностную схему на консервативность.

Проинтегрируем исходное дифференциальное уравнение по всей области:

JJJк; ■ f)= JJJdiv( ■ grad (И ))dV .

Применяя к правой части формулу Остроградского-Гаусса, получаем

JJJ(К ■ f] = JJK ■ grad (h)dw

dD

(2)

Шк; ■f )V=-ffvdw

D dD

ИЛИ

ШК' ■ f )V = -ffvndw .

D dD

,

происходит только за счёт потока через стенки области.

, . Для этого понадобится вспомогательное соотношение:

h - h N2 -1

k 1+1 + х

N- +_ Ъ ^

2

h

k hi+- - hi - k hi - hi-i

. 1 7 2 . -

i=N +1\, i+ 2 hx

1 hX у

Их - k„ 1

h - h

N2 N2 -1

= X

i=N1

2 hi+1 - hi - hi - hi-1

k -

i+— И

V. 2 '/.

--k -

, , hN_ hN_ -1 ,

hx + k — —1------------------------------------------1-k —

X N_ -_ h N2 +_

1 2 x 2 2

h

И - h

"n2+1 "n2

h

= X

i=N1

N2 h - h

A/i+1 ni

k -

= X

i+1 h

V. 2 x у

N2 h - h

k i+- i

hx - X2

N2 h - h k hi hi--

k 1 7.2

i=N-V i+2 h;t у

i=N1

n, -1

i— h

\ 2 nx J

h - h

N2+1 N2

, , hN_ hN_ -1 ,

h. + k --------------------k -

Nl-2 hx N2 + і h

hx - X

k hi+1- hi

i=N- -^ i+ 2 hi у

h - h

N2 +1 "n.

= k -

N2 +22

h - h

N2+1 N2

hN - hN -1 hN - hN -1

- k N N--1 + k 1 1-1 - k

/V і І /V і /V і

N-— — N-— — n, +—

1 2 x 1 2 x 2 2

h - h

N2+1 N2

= О.

D

D

D

2

h

h

X

X

Обозначим

Х

дк

дп

дк

дп

тогда

к - к

Х "ц-и п» - х

/V 1 /V 1

»1+2 кх Ы1-~2

к - к

пщ пщ-- + 'у

к2 £

>1.. І-КЇ.4.

^-1 (к Х+1 - Х - Х Х - кг -1

+1 к2

V 2 х

1 к2

2 х у

к+

, кN2 +1 Х»2 ; Х»2 2-1 / »

+Х .1—^Тг ^- К 1 , = К 1 1

к - к

2

к2

к

N +—

12

к

- + Ц +

+

2

£

і к+ - кг кг - кг ^

Х г+1 г - х г 1-1

.+1 к2

г=»1+1 V 2 "х

1 к2

2 х у

7 7 к»2 к»2-1 _

кх + ^2 - Х 1 ---------------)------ = ^1 + ^2-

»2-т Х

То есть для сеточного уравнения выполняются те же законы сохранения, что и для непрерывного, значит, полученная разностная схема консервативна.

Найдём порядок аппроксимации исходного дифференциального уравнения .

Хорошо известны оценки вида:

,, кп+ - кп

к'----------------------

' к

0{к,),

{Хкх ух -(хху ) -

{ХК )х -

^ к+ , - кг , кг , - кг ^ ^

г+1, — г, — и г, — г-1, —

----------2-----------Х 1 ------------------------

Х 1

г+—,, к

V 2 ,1х

кх ;

\

кг+1 — - кг - кг - - кг ^ -

Х г+1, —___________і-____х г— г-1, —

г'+-,,

V 2

к2

1

г-Т, —

Х кг+1,— ^ к,— - х

к2

кг,, - кг-1

1

и—, — V 2

х /

Л

— "г-1, —

к

г-?'

к; ,

= 0(к’), = Ф; ^ = о(к?).

Рассмотрим теперь нелинейный коэффициент в левой части уравнения.

/п*1 = /п + 0(А,),

а так как

к'=+о(к), к

то

к'г+ =

кп+1 - кп кп+1 - кп

к

/п +■

к

О (к') + /п0(к,)+ о(),

п+1 п

'-'^п+1 - к -— /п + О (к').

к'

к'Г+ =

г

7:

,

уравнения полученной разностной схемой составляет о(к{ + И2Х + И2у + И22).

В результате аппроксимации систем линейных алгебраических уравнений была получена СЛАУ большой размерности. Матрица СЛАУ является симметричной, положительно определённой и имеет разреженную структуру, поэтому в качестве экономичного метода решения данной системы решено было использовать адаптивный попеременно-треугольный метод скорейшего спуска [4].

Решение исходной системы

Au = f

заменяется итерационным процессом

yn+1 - yn ВпУ---------------------^ + Ay” = f,

Tn+1

У0 = иО.

где в качестве матрицы переобуславливателя берётся следующая матрица:

Bn =( +®nA +«nA2 )> A = A1 + A2 , A1 = A2 ’

=

y

Tn+1 = 2^n

A2 ,yn

где A1 и A2 соответственно верхне- и нижнетреугольные матрицы, составленные из элементов матрицы A .

Результаты

По результатам работы были составлены параллельная и последовательная .

распределённых вычислений под управлением MPI.

Результаты моделирования представляют собой трёхмерные сеточные функции напора и скоростей движения флюида и требуют адекватной визуализации. Эффективным средством представления трёхмерных объектов является программно-аппаратный комплекс трёхмерной визуализации. Программное обеспечение Visual Instruments, установленное на комплексе представляет собой постпроцессор и предназначено для удобного представления трёхмерных сеточных функций и

. Visual Instruments

функции построения срезов, изоповерхностей, линий тока, траекторий движения ,

способна выводить все вышеперечисленные объекты в виде, воспринимаемом человеком как трёхмерный. Комплекс способен читать некоторые промышленные форматы данных, такие, как VIB и encase. В связи с этим были изучены структуры данных форматов и реализованы специальные процедуры преобразования внут-

VIB encase.

а

Рис. 3. Изоповерхности пьезометрической высоты (слева) и линии тока флюида и изолинии пьезометрической высоты на поверхности области (справа)

Расчет был проведён для модельной задачи. За отсутствием реальных данных о свойствах почв для грунта первого типа были взяты характеристики

М / \

k1 = к2 = к3 = 1000------, m = 0,3, S(p) = 0,2, для грунта второго типа

мес

к1 = к2 = к3 = 700^-, m = 0,03 , S(p) = 0,2, т.е. грунт второго типа малопрони-

мес

. . взяты шаги hx = hy = hz = 1, ht = 0,1, а количество узлов n = m = к = 8 . Качест-

. 3.

БИБЛИОГРЛФИЧЕСКИЙ СПИСОК

1. Афо нин АЛ. Задачи фильтрации жидкостей в пор истых средах со свободной границей.

- Ростов-на-Дону. 2008. - С. 6.

2. Величко О.М., Горев В.В., Горев КВ., Дерюгин ЮМ., Зеленский Д.К., Панов А.И., АЛ. Савельев, Селин В.В., Субботин АТ., Тихомиров Б.П., Чекалин AM. Математическое моделирование. 2001. Т. 13. №2. Методика «Триада» численного решения трёхмерных задач переноса примесей подземными водами.

3. Принцыпар Г.В., Сухинов AM. / Математическое моделирование задач двухфазной фильтрации па кластере распределённых вычислений. 2008. - С. 17.

4. Коновалов AM. // Сибирский математический журнал, май-июнь 2002. Т. 43, - №3.

- С. 552-572. К теории попеременно-треугольного итерационного метода.

Шишеня Александр Владимирович

Технологический институт федерального государственного образовательного учреждения высшего профессионального образования «Южный федеральный университет» в г. Таганроге.

E-mail: primat-55-alex@yandex.ru.

347928, г. Таганрог, пер. Некрасовский, 44.

Тел.: 8(928)322-282, 8(908)176-18-37.

Кафедра высшей математики; студент.

Shishenya Alexander Vladimirovich

Taganrog Institute of Technology - Federal State-Owned Educational Establishment of Higher Vocational Education “Southern Federal University”.

E-mail: primat-55-alex@yandex.ru.

44, Nekrasovskiy, Taganrog, 34792В, Russia.

Phone: 8(928)322-282, 8(908)176-18-37.

The Department of Higher Mathematics; student.

Афонин Анатолий Андреевич

Технологический институт федерального государственного образовательного учреждения высшего профессионального образования «Южный федеральный университет» в г. Таганроге.

E-mail: afonin_920@gmail.com.

347928, . , . , 44.

Тел.: 8(8634)371-606

Кафедра высшей математики; профессор.

Afonin Anatoliy Andreevicy

Taganrog Institute of Technology - Federal State-Owned Educational Establishment of Higher Vocational Education “Southern Federal University”.

E-mail: afonin_920@gmail.com.

44, Nekrasovskiy, Taganrog, 347928, Russia.

Phone: 8(8634)371-606

The Department of Higher Mathematics; professor.

Сухинов Александр Иванович

Технологический институт федерального государственного образовательного учреждения высшего профессионального образования «Южный федеральный университет» в г. Таганроге.

E-mail: sukhinov@gmail.com.

347928, . , . , 44.

Тел.: 8(8634)310-599; 7(928)102-11-06.

Руководитель ТТИ ЮФУ; профессор.

Sukhinov Alexander Ivanovich

Taganrog Institute of Technology - Federal State-Owned Educational Establishment of Higher Vocational Education “Southern Federal University”.

E-mail: sukhinov@gmail.com.

44, Nekrasovskiy, Taganrog, 347928, Russia.

Phone: 8(8634)310-599; 7(928)102-11-06.

Chief of TIT SFedU; professor.

УДК 534.29:551.594.25

А.А. Афонин, АЛ. Сухинов

МАТЕМАТИЧЕСКИЕ МОДЕЛИ ГЕОФИЛЬТРАЦИИ И ГЕОМИГРАЦИИ В ПОРИСТЫХ СРЕДАХ, ОБЛАДАЮЩИХ ФРАКТАЛЬНОЙ СТРУКТУРОЙ

В данной статье представлены двумерные линейные модели, вытекающие из обобщенного уравнения Бусстеска, описывающего геофильтрацию в почвах с фрактальной структурой. Кроме того представлены математические модели геомиграции загрязнений в грунтовых водах в классической постановке, а также в почвах, обладающих фракталь-.

Геофильтрация; уравнение Буссинеска; фрактальные структуры; геомиграция, грунтовые воды; почвы.