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

УДК 539.219.3 ББК В375.6

Анна Николаевна Корчагина

аспирант,

Институт гидродинамики им. М. А. Лаврентьева Сибирского отделения Российской академии наук (Новосибирск, Россия), e-mail: anchouse@ngs.ru

Лев Алексеевич Мержиевский, доктор физико-математических наук, профессор, Институт гидродинамики им. М. А. Лаврентьева Сибирского отделения Российской академии наук (Новосибирск, Россия), e-mail: merzh@hydro.nsc.ru

Численное моделирование диффузионных процессов в фрактальных средах1

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

Ключевые слова: аномальная диффузия, дробные производные, фрактальная среда.

Anna Nikolaevna Korchagina

Postgraduate Student, Lavrent ’ev Institute of Hydrodynamics, Siberian Branch, Russian Academy of Sciences (Novosibirsk, Russia), e-mail: anchouse@ngs.ru Lev Alekseevich Merzhievskiy Doctor of Physics and Mathematics, Professor, Lavrent ’yev Institute of Hydrodynamics, Siberian Branch, Russian Academy of Sciences, (Novosibirsk, Russia), e-mail: merzh@hydro.nsc.ru

Numerical Modeling of Diffusive Processes in Fractal Media

Derivatives of a fractional order are used for modeling of anomalous diffusion. Various definitions of fractional derivatives are considered, comparison of numerical solutions of a number of problems of diffusion by various numerical methods is carried out. The most perspective definitions and methods of the numerical decision are specified.

Keywords: anomalous diffusion, fractional derivatives, fractal media.

Введение. Классическое описание процессов диффузии базируется на законах Фика. Следствием из второго закона является классическое дифференциальное уравнение диффузии. В последние годы сформировался повышенный интерес к исследованию диффузионных процессов, не подчиняющихся законам Фика и не описывающихся классическим уравнением. Явления переноса, не укладывающиеся в классические представления, наблюдаются, например, в турбулентных потоках, в аморфных полупроводниках, высокоэнергетической плазме, пористых средах. Эти явления получили название «аномальная диффузия». Довольно полное представление о состоянии развития исследований аномальной диффузии применительно к различным задачам физики дано, например, в [1; 2]. Одним из проявлений «аномальности» является диффузия в гетерогенных, в частности во фрактальных, средах.

Для описания таких процессов используется модифицированный закон Фика [3], что требует привлечения математического аппарата дробного интегро-дифференциального исчисления [4]. В классическое уравнение диффузии вводятся производные дробного порядка как по пространству, так и по времени. Возникают начально-краевые задачи для дифференциальных уравнений с дробными производными. Развиваются аналитические методы решения задач, однако наибольшее распространение получили численные методы [5-10]. Это связано, в первую очередь, с тем, что аналитические решения удается получить только в редких частных случаях.

1Работа выполнялась при поддержке Интеграционного проекта СО РАН № 64 и гранта РФФИ № 12-01-00726-а.

© Корчагина А. Н., Мержиевский Л. А., 2013

53

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

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

Определения дробных производных. Существует ряд различных подходов к определению понятия производной дробного порядка, отражающих особенности становления дробного исчисления. Наиболее широким и часто используемым является определение Римана-Лиувилля, основанное на обобщении уравнения Абеля [4]:

X

