NUMERICAL CALCULATION OF NONSTATIONARY FRACTIONAL DIFFERENTIAL EQUATION IN PROBLEMS OF MODELING TOXIC SUBSTANCES DISTRIBUTION IN GROUND WATERS
Abstract and keywords
Abstract (English):
At present, the theory of fractional calculus is widely used in many fields of science for modeling various processes. Differential equations with fractional derivatives are used to model the migration of pollutants in porous inhomogeneous media and allow a more correct description of the behavior of pollutants at large distances from the source. The analytical solution of differential equations with fractional order derivatives is often very complicated or even impossible. There has been proposed a numerical method for solving fractional differential equations in partial derivatives with respect to time to describe the migration of pollutants in groundwater. An implicit difference scheme is developed for the numerical solution of a non-stationary fractional differential equation, which is an analogue of the well-known implicit Crank-Nicholson difference scheme. The system of difference equations is presented in matrix form. The solution of the problem is reduced to the multiple solution of a tridiagonal system of linear algebraic equations by the tridiagonal matrix algorithm. The results of evaluating the spread of pollutant in groundwater based on the numerical method for model examples are presented. The concentrations of the substance obtained on the basis of the analytical and numerical solutions of the unsteady one-dimensional fractional differential equation are compared. The results obtained using the proposed method and on the basis of the well-known analytical solution of the fractional differential equation are in fairly good agreement with each other. The relative error is on average 9%. In contrast to the well-known analytical solution, the developed numerical method can be used to model the spread of pollutants in groundwater, taking into account their biodegradation.

Keywords:
fractional derivative, time derivative, fractional differential equation, difference scheme, groundwater
Text
Publication text (PDF): Read Download

