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

10. Belotserkovskiy O.M. Turbulentnost': novye podkhody [The theory of difference schemes]. Moscow: Nauka, 2003, 286 p.

11. Samarskiy A.A., Nikolaev E.S. Metody resheniya setochnykh uravneniy [Methods for solving grid equations]. Moscow: Nauka, 1978, 592 p.

12. Konovalov A.N. K teorii poperemenno-treugol'nogo iteratsionnogo metoda [The theory of alternating triangular iterative method], Sibirskiy matematicheskiy zhurnal [Siberian mathematical journal], 2002, No. 43:3, pp. 552-572.

13. Sukhinov A.I., Chistyakov A.E. Adaptivnyy modifitsirovannyy poperemenno-treugol'nyy iteratsionnyy metod dlya resheniya setochnykh uravneniy s nesamosopryazhennym operatorom [Adaptive modified alternating triangular iterative method for solving grid equations with non-self-adjoint operator], Matematicheskoe modelirovanie [Mathematical modeling], 2012, Vol. 24, No. 1, pp. 3-21.

14. Voevodin V.V., Voevodin Vl.V. Parallel'nye vychisleniya [Parallel computing]. St. Petersburg: BKhV-Peterburg, 2004, pp. 134-154.

15. Sukhinov A.A. Rekonstruktsiya ekologicheskoy katastrofy v Azovskom more na osnove matematicheskikh modeley // Matematicheskoe modelirovanie [Mathematical modeling], 2008, Vol. 20, No. 6, pp. 15-22.

Статью рекомендовал к опубликованию д.ф.-м.н., профессор А.А. Илюхин.

Лапин Дмитрий Вадимович - Южный федеральный университет; e-mail: dmitri.lapin@gmail.com; 347928, г. Таганрог, пер. Некрасовский, 44; тел.: 88634371606; кафедра математического обеспечения суперкомпьютеров.

Чистяков Александр Евгеньевич - e-mail: cheese_05@mail.ru; кафедра математического обеспечения суперкомпьютеров; доцент.

Сухинов Антон Александрович - e-mail: soukhinov@gmail.com. 347932, г. Таганрог, ул. Пархоменко, 58/1; тел.: 89296796206; научный сотрудник.

Lapin Dmitry Vadimovich - Southern Federal University; e-mail: dmitri.lapin@gmail.com; 44, Nekrasovskiy, Taganrog, 347928, Russia; phone: +78634371606; the department of mathematical software supercomputers.

Chistyakov Alexander Evgenjevich - e-mail: cheese_05@mail.ru; the department of mathematical software supercomputers; associate professor.

Sukhinov Anton Alexandrovich - e-mail: soukhinov@gmail.com; 58/1, Parkhomenko, Taganrog, 347932, Russia; phone: +79296796206; research scientist.

УДК 519.684.6

М.Д. Чекина

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ЗАДАЧ ГЕОФИЛЬТРАЦИИ В ПОЧВОГРУНТАХ С ФРАКТАЛЬНОЙ СТРУКТУРОЙ НА МНОГОПРОЦЕССОРНЫХ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМАХ*

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

модели на основе уравнения Букингема-Ричардса, в котором оператор дифференцирования по

*

Работа выполнена при поддержке Министерства образования и науки РФ по Соглашению с уникальным идентификатором КГМЕП57814Х0006.

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

Геофильтрация; аномальная диффузия; суперЭВМ; численный метод; МПТМ; дробные производные; фрактальные структуры.

M.D. Chekina

MATHEMATICAL MODELING OF GEOFILTRATION IN SOILS WITH FRACTAL STRUCTURE ON MULTIPROCESSOR COMPUTER SYSTEMS

When modeling the geofiltration in soils it is necessary to consider the fractal structure of the soil. The complex geometry of the pores and capillaries provides this structure. Anomalous diffusion describes the fractal diffusion mathematically, which is described by equations in partial fractional derivatives. This paper describes the construction of mathematical models based on equations of Buckingham-Richards, in which the differentiation operator time was replaced by the private fractional derivative Riemann-Louiville. For the continuous model was obtained the discrete analogue using the integro-interpolation method. In the numerical solution of problems of this type is to process large amounts of data, hence there is the need to use in the calculations the super computer To solve this problem was developed parallel implementation of the modified alternating triangular method, which allowed to increase several times the speed of the program complex (PC). High performance PC needed for rapid simulation results. For example, this will allow the rise of the groundwater table in real time, and thus to minimize damage from flooding.