D“u(x) = —------- —- [ ---[fp———г, х > a, 0 < то — 1 < a < то. (1)

x v ’ Г(т - a) dxm J (x - £)«-m+i’ - - w

a

Здесь использованы стандартные обозначения оператора дифференцирования и Г-функции.

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

X

D^uix) = —------- [(х — dt;, х > а, 0 < то — 1 < a < то. (2)

Г (m - a) J

a

Операция дробного дифференцирования Римана-Лиувилля обратна операции дробного интегрирования. Дробная производная в форме Капуто этим свойством не обладает. Развивая идею Лиувилля, А. Грюнвальд и - независимо - А. В. Летников ввели понятие дробной производной, как предела разностных отношений. Согласно определению Грюнвальда-Летникова, правая дробная производная определяется выражением

= lim д= Уурф - (к -ш (3)

dxa h^o hah^k

k=0

Биномиальные коэффициенты имеют вид

nfc Г(а + 1) _ Г {к-а)

к ( ’ \к) ( ’ Т(к + 1)Г(а — к + 1) Г(-а)Г(*+1)' ( '

Если u(x) непрерывна, а du/dx интегрируема на отрезке [a, x], то производные Римана-Лиувилля, Капуто и Грюнвальда-Летникова существуют и совпадают.

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

Уравнение диффузии с дробными производными. Для вывода уравнения диффузии с дробными производными используется соответствующий вариант модифицированного закона Фика [3], тогда оно принимает вид (здесь Ko = const):

D] u(x,t) = KoD^u(x,t),

1 da 1 da

^ = 2<1 + «a? + 2(1-«5FiF- <5>

0 <Y < 2,1 < a < 2, -1 < ß < 1.

Здесь a - дробный порядок дифференцирования по пространству; ß - «коэффициент скошенности», который характеризует направление переноса вещества при a ^ 1; y - дробный порядок

дифференцирования по времени. Дробная производная по пространству возникает в случае фрак-тальности среды, параметр дифференцирования зависит от хаусдорфовой размерности фрактала. Дробная производная по времени возникает при учёте нелокальности по времени, которая связана с прилипанием диффундирующих атомов к стенкам пор [7]. При a = 2 получаем уравнение классической диффузии. Случай 1 < a < 2 отвечает «быстрой» диффузии (super-diffusion), когда частицы распространяются быстрее, чем предсказывает классическая модель. И случай a =1 - это классический перенос. Для коэффициента a рассматриваются следующие области значений:

0 < y < 1 — «медленная» диффузия (slow diffusion, sub-diffusion);

1 < y < 2 — «быстрая» диффузия (fast diffusion, hyper-diffusion);

Y =1 — обычная, классическая диффузия .

В режиме субдиффузии скорость роста среднеквадратичного смещения частиц монотонно убывает со временем, тогда как в режиме супердиффузии скорость со временем возрастает [8]. При Y ^ 1 рассматриваемое уравнение переходит в классическое уравнение диффузии с экспоненциальным затуханием решения на бесконечности. При y ^ 2 получаем волновое уравнение. Как варианты, могут рассматриваться уравнения, в которых только одна из производных заменяется на дробную.

Методы численного решения. Методы численной аппроксимации дробных производных напрямую связаны с их определениями. Детальный анализ применимости разных методов аппроксимации и существующих разностных схем решения уравнений с дробными производными осуществлен в [6; 9; 10]. Поясним основные идеи использованных методов на некоторых примерах. Для упрощения анализа результатов по разным методам рассмотрим случай, когда в уравнения вводится дробная производная только по времени. Представим уравнение теплопроводности в виде:

did7 1u(x,t) dt 1

д 2u(x, t)

(6)

u?+11-2u?+1 + u£11

д^-1 ) дх2 ’

1 < 7 < 2. Этому уравнению сопоставим следующий разностный аналог:

7-1Ьм”+1 -7-1 ,

т Ъ?

где 1-1Ь - численная аппроксимация оператора дробной производной порядка 7 — 1. Воспользуемся определением производной дробного порядка в смысле Римана-Лиувилля. Представим оператор 7-1Ь в виде конечной суммы интегралов по отрезкам < £ < ¿*+1, расположенным между узлами расчётной сетки:

(7)

7-1

Lu” =

1

d

г(2 — Y) dtn

^k + l

Е

k=0

i(£) d£

(tn — e)Y-1

+

i(£) d£

(tn — e)^-1

tk = кт

(8)

Функция м(£) на отрезках < £ < ¿*+1 аппроксимируется линейно:

«(£) = А е + в* . (9)

Тогда, вычисляя аналитически интегралы в скобках, получаем выражение для оператора 7-1Ь:

Y-1L

^0\ 0 0 0 0 /u0^

u1 0 0 0 0 01 u1

u2 0 0 2 0 u2

u3 0 со V to v V 0 u3

Ку \0n Vn-1 Vn-2 . . . V0y \u”)

(10)

в0 ^к1^7 - Ь2 7 ^ ^ , к = 1, 2,..., Т.

1

9

0

0

k

1 Хк = Ао [(к + 1)2~7 - 2 к2-~< + (к- 1)2~7] , к = 1,2,...,Т.

0 т7-1(2 — 7 )Г(2 — 7)

Выделим из левой части уравнения слагаемое, относящееся к слою п + 1:

и”+1 . 7-1Ьи”+1 -7-1 Ьи” <+1 - 2и”+1 + м”—-,1

і ____г_______г_ _____________ ^і~і_£_ г— і М 1 1

гТ-1(2-7)Г(2-7) г " К2 ’ ^ ;

где 7-1Ь - то же самое, что и 7-1Ь, но с коэффициентом Ао .Теперь уравнение можно представить в виде:

Аи^+1 - Сип+1 + Ви^+1 + В = 0,А = В =-1,

2 7-1Ьм"-7-1/^м"+1

С’=4 + Ло,Д =--^(12)

д2 т

и решить методом трехточечной прогонки:

иі = иі+1аі+1 + ^¿+Ъ а*+1 = 7^ д і ^¿+1 = 7^ д • (13)

С — Ва* С — Ва*

Коэффициенты а1, 61 и находятся из краевых условий. Порядок аппроксимации данной схемы: 0(т + Д2).

Для производной Капуто оператор 7-1Ь представим в следующем виде:

п—1 Ік + 1

Е / (<" -о1'7^(»)

Как и раньше, аппроксимируем искомую функцию линейной на каждом из отрезков < і < і&+1 и вычислим соответствующие интегралы. После преобразований, оператор приводится к тому же виду, что и в случае производной Римана-Лиувилля с коэффициентами 9^:

0Сар«4О = А0 ((к — 1)2-7 — к2-7) , к = 1, 2, . . . , Т. (15)

Для построения численной схемы в случае производной Грюнвальда-Летникова возьмем от

обеих частей уравнения теплопроводности дробную производную порядка 2 — 7:

&и2^сРи

дЬ2 4 дх2' ( ;

Полученное уравнение аппроксимируем конечно-разностным с использованием определения производной Грюнвальда-Летникова:

.”+1 О*,” I ^,П— 1 1 ____

\ л ^,2—7 (ап—\ _________ Ог11п~1 г11п~1

^=0

Е"*"’ ("« - 2»г‘ + “Г-11) ■ <17>

порядок аппроксимации которой также 0(т + Л2). Схема является условно устойчивой, достаточное условие устойчивости: рг < ^тЦг [8]. Кроме перечисленных авторами рассматривались и другие варианты определений дробных производных и их численные аппроксимации.

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

Результаты решения задач. Приведём результаты решения по описанным методикам двух модельных задач. Решение в безразмерных величинах проводилось на отрезке 0 < х < 2п.

Задача 1. Диффузия с правой границей области. Начально-краевые условия:

дм

и(х, 0) = 0.02; — v ’ 7 ’ ді

= 0; м(0,і) = 0. 02; м(2п,і) = 2. (18)

(х,0)

3 4

Рис. 1

і 4

Рис. 2

3 4

Рис. 3

Рис. 4

-0.8 Рис. 5

Рис. 6.

Результаты решения задачи для двух значений параметра y приведены на рис. 1, где y = 1 (классическое уравнение), и рис. 2, где y = 0.6. Влияние на решение величины y прослеживается на рис. 3 и 4, где приведены рассчитанные значения функции u(x) для различных значений y в фиксированные моменты времени.

Задача 2. Эволюция примеси из локальной области. Начально-краевые условия:

du

u(x, 0) = S(x — 7г); — = 0;«(0,i) = u(27r,i) = 0. (19)

Результаты решения для линейного уравнения диффузии с дробной производной по времени при y = 0.8; 1; 1.2 на фиксированный момент времени показаны на рис. 5. Решение этой же задачи в случае дробной производной по пространству при a = 1,15 на разные моменты времени приведены на рис. 6.

Обсуждение результатов, выводы. На основе полученных решений уравнений с дробны-

ми производными по времени либо пространству проанализирована зависимость поведения решений от параметров порядка дифференцирования и кососимметричности. При 7 < 1 скорость протекания процесса вначале больше скорости классической диффузии, но с течением времени наблюдается замедление, характерное для субдиффузии. При 7 > 1 скорость процесса выше, чем в классическом случае, и процесс с течением времени ускоряется. В этом случае проявляются «волновые» свойства решения.

Решения уравнений с дробной производной по пространству показывают, что зависимость скорости диффузии от порядка дробной производной, оказывающейся большей, чем предсказывает классическая модель. При приближении параметра дифференцирования к 1 наблюдается явно выраженный процесс переноса (рис. 6).

Сравнение результатов, полученных при использовании разных определений дробных производных (Римана-Лиувилля, Капуто, Грюнвальда-Летникова) и соответствующих разностных аппроксимаций, показало, что получаемые для рассмотренных задач данные практически совпадают. Это означает, что для решения данного класса конкретных краевых задач эти методы равноценны и дают решения, достаточно близкие к полученным в некоторых случаях аналитическим.

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

Список литературы

1. Metzler R., Klafter J. The random walk’s guide to anomalous diffusion: a fractional dynamics approach // Phys. Rep. 2000. V. 339 P. 1.-77.

2. Учайкин В. В. Автомодельная аномальная диффузия и устойчивые законы //

УФН 2003. Т. 173, № 8. С. 847-876.

3. Paradisi P., Cesari R., Mainardi F., Tampieri F. The fractional Fick’s law for non-local transport processes // Physica A. 2001. Т. 293 P. 130-142.

4. Самко С. Г., Килбас А. А., Маричев О. И. Интегралы и производные дробного порядка и некоторые их приложения. Минск: Наука и техника, 1987.

5. Gorenflo R. Fractional calculus: some numerical methods // Fractals and Fractional Calculus in Continuum Mechanics, eds. A. Carpinteri and F. Mainardi. Springer Verlag, Wien,

1997. №. 378. P. 277-290.

6. Мержиевский Л. А., Корчагина А. Н. Сравнение методов численного решения задач для уравнения теплопроводности дробного порядка // X Международный семинар «Супервычисления и математическое моделирование». Саров, 2008. С. 85-86.

7. Головизнин В. М., Киселёв В. П., Короткин И. А. Численные методы решения уравнения дробной диффузии в одномерном случае. М., 2002 (Препринт / ИБРАЭ РАН: IBRAE-2002-01).

8. Лукащук С. Ю., Костригин И. В. Численное решение диффузионно-волновых уравнений дробного порядка на кластерных системах // Труды VI Всероссийской конференции молодых ученых по мат. моделированию и информ. технологиям. Кемерово,

2005. С. 19.

9. Мержиевский Л. А., Корчагина А. Н. Моделирование распространения теплового импульса во фрактальной среде // Экстремальные состояния вещества. Детонация. Ударные волны. Труды международной конференции «XI Харитоновские тематические научные чтения». Саров, 2009. С. 250-254.

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

References

1. Metzler R., Klafter J. The random walk’s guide to anomalous diffusion: a fractional dynamics approach // Phys. Rep. 2000. V. 339 P. 1.-77.

2. Uchayykin V. V. Avtomodelnaya anomaliya diffuziya i ustoychivye zakony // UFN 2003. T. 173. № 8. S. 847-876.

3. Paradisi P., Cesari R., Mainardi F., Tampieri F. The fractional Fick’s law for non-local transport processes // Physica A. 2001. Т. 293 P. 130-142.

4. Samko S. G., Kilbas A. A., Marichev O. I. Integraly i proizvodnye drobnogo poryadka

i nekotorye ikh prilozheniya. Minsk: Nauka i tekhnika, 1987.

5. Gorenflo R. Fractional calculus: some numerical methods // Fractals and Fractional Calculus in Continuum Mechanics, eds. A. Carpinteri and F. Mainardi. Springer Verlag, Wien,

1997. №. 378. P. 277-290.

6. Merzhiyevsky L. A., Korchagina A. N. Sravneniye metodov chislennogo resheniya zadach dlya uravneniya teploprovodnosti drobnogo poryadka // X Mezhdunarodny seminar «Supervychisleniya i matematicheskoye modelirovaniye». Saratov, 2008. S. 85-86.

7. Golovizin V. M., Kiselyov V. P., Korotkin I. A. Chislennye metody uravneniya drobnoy diffuzii v odnomernom sluchaye. M., 2002 (pereprint /IBR AE RAN. IBRAE-2002-01)

8. Lukashchuk S. Yu. Kostrigin I. V. Chislennoye resheniye diffuzno-volnovykh uravneny drobnogo poryadka na klasternykh sistemakh // Trudy VI Vserossyskaya konferentsiya molodykh uchenykh po mat. modelirovaniyu i inform. tekhnologiyam. Kemerovo, 2005. S.

19.

9. Merzhiyevsky L. A., Korchagina A. N. Modelirovaniye raspredeleniya teplovogo impulsa vo fraktalnoy srede // Ekstremalnye sostoyaniya veshchestva. Detonatsiya. Udarnye volny. Trudy mezhdunarodnoy konferentsii «XI Kharitonovskiye tematicheskiye nauchnye chteniya». Sarov, 2009. S. 250-254.

10. Taukenova F. I., Shkhanukov-Lafishev M. Kh. Raznostnye metody resheniya krayevykh zadach dlya differentsialnykh uravneny drobnogo poryadka / / Zhurnal vychislitelnoy matematiki i matematicheskoy fiziki. 2006. T. 46. № 10. S. 1871-1881.

Статья поступила в редакцию 25.04-2013