В настоящее время многих исследователей интересует использование математических моделей с дробными производными в силу точности моделей при описании объектов в различных областях науки. Причинами возросшего интереса к дробным производным и дробным дифференциальным уравнениям являются следующие: 1) определение дробной производной нелокально и может быть использовано для математического моделирования сред с памятью; 2) дробные производные и дифференциальные уравнения дробных порядков точнее описывают рассматриваемые явления во многих задачах физики, биологии, геологии; 3) дробными дифференциальными уравнениями можно описывать немарковские процессы. Повышенный интерес к дробному исчислению, наблюдаемый в последнее время, позволяет говорить о нем как о новой области науки [1-10]. Особенный интерес к дробным производным проявляют гидрогеологи в связи с проблемой обоснования безопасности хранения опасных химических и радиоактивных отходов в геологических формациях [1, 2]. В соответствии с имеющимися экспериментальными данными изменение концентрации вещества на больших расстояниях от источника не может быть описано законами классической диффузии, поскольку подчиняется степенному закону убывания, а не экспоненциальному, как это считалось ранее. Данное явление связывают с тем, что область распространения подземных вод представляет собой неоднородную гетерогенную среду, которая включает различные по литологическому составу (по водопроницаемости) породы [1, 2]. Как указано в [1, 2], в таких случаях имеет место так называемая аномальная диффузия, которую подразделяют на субдиффузию и супердиффузию. Субдиффузия протекает намного медленнее классической диффузии, поскольку диффузант захватывается ловушками (дефектами, адсорбционно- или химически активными центрами), уходит в боковые, тупиковые пути и на некоторое время или навсегда выводится из миграционного процесса. Наоборот, супердиффузия осуществляется со скоростями, существенно превышающими классическую диффузию. Этот режим наблюдается, если в системе есть облегченные пути (например, трещины) или присутствуют процессы случайной или направленной адвекции (увлечение диффузанта потоками флюидов). Эти явления невозможно корректно описать классическим уравнением турбулентной диффузии. С учетом этого обстоятельства для более точного описания данных процессов необходимо использовать уравнения диффузии в частных производных дробных порядков [1, 2, 10]. Чтобы учесть в математической модели структуру среды или изменения внешнего поля со временем, нужно заменить обычный оператор производной оператором производной вариационного порядка. Дифференциальные уравнения в частных производных дробных порядков подразделяют на два вида: дифференциальные уравнения с дробной производной по пространству и дробной производной по времени. Чаще всего получить решение для таких уравнений в явном виде не удается. В настоящее время известны аналитические решения дробных дифференциальных уравнений для ограниченного круга задач. Например, в работе [11] приведено аналитическое решение дробного дифференциального уравнения, которое описывает распространение вещества в подземных водах без учета его возможных химических превращений. Поэтому для решения дробных дифференциальных уравнений в частных производных очень эффективными являются приближенные методы [3]. Несмотря на развитый математический аппарат дробного интегрирования и дифференцирования, вопросы численного решения уравнений с дробными производными практически не нашли отражения в существующей литературе [1, 2, 12-15]. В данной работе представлены результаты исследований по разработке численного метода и разностной схемы для решения одномерного дифференциального уравнения с дробной производной по времени. Наибольшее распространение в практике гидрогеологических расчетов получил метод конечных разностей. Однако данный метод обладает некоторыми недостатками [4, 5]. Избавиться от недостатков данного метода можно, используя различные подходы к построению разностной схемы [7]. В данной статье рассматривается уравнение параболического типа, описывающее процесс диффузии, и один из подходов к построению разностных схем для такого типа уравнений - схема Кранка - Николсона, называемая также схемой с весом 1/2. Устойчивость и сходимость решения конечно-разностных уравнений на основе схемы Кранка - Николсона доказана во многих работах [1, 2]. Применение данной схемы позволяет без существенных вычислительных затрат повысить (по сравнению с другими подходами) порядок аппроксимации разностной схемы [13-19]. Таким образом, наиболее целесообразной является разработка неявной разностной схемы для решения дробного дифференциального уравнения на основе известной схемы Кранка - Николсона. Неявная разностная схема для численного решения нестационарного одномерного дробного дифференциального уравнения Рассмотрим общее одномерное уравнение для моделирования миграции токсичных веществ в подземных водах с начальными и граничными условиями вида (1) начальные условия: c(x, 0) = c0, 0 < x < L; граничные условия: где c - концентрация токсичного вещества, мг/м3; c0- начальная концентрация токсичного вещества, мг/м3; D - коэффициент дисперсии в направлении оси х, м2/с; λ - константа скорости превращения токсичного вещества, 1/с; ν - скорость движения воды в порах почвы в направлении оси x, м/с; t - время, с; L - расстояние, на котором концентрация принимает значение с0, м; α - порядок производной. В 1967 г. М. Капуто ввел новое и полезное определение дробной производной, которое сегодня широко известно. Приведем это определение для дробной производной по времени порядка α [20, 21]: (2) где Г (.) - Гамма-функция. Разработаем неявную разностную схему, являющуюся аналогом схемы Кранка - Николсона, для численного решения дробного дифференциального уравнения. Схема Кранка - Николсона может быть записана в виде (3) В качестве сетки возьмем совокупность точек , координаты которых ; ; ; ; ; ; . На данной сетке будем строить разностную аппроксимацию дифференциального уравнения (1), используя шеститочечный шаблон: Ш . Следовательно, схема Кранка - Николсона для модели (1) и шеститочечного шаблона имеет вид (4) (5) (6) Приближенное решение задачи (1) получаем, заменяя значение исходной непрерывной функции на функцию, определяемую в узлах сетки (k - номер узла по времени и l - номер узла по пространству): и означает численную аппроксимацию точного решения (или обозначает численное приближение к точному решению). Аналог неявной разностной схемы Кранка - Николсона для численного решения дробного дифференциального уравнения (1) с учетом выражения для дробной производной Капуто по времени (2) можно сформулировать следующим образом: (7) Выполним последовательно преобразования правой части уравнения (7): (8) Запишем уравнение (1) с учетом (4)-(6), (8): (9) где ; . Проведем упрощающие преобразования. Примем, что коэффициент перед выражением в левой части равен единице, и введем следующие обозначения: (10) После преобразований (8)-(10) система разностных уравнений для первоначальной краевой задачи с дробной производной (1) имеет вид (11) Отметим, что при k = 0 слагаемое с суммой исчезает, что упрощает реализацию метода при проведении вычислений. Выполним преобразование уравнения (11): или (12) Коэффициенты в узлах сетки вычисляются по следующим формулам: Систему разностных уравнений (12) с учетом преобразований можно представить в следующем виде: или Первоначальная краевая задача с дробными производными (2) может быть представлена в матричном виде: для k = 0; (13) для , (14) где ; начальные условия и граничные условия: Матрицы А и В имеют следующий вид: ; 1 Решение задачи (13), (14) в конечном счете сводится к многократному решению системы линейных алгебраических уравнений (СЛАУ) с трехдиагональной матрицей. Для решения таких СЛАУ известны специальные методы, например метод прогонки и метод циклической редукции [9]. Наиболее простым и эффективным методом является метод прогонки, который привлекателен своей экономичностью и простотой реализации в последовательном режиме [8]. Именно поэтому метод прогонки использовался в данном случае для решения системы линейных алгебраических уравнений. Разработанный численный метод решения дробного дифференциального уравнения реализован в программном комплексе [22] для решения задачи распространения токсичных веществ в условиях аномальной диффузии. Алгоритм реализации численного метода представлен на рис. 1. Рис. 1. Алгоритм реализации разработанного численного метода Данный алгоритм предусматривает проверку корректности исходных и расчетных данных на всех этапах вычислений. Новый численный метод реализован на языке программирования Microsoft Visual Basic 6.5, в качестве системы управления базой данных использовался Microsoft Access. Оценка распространения токсичного вещества в подземных водах на основе численного метода В данном разделе приведены результаты численного эксперимента, иллюстрирующие вычислительные свойства нового метода. Выполнено сравнение концентраций вещества, полученных на основе решения нестационарного одномерного дробного дифференциального уравнения разработанным методом и на основании подхода, представленного в работе [11], для двух значений порядка производной по времени - α. Результаты приведены для модельного примера, в котором химические превращения токсичного вещества не рассматривались. Параметры расчетной сетки: h = 0,10; τ = 0,1. Число расчетных точек по осям - 100. На рис. 2 представлены графики зависимости концентрации вещества от расстояния, полученные с помощью разработанного численного метода и на основании аналитического решения, приведенного в работе [11], при α = 0,2. Рис. 2. Зависимость концентрации токсичного вещества от расстояния (α = 0,2) На рис. 3 представлены графики зависимости концентрации вещества от расстояния, полученные с помощью разработанного численного метода и на основании аналитического решения, приведенного в работе [11], при α = 0,6. Рис. 3. Зависимость концентрации токсичного вещества от расстояния (α = 0,6) Из рис. 2 и 3 видно, что графики зависимостей концентрации токсичного вещества от расстояния изменяются симбатно. Результаты, полученные с помощью предлагаемого метода и известного точного решения дробного дифференциального уравнения, достаточно хорошо согласуются между собой. Относительная ошибка составляет в среднем 9 %. В отличие от известного аналитического решения разработанный численный метод может использоваться при моделировании распространения токсичных веществ в подземных водах с учетом их биодеградации. Заключение В результате проведенных исследований разработан новый численный метод решения дифференциального уравнения с производной дробного порядка по времени, а также неявная разностная схема, являющаяся аналогом схемы Кранка - Николсона. Система разностных уравнений представлена в матричном виде. Решение задачи сводится к многократному решению трехдиагональной системы линейных алгебраических уравнений методом прогонки. Представлены результаты оценки распространения токсичных веществ в подземных водах на основе численного метода при решении модельных задач. Полученные результаты позволяют рекомендовать разработанный метод численного решения нестационарного одномерного дробного дифференциального уравнения, основанный на использовании линейных интегро-дифференциальных уравнений достаточно общей структуры, для использования на практике при моделировании процессов тепло- и массопереноса в сильно неоднородных средах.
References