Gelfiltration; anomalous diffusion; supercomputers; numerical method; the modified alternating triangular method; fractional derivatives; fractal patterns.

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

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

При численном решении задач такого типа возникает необходимость в обработке большого количества данных, поэтому целесообразно подобные расчеты проводить на суперЭВМ.

Моделируется процесс просачивания воды в почву. Растекание жидкости по порам и капиллярам почвы происходит фрактально, как показано на рис. 1.

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

Рис. 1. Фрактальная картина растекания жидкости в пористой среде

Непрерывная модель. Для моделирования процессов фильтрации в пористых средах используется уравнение Букингема-Ричардса [1]

С• Ь*^ = {к(И)И')х' + (А (/;)/;;) ' + {к(И)И')_' О. (1)

В данном уравнении x, у - горизонтальные координаты; г - вертикальная; в - влажность; к - потенциал почвенной влаги; к (к) - коэффициент влагопровод-

ности; - толщина распределения источников и потерь; ¥ - источники (приток, осадки и др.); Q - потери (сток, испарение и др.); ъ* 3(?) = -

1 \rndT -

' ' ^3 (1) Г(1 -д) { -г)*

регуляризованная дробная производная Римана-Луивилля порядка х (0 < * < 1 )• Зависимость коэффициента влагопроводности к ( к ) от потенциала к имеет

вид

^ак

к < 0

|к (к ) = к} [к (к ) = к,, к > 0, где ку - коэффициент фильтрации; X - эмпирический параметр.

Для влажности в будет справедлива следующая зависимость от потенциала к:

\в(к) = (в,-вГ)е(к + вг, к <0

|в( к ) = в,, к > 0, где в - максимальная влажность; в - минимальная влажность.

5 Г

Дополним уравнение (1) граничными условиями, для этого воспользуемся законом Дарси:

V = -к^аб ( к ), (2)

где V - скорость фильтрации флюида.

В случае непроницаемой стенки должно выполняться требование отсутствия

потока флюида через эту границу V ■ П = 0 . Это требование в совокупности с законом Дарси (2) дает следующее граничное условие:

д±-= 0. Эп

В случае, если через границу задан поток жидкости 5", то, исходя из тех же соображений, будем иметь следующее граничное условие:

д-К = 8.

дп

Во многих случаях логично предположить, что поток жидкости зависит от потенциала линейно, тогда получим следующее граничное условие:

дН , о — = ак + В. дп н

Разностная схема. Построим разностный аналог уравнения (1). Для этого покроем расчетную область равномерной прямоугольной сеткой с шагами кх, ку, по соответствующим трем пространственным направлениям:

3 = {(х1,у],гк ) ,хг=г ■ К,УГ] ■ Иу,гк = к ■ кг,1 = 1..п1а/ = 1..п2, к = 1..П3}, где И , Н , Н .

П

п

Обозначим внутренние узлы сетки через Ш,

з = {(х,ур2к), х = г ■ К' У/ = ■ ■ Ну,?к=к ■ Нг,г= 1..п -1, ■ = 1..п2 - 1, к = 1..пъ -1},

а границу - у,

у = {з \ з}.

Проведем аппроксимацию уравнения (1). Дробную производную ^Н

заменим дискретным аналогом [2]:

Цх К = -

1

п

-У^1-* - (1-*) И

I /Л 1п-5+1 1п-$

Ч ,г,],к

М Г ( 2

Заменим частные производные в уравнении (1) на их конечно-разностные аналоги, которые были получены при помощи интегроинтерполяционного метода.

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

( \

(К )х = ^ |\{кн\)'хах=±

К:х 1 К

кК

- кК

1

: к — (Нп+1 - Нп+1) - к —(Нп+1 - Нп+1 ) = ■к. 1 , и2 (Н/+1,¡к Нг,¡к) к 1 , /„2 (Нг, ¡к Н1-1,_¡к)

к , Нп+1

,к г+1,]к

■ I . , 7„2 ^ г,¡,1

Г

Л

к | + к |

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

V 2 2 У

Нп+1 -к 1 Нп+1 ¡к

¡,]к г-1,¡к г-^ ,к

Н

2

(кИу )У = £ ^'у )У^у=тт

у Ч

к.

кК ,

- кН\

у и. 1 ' + 2

у. 1 '- 2 У

к —(Нп+1 - Нп+1) - к — (Нп+1 - Нп+1 V ■ ки+1, Н2 (',./+1,к ) \¡+1, Н2 1,к)

Л.М* Н2

V 2 Л*у

f

n+1

к , к

\

к 1 + к |

i, /+1,к i, j-—,k V ,j 2' ,j 2' У

кй+1 -к 1 hn+\к

i,j,к i'j-1'к i'J-1'к 'J 2

h>2

(кк * )2 - /Г |Гг+1(кк )*ö2 --!

h,Jz, 1

' 2

;,n+1 7„n+1

:+1

f

кк'

-кк

'+2

1

2 У

- к — (hf+1 -hn+1)-к — (hn+1 -hn+1 ) =

к h

; U . 1

n+1

i^^1 i'j'k

+1

2 "2

Л

к 1 + к 1 hn+1 -к 1 hj 1

V 'j' 2 'j' 2 У _J 2_

к2

Запишем аппроксимацию уравнения (1) по неявной схеме в центральных узлах расчетной сетки:

1 п

С__1-У(/-1-х -)к* =

С ' "-+ Ч' '¿¿к

г (2 - X) %

1

h

к 1 h

V

'i+Ик

Л

Л

к 1 , + к 1 , hi, j к к 1 ,hi-1' у,к

i+—, у,к i—, j .к i—, j ,к

/

+ -

hl

V 2 /

к 1 , hi, j+1'к V 2

2 У Л

к y + к y

i, j—,к i, j—,к V 2 2 у

hi, jt к 1 , hi, j-1,к i, j—,к 2

+

+ -

h

с

1 hi,j,k +1

V 2

Л

к j + к j

i, j ,к+— i, /,к— V 2 2 У

Л

* 1 hi,j ,к-1

i,j^+

+К, к + öi,j ,к •

2 У

Рассмотрим аппроксимацию выражений (к ( к) к') ', (к ( к) к') ' и

(к (к) к') ' на границах расчетной области. Получим выражения в случае граничных условий третьего рода, которые являются наиболее универсальными, из них можно легко получить условия первого и второго рода, меняя параметры:

Эк , 0

— = ак + р. дп

Сначала рассмотрим узлы (7, у,к) е у, направим вектор нормали И параллельно оси Ох, (г +1, j, к) е О, где (г — 1, к) находится вне расчетной области, тогда

dh ön

dh öx

= ah + ß.

Тогда в узле (0, j, к) е/ выражение (к (к) к') ' примет вид:

1

1

Г

(кН\)Х = 1Г Х'+2(ШХ )'xdx =1 к К -Ч 1 К

кИ I - кН

1

2 У

1 К

\

Шх\ + к (аИ + В)

=% - ^)+к^К+в).

Теперь рассмотрим случай, когда (г +1, ¡, к) находится вне расчетной сет-

ки, тогда

дИ_ дп

г

ди

дх

= аН + В.

г

В этом случае выражение ( к ( Н) Н' ) ' примет вид

( \

кН'\ - кН'\

(кИ х )Х = 1 |;+2(кН'х=1

1 К

К Jx. 1

( \

к (аН + В) - кИ \

'х. 1

г-2

= к ,, 1(аНп+1, +В) - к , 1(Нп+\ - Нп+\1Л.

""¡'к Н 1,1,к 1 -■ >- 7,2 4 п,7,к п-1,1,к'

1 1 7,2^ п,7,к п-1,7,к* Ч-^'Зк И 1 ,¡,

Пусть теперь п \ \ Оу, (/, / +1, к) принадлежит расчетной сетке, а (г, 1 -1,к) - нет, тогда аналогично получим выражение для (к(Н)Н') ':

( \

(кИу )у = 1 ¡^(кН'у Уу*у=}

К

кИ

у

у+1

'+2

кН

- кН

у.., 1 2

у. 1

2

\

- к (аН + В)

=к 1 * - и& )+(ак+к+в)

V Ну Ну

Теперь (г, 1 -1, к) принадлежит расчетной сетке, а (г, 1 +1, к) - нет, тогда

получим выражение

(к(Н)Ну)у' '

(кН'у )у = ¥ ¡у+^'у )'у^=Н

которое примет вид:

у ,

( \ кН\\ -Шу\

. 2 2 У

И,

к (аИ + В) - кК

'ЧУ

= кЛк Н^аНщк + В) -к,п2-1,к ¿(ИЩк -Н!;;;^)

2

Теперь рассмотрим границы по направлению О/. Пусть и 11 Ог, (/, /, к +1) принадлежит расчетной сетке, а (г,¡,к -1) находится вне расчетной области, тогда для (к (И ) И' ) ' получим

2

1

2

1

(kh\)z - 1 \4{kh\)'zdz -±

w

- kh

z - 2,

1 К

kh

+ k (ah + ß)

Z 1 '+2

-kj j \ h2(hn+i - J+kj 0 hh+ß)-

Если же, наоборот, (/, к — 1) принадлежит расчетной сетке, а к +1) находится вне расчетной области, тогда для выражения ^к(к)к') ' получим

( \

(khz)Z - 1 \4(khz)zdz =±

kh'\ - kh\

z. i

'+2

z' - 2 у

1 К

k {ah + ß) - kh

Z'-2 у

= k 1(ahn+1

ß) - k 1тЖП - hj-i).

ujn-- h J n

Описание параллельного алгоритма. Для повышения оперативности получения результата моделирования целесообразно распараллеливать вычисления на сетках с большим количеством ячеек.

Параллельная версия программного комплекса разработана на языке высокого уровня С++ с использованием библиотеки передачи сообщений MPI.

В качестве метода решения СЛАУ был взят модифицированный попеременно-треугольный метод (МПТМ) [6]. Идея параллельного алгоритма для данного метода заключается в следующем: расчетная область разбивается на части по двум координатным направлениям, как показано на рис. 2; каждый процессор получает свою расчетную область по направлению, перпендикулярному плоскости разбиения, смежные области перекрываются двумя слоями узлов.

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

Рис. 2. Декомпозиция расчетной области

Рис. 3. Передача элементов после расчета двух слоев первым процессором

Эффективность параллельного алгоритма. Проведенные численные эксперименты показали, что максимальное ускорение для задачи размерности 351^251x14 достигалось на 128 процессорах и было равно 43.6. На рис. 4 показаны графики зависимости ускорения от количества процессоров, которые были рассчитаны теоретически и экспериментально. Сплошная кривая - теоретическая зависимость, а ломаная - экспериментальная.

60-:-1-1-1-1-

О 10 40 60 80 100 ¡20 140

Рис. 4. Зависимость ускорения от числа процессоров

В теоретических оценках ускорения рассматривается случай модельной задачи с прямоугольной областью. При решении задачи для реального объекта расчетная область может иметь довольно сложную форму, при этом реальное ускорение меньше его теоретической оценки.

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

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

БИБЛИОГРАФИЧЕСКИЙ СПИСОК 1. Кащенко Н.М., Никитин М.А. Моделирование эффектов аномальной диффузии для дренажных схем // Вестник Балтийского федерального университета им. И. Канта. - 2010. - № 4.

2. Шхануков-Лафишев М.Х., Такенова Ф.И. Разностные методы решения краевых задач для дифференциальных уравнений дробного порядка // Журнал вычислительной математики и математической физики. - 2006. - Т. 46, № 10. - С. 1871-1881.

3. Самарский А.А. Теория разностных схем. - М.: Наука, 1977. - 656 с.

4. Сухинов А.И., Чекина М.Д. Математическая модель и численный метод для задачи динамики выпадения осадков и затопления // Известия ЮФУ. Технические науки. - 2009.

- № 8 (97). - С. 42-52.

5. Сухинов А.И., Чекина М.Д. Математическое моделирование процессов накопления и фильтрации осадков с помощью супервычислительных систем // Известия ЮФУ. Технические науки. - 2010. - № 6 (107). - С. 103-113.

6. Сухинов А.И., Чистяков А.Е. Адаптивный модифицированный попеременно-треугольный итерационный метод для решения сеточных уравнений с несамосопряженным оператором // Математическое моделирование. - 2012. - Т. 24, № 1. - С. 3-20.

7. Коновалов А.Н. К теории попеременно-треугольного итерационного метода // Сибирский математический журнал. - 2002. - Т. 43, № 3. - С. 552.

8. Сухинов А.И., Чистяков А.Е., Тимофеева Е.Ф., Шишеня А.В. Математическая модель расчета прибрежных волновых процессов // Математическое моделирование. - 2012.

- Т. 24, № 8. - С. 32-44.

9. Сухинов А.И., Чистяков А.Е., Шишеня А.В. Оценка погрешности решения уравнения диффузии на основе схем с весами // Математическое моделирование. - 2013. - Т. 25, № 11. - С. 53-64.

REFERENCES

1. Kashchenko N.M., Nikitin M.A. Modelirovanie effektov anomal'noy diffuzii dlya drenazhnykh skhem [Simulation of the effects of anomalous diffusion for drainage schemes], Vestnik Baltiyskogo federal'nogo universiteta im. I. Kanta [Bulletin of Baltic Federal University I. Kant], 2010, No. 4.

2. Shkhanukov-Lafishev M.Kh., Takenova F.I. Raznostnye metody resheniya kraevykh zadach dlya differentsial'nykh uravneniy drobnogo poryadka [Difference methods for solving boundary value problems for differential equations of fractional order], Zhurnal vychislitel'noy matematiki i matematicheskoy fiziki [Journal of computational mathematics and mathematical physics], 2006, Vol. 46, No. 10, pp. 1871-1881.

3. Samarskiy A.A. Teoriya raznostnykh skhem [The theory of difference schemes]. Moscow: Nauka, 1977, 656 p.

4. Sukhinov A.I., Chekina M.D. Matematicheskaya model' i chislennyy metod dlya zadachi dinamiki vypadeniya osadkov i zatopleniya [Mathematical model and numerical method for the problem of the dynamics of rainfall and flooding], Izvestiya YuFU. Tekhnicheskie nauki [Izvestiya SFedU. Engineering Sciences], 2009, No 8 (97), pp. 42-52.

5. Sukhinov A.I., Chekina M.D. Matematicheskoe modelirovanie protsessov nakopleniya i fil'tratsii osadkov s pomoshch'yu supervychislitel'nykh sistem [Mathematical modelling of accumulation and filtration deposits processes by means of supercomputing systems], Izvestiya YuFU. Tekhnicheskie nauki [Izvestiya SFedU. Engineering Sciences], 2010, No. 6 (107), pp. 103-113.

6. Sukhinov A.I., Chistyakov A.E. Adaptivnyy modifitsirovannyy poperemenno-treugol'nyy iteratsionnyy metod dlya resheniya setochnykh uravneniy s nesamosopryazhennym operatorom [Adaptive modified alternating triangular iterative method for solving grid equations with non-self-adjoint operator], Matematicheskoe modelirovanie [Mathematical modeling], 2012, Vol. 24, No. 1, pp. 3-20.

7. Konovalov A.N. K teorii poperemenno-treugol'nogo iteratsionnogo metoda [The theory of alternating triangular iterative method], Sibirskiy matematicheskiy zhurnal [Siberian mathematical journal], 2002, Vol. 43, No. 3, pp. 552.

8. Sukhinov A.I., Chistyakov A.E., Timofeeva E.F., Shishenya A.V. Matematicheskaya model' rascheta pribrezhnykh volnovykh protsessov [The mathematical model of calculation of coastal wave processes], Matematicheskoe modelirovanie [Mathematical modeling], 2012, Vol. 24, No. 8, pp. 32-44.

9. Sukhinov A.I., Chistyakov A.E., Shishenya A.V. Otsenka pogreshnosti resheniya uravneniya diffuzii na osnove skhem s vesami [Error estimation for the solution of the diffusion equation-based schemes with weights], Matematicheskoe modelirovanie [Mathematical modeling], 2013, Vol. 25, No. 11, pp. 53-64.

Статью рекомендовал к опубликованию д.т.н., профессор Я.Е. Ромм.

Чекина Мария Дмитриевна - Научно-образовательный центр «Комплексные исследования и математическое моделирование» НИИ МВС ЮФУ; e-mail: elfik55@jmail.com; 347928, г. Таганрог, ул. Чехова, 2; тел.: 88634315491; программист.

Chekina Maria Dmitriyevna - Scientific Research Centre "Complex research and mathematic modelling" at SRI MCS SFU; e-mail: elfik55@jmail.com; 2, Chekhov street, Taganrog, 347928, Russia; phone: +78634315491; software developer.