1. Goloviznin V. M., Kiselev V. P., Korotkin I. A., Yurkov Yu. I. Nekotorye osobennosti vychislitel'nyh algoritmov dlya uravneniy drobnoy diffuzii. M.: Izd-vo In-ta problem bezopasnogo razvitiya atomnoy energetiki RAN, 2002. 57 s.

2. Goloviznin V. M., Kiselev V. P., Korotkin I. A. Chislennye metody resheniya uravneniy drobnoy diffuzii v odnomernom sluchae. M.: Izd-vo In-ta problem bezopasnogo razvitiya atomnoy energetiki RAN, 2002. 35 s.

3. Merkulova N. N., Mihaylov M. D. Metody priblizhennyh vychisleniy: ucheb. posobie. Tomsk: Izd-vo TGU, 2011. Ch. III. 184 s.

4. Samarskiy A. A., Gulin A. V. Chislennye metody. M.: Nauka, 1989. 432 s.

5. Verzhbickiy V. M. Chislennye metody. M.: Vyssh. shk., 2001. 382 s.

6. Baykov V. A., Zhiber A. V. Uravneniya matematicheskoy fiziki. M.; Izhevsk: Izd-vo In-ta komp'yuternyh issledovaniy, 2003. 252 s.

7. Samarskiy A. A., Nikolaev E. S. Metody resheniya setochnyh uravneniy. M.: Nauka, 1987. 130 s.

8. Tihonov A. N., Samarskiy A. A. Uravneniya matematicheskoy fiziki. M.: Nauka, 1977. 735 s.

9. Kustikova V. D. Differencial'nye uravneniya v chastnyh proizvodnyh. N. Novgorod: Izd-vo NGU im. N. I. Lobachevskogo, 2011. 54 s.

10. Hendi Ahmed Sand Abdelaziz. Chislennye algoritmy resheniya drobnyh differencial'nyh uravneniy s zapazdyvaniem: dis. … kand. fiz.-mat. nauk. Ekaterinburg, 2017. 124 s.

11. Shukla H. S., Tamsir M., Srivastava V. K., Kumar J. Approximate Analytical Solution of Time-fractional order Cauchy-Reaction Diffusion equation // Computer Modeling in Engineering & Sciences. 2014. V. 103. N. 1. R. 1-17.

12. Samko S. G., Kilbas A. A., Marichev O. I. Integraly i proizvodnye drobnogo poryadka i nekotorye ih prilozheniya. Minsk: Nauka i tehnika, 1987. 687 s.

13. Shtumpf C. A., Bahtin M. A. Matematicheskie metody komp'yuternyh tehnologiy v nauchnyh issledovaniyah: ucheb. posobie. SPb.: Izd-vo NIU ITMO, 2012. 148 s.

14. Fadugba S. E., Nwozo C. R. Crank-Nicolson Finite Difference Method for the Valuation of Options // The Pacific Journal of Science and Technology. 2013. V. 14. N. 2. R. 136-145.

15. Wani S. S., Thakar S. H. Crank-Nicolson Type Method for Burgers Equation // International Journal of Applied Physics and Mathematics. 2013. V. 3. N. 5. R. 324-328.

16. Kilbas A. A., Trujillo J. J. Differential equations of fractional order: methods, results and problems // Applicable Analysis. 2001. V. 78. N. 1. R. 153-192.

17. Atangana1 A., Kilicman A. On the Generalized Mass Transport Equation to the Concept of Variable Fractional Derivative // Mathematical Problems in Engineering. 2014. R. 1-9.

18. Erohin S. V. Matematicheskoe modelirovanie napryazhenno-deformirovannogo sostoyaniya vyazkouprugih tel s ispol'zovaniem metodov drobnogo ischisleniya: dis. … kand. tehn. nauk. M., 2016. 130 s.

19. Begli R. L., Torvik P. Dzh. Differencial'noe ischislenie, osnovannoe na proizvodnyh drobnogo poryadka - novyy podhod k raschetu konstrukciy s vyazkouprugim dempfirovaniem // Aerokosmicheskaya tehnika. 1984. T. 2. № 2. S. 84-91.

20. Caputo M. Linear model of dissipation whose Q is almost frequency independent // Geophis. J. R. Astr. Soc. 1967. V. 13. P. 529-539.

21. Caputo M., Mcinardi F. Linear models of dissipation in an elastic solids // Riv. Nuovo Cimento. 1971. V. 1. R. 161-196.

22. Svidetel'stvo o gosudarstvennoy registracii programmy dlya EVM № 2019614483. Programmnoe obespechenie dlya modelirovaniya rasprostraneniya zagryaznyayuschih veschestv v podzemnyh vodah / Afanas'eva A. A., Nazarenko D. I., Shvecova-Shilovskaya T. N., Kazarezova E. V.; zareg. 26.03.2019; opubl. 05.04.2019.


Login or Create
* Forgot password?