Содержание к диссертации
Введение
Глава 1. Элементы дробного интегро-дифференциального исчисления. Дробная диффузия 40
1.1. Понятие дробной производной 40
Определение Вейля 41
Определение Грюнвальда-Летникова 41
Определение Римана-Лиувилля 42
Примеры производных дробного порядка 43
Дробный Лапласиан 44
1.2. Модель дробной диффузии 48
1.3. Вывод уравнения дробной диффузии 52
Способ Эйлера 52
Распределение Леви 57
Глава 2. Методы численного решения уравнения диффузии с дробной производной по пространству 60
2.1. Постановка задачи 60
2.2. Фундаментальные решения 61
2.3. Метод Фурье 62
2.4. Конечно-разностный метод 65
2.5. Метод сплайн-аппроксимаций 68
2.6. Анализ численных методов 74
Устойчивость явной разностной схемы 74
Стационарная задача. Аппроксимация 78
2.7. Двумерные уравнения дробной диффузии и метод быстрого преобразования Фурье для их численного решения 85
Глава 3. Методы численного решения уравнения диффузии с дробной производной по времени 88
3.1. Квазиволновое представление уравнения диффузии с дробной производной по времени 89
3.2. Фундаментальные решения 91
Задача Коши (Cauchy Problem) 94
Краевая задача (Signalling Problem) 100
3.3. Постановка вычислительной задачи 101
3.4. Методы, основанные на определении Римана-Лиувилля 103
Начальное условие на первом временном слое 104
Вычислительные особенности 106
Квазиволновое представление 110
3.5. Методы, основанные на определении Капуто 115
Начальное условие на первом временном слое 115
Вычислительные особенности 116
Квазиволновое представление 119
Глава 4. Численное моделирование стохастической адвекции радионуклидов в сильно неоднородных средах 126
4.1. Генерация случайных соленоидальных полей скорости с заданной корреляцией в случае двух пространственных измерений 128
Результаты работы алгоритма 130
4.2. Алгоритм численного решения двумерного уравнения конвективного переноса 134
Постановка и дискретизация задачи 134
Вычисление новых значений консервативных переменных 135
Коррекция консервативных переменных 136
Вычисление новых значений потоковых переменных 137
Тестирование алгоритма 141
4.3. Задача о стохастической адвекции примеси 145
Случайные поля скоростей 147
Распространение примеси с течением времени 153
«Тяжелые хвосты» 163
Заключение 170
Список литературы
- Определение Грюнвальда-Летникова
- Конечно-разностный метод
- Постановка вычислительной задачи
- Алгоритм численного решения двумерного уравнения конвективного переноса
Введение к работе
Состояние вопроса. Проблема обращения с отработанным ядерным топливом и высокоактивными отходами приобрела актуальность еще к середине 60-х годов прошлого века и является одной из наиболее серьезных проблем во всех странах, использующих ядерные энергетические технологии. В ближайшие десятилетия наиболее реальным технически достижимым методом, способным обеспечить изоляцию отработанного ядерного топлива и остеклованных высокоактивных отходов от окружающей среды в течение 10 тыс. или более лет, является удаление этих ядерных отходов в глубокозалегающие породы — гранит, глина, соль, туф.
Решение о переработке или окончательном захоронении отработанного ядерного топлива зависит в основном от затрат на ядерный топливный цикл. Результаты экономического анализа замкнутого (с переработкой отработавшего топлива с последующим возвращением в цикл) и открытого (прямое удаление отработавшего топлива без переработки) ядерных топливных циклов в Европе показали некоторое превышение затрат в замкнутом цикле по сравнению с открытым циклом. Тем не менее, ряд стран (Россия, Великобритания, Франция, Япония), учитывая экологическое преимущество (снижение объема и радиоактивности отходов при захоронении после переработки) и возможность рециклирования извлеченных плутония и урана, приняли стратегию замкнутого цикла. США, Канада, Швеция и другие страны предпочли открытый цикл и имеют программы строительства хранилищ для удаления отработанного ядерного топлива в геологических формациях. Начало эксплуатации геологических хранилищ ожидается в США (Юкка-Маунтин) не ранее 2010 г., Швеции — в 2015 г., Финляндии (Олкилуото) — в 2020 г.
Накопление больших объемов радиоактивных отходов и отработанного ядерного топлива с тепловых реакторов поставило вопрос о выработке и реализации долгосрочной стратегии обращения с отходами и в России. Ближайшее время предполагается определить места сооружения крупных могильников для централизованного хранения долгоживущих высокоактивных отходов на территории РФ. Их создание позволило бы решить на самом современном уровне проблему обращения с отходами ядерной энергетики в стране.
Несмотря на интенсивные исследования и большие капиталовложения в разработку и строительство хранилищ, проблема высокоактивных отходов считается сейчас одной из сгещр слоящых, в атомлой
БИБЛИОТЕКА С.-Пет«р9ург
ОЭ ЖйнтЦ?^
энергетике. При анализе возможности строительства долговременных хранилищ радиоактивных отходов в геологических формациях принципиальным является вопрос экологической безопасности таких захоронений. Необходимы надежные оценки скорости миграции ра-диоігуклидов в сильно неоднородных средах с полным и неполным влагонасыщением. Многие из нуклидов, которые содержатся в продуктах переработки ядерного топлива, имеют очень большой период полураспада, исчисляемый сотнями и тысячами лет, а для распада трансурановых нуклидов (актинидов) нужно время порядка миллиона лет. Понятно, что сегодня у людей нет опыта контроля над какими-либо техническими устройствами в течение столь длинного периода времени. Решить данную проблему можно с помощью надежного прогнозирования поведения радионуклидов в подземном хранилище. Один из подходов к решению данной задачи заключается в построении математических моделей, адекватно описывающих перенос радионуклидов в сильно неоднородных, неупорядоченных геологических структурах.
Отличительной особенностью классических моделей диффузионного переноса является то, что согласно им размер облака частиц
примеси на больших временах растет как а убывание кон-
центрации в области «хвостов» (то есть на больших расстояниях, при г » R(t) ) происходит по гауссову закону. Вместе с тем, результаты
лабораторных и полевых экспериментов показывают, что для сильно неоднородных неупорядоченных сред, какими являются трещиноватые горные породы, неклассический характер переноса представляется скорее правилом, чем исключением. При этом рост облака частиц примеси со временем ( #(0) может идти быстрее, чем R~4t ,& убывание концентрации при г » R(t) — медленнее, чем по гауссову закону. Если «хвосты» убывают по степенному закону, то их называют «тяжелыми». Отсюда следует, что предельно допустимые концентрации радионуклидов могут распространяться на расстояния, которые на много порядков превышают найденные из классических представлений. Все это означает, что для получения адекватных оценок надежности радиоактивных захоронений, требуется разработка новых подходов.
Простейшими математическими моделями переноса в сильно неоднородных средах, характеризующихся отсутствием пространственного масштаба, являются модели случайного блуждания частиц, в
которых среда представляется однородной, а ее стохастические свойства проявляются в выборе функции распределения приращений координат блуждающей частицы. Если приращения координат происходят через одинаковые промежутки времени, и функция приращений имеет конечную дисперсию, то изменение со временем плотности пространственного распределения частиц описывается классическим уравнением диффузии.
Если дисперсия приращений координат становится бесконечной, то в одномерном случае средняя плотность описывается уравнением дробной диффузии с оператором дробного пространственного дифференцирования, определяемым двумя параметрами, характеризующими функцию распределения. Это параметр а, задающий степенное убывание «хвостов», и параметр f$, задающий степень асимметрии распределения. При а = 2 распределение стремится к нормальному, при а = 1 — переходит в распределение Коши. При а < 1 функция распределения приращений координат блуждающей частицы характеризуется не только бесконечной дисперсией, но и бесконечно большой средней величиной. И хотя, в этом случае, мы не можем получить устойчивых средних для каждой из частиц в отдельности, средние значения плотности частиц в заданной точке пространства остаются устойчивыми переменными и имеют привычный физический смысл.
Если перемещения частиц осуществляются не через каждый фиксированный промежуток времени, а иногда происходят задержки в перемещениях, то можно ввести среднее время ожидания частицей акта перемещения. Если это среднее время конечно, то уравнение дробной диффузии будет содержать производную по времени первого порядка. При бесконечном среднем времени ожидания, производная по времени в уравнении дробной диффузии также становится дробной, с показателем, меньшим единицы.
Описанный подход к моделированию переноса консервативной примеси в сильно неоднородных средах по своей сути является феноменологическим. Основным при таком подходе является вопрос об определении показателей дробных производных, или параметров функции распределения приращений координат частиц в их случайном блуждании, по результатам лабораторных измерений или натурных наблюдений. Такую задачу естественно назвать обратной. Для отработки методик решения обратной задачи необходимо использовать расчеты прямой задачи о временной эволюции заданного на-
чального распределения переносимой субстанции. Для этого целесообразно использовать алгоритмы численного решения интегро-дифференциальных уравнений дробной диффузии, либо алгоритмы прямого стохастического моделирования.
Другим возможным подходом к моделированию распространения примеси в зоне полного влагонасыщения является применение модели стохастической адвекции. Предложенная модель трактует перенос радионуклидов как адвекцию на случайном поле скоростей с дально-действующими корреляциями скорости.
Модель стохастической адвекции представляет собой яркий пример того, как во фрактальной, сильно неупорядоченной среде может реализоваться режим дробной диффузии (супердиффузии) на умеренных расстояниях с одновременным отсутствием «тяжелых» степенных «хвостов» в распределении концентрации на больших расстояниях.
Цель работы. Целью работы является разработка вычислительных методов решения уравнений диффузии с дробной производной по пространству и по времени в одномерном случае, разработка алгоритма для проведения прямого численного моделирования стохастической адвекции примеси в сильно неоднородных геологических средах, создание и отладка соответствующих программ. Результаты проделанной работы предполагается использовать для оценки безопасности подземных захоронений и хранилищ высокоактивных ядерных отходов.
Научная новизна. Общепринятый подход при анализе безопасности захоронений долгоживущих радионуклидов в геологических формациях заключается в использовании математической модели фильтрации в трещиновато-пористой среде, сводящейся к стандартному уравнению конвекции-диффузии относительно концентрации переносимой субстанции. Однако имеются экспериментальные данные, свидетельствующие о том, что в сильно неупорядоченных средах классические закономерности переноса примеси могут нарушаться. В данной работе предложены новые методики, претендующие на адекватное описание процессов переноса радионуклидов в сильно неоднородных средах. Это модель так называемой дробной диффузии и модель стохастической адвекции с дальнодействующи-ми корреляциями скорости.
Для модели дробной диффузии разработаны новые вычислительные алгоритмы решения уравнений диффузии с дробной производной по пространству и с дробной производной по времени в одномерном случае. Предложены способы обобщения алгоритмов на многомерный случай.
Для модели стохастической адвекции разработан численный алгоритм генерации случайных соленоидальных двумерных полей скоростей и создан новый численный метод решения задачи конвективного переноса примеси в прямоугольной области. Данный метод обладает минимальной сеточной диффузией, улучшенными дисперсионными свойствами и является модернизацией и обобщением одномерного алгоритма «КАБАРЕ». В ходе совершенствования двумерной схемы «КАБАРЕ» разработаны и апробированы на тестовых задачах новые алгоритмы коррекции консервативных и потоковых величин.
Практическая ценность работы. В рамках описанного в работе подхода при оценке безопасности подземных хранилищ радиоактивных отходов наиболее важной является задача нахождения неизвестных параметров дробной диффузии, характеризующих особенности неклассического переноса в данной геологической формации, то есть «обратная задача неклассической диффузии». Решение обратных задач возможно только при наличии достаточно эффективной техники решения прямых задач. Несмотря на достаточно длинную историю развития теории дробных производных, как одной из ветвей теории аналитических функций, теория численных методов решения уравнений в дробных производных до последнего времени систематически не развивалась. В настоящее время известны только отдельные, довольно частные результаты, имеющие ограниченную область применимости, поэтому необходимо создание и развитие численных методов.
Численные алгоритмы решения уравнений с дробными производными, предлагаемые в настоящей работе, используются, в частности, при обучении и тестировании нейронных сетей, применяющихся для решения обратной задачи дробной диффузии. В целом, разработка высокоэффективных методов численного решения прямых задач неклассического переноса радионуклидов в сильно неоднородных средах является важнейшей составной частью технологии оценки безопасности подземных захоронений радиоактивных отходов.
Кроме того, вычислительный алгоритм решения двумерного уравнения конвективного переноса, разработанный для прямого численного моделирования стохастической адвекции, показал свою эффективность при решении различных тестовых задач, и пригоден для практического использования.
Личный вклад автора. Соискателем лично:
-
Разработаны и программно реализованы вычислительные алгоритмы решения уравнения диффузии с дробной производной по пространству в одномерном случае. Проведен сравнительный анализ полученных численных методов. Показана аппроксимация и устойчивость. Исследован характер убывания решения на больших расстояниях от источника начальных данных (область «тяжелых хвостов»).
-
Разработаны и программно реализованы вычислительные алгоритмы решения уравнения диффузии с дробной производной по времени в одномерном случае. Проведен сравнительный анализ полученных численных решений с аналитическими фундаментальными решениями и асимптотиками. Исследован характер убывания решения в области далеких «хвостов».
-
Разработан и программно реализован численный метод генерации случайных соленоидальных полей скоростей с заданной корреляцией в двумерном случае.
-
Разработан алгоритм решения прямой задачи конвективного переноса в прямоугольной области, обладающий минимальной сеточной диффузией. Создана компьютерная программа и протестирована с использованием специальных методик (например, с помощью задачи на равномерное вращение). Показано, что предложенный метод не уступает другим высокоточным методикам при меньших вычислительных затратах.
-
Проведено прямое численное моделирование стохастической адвекции радионуклидов в сильно неоднородных геологических формациях. Показано, что результаты прямого численного моделирования стохастической адвекции примеси хорошо согласуются с теоретическими оценками.
Защищаемые положения. На защиту выносятся:
1. Вычислительные алгоритмы решения уравнения диффузии с дробной производной по пространству в одномерном случае с
применением различных подходов к заданию производной дробного порядка.
-
Вычислительные алгоритмы решения уравнения диффузии с дробной производной по времени в одномерном случае.
-
Прямое численное моделирование стохастической адвекции радионуклидов в сильно неоднородных геологических формациях. Вычислительный алгоритм решения прямой задачи конвективного переноса в прямоугольной области, обладающий минимальной сеточной диффузией.
Апробация работы. Основные результаты были представлены на следующих конференциях и семинарах:
-
XIV Всероссийская конференция «Теоретические основы и конструирование численных алгоритмов решения задач математической физики», г. Новороссийск, Абрау-Дюрсо, 19-21 сентября 2002 г.
-
International Conference on "Groundwater in Fractured Rocks", Prague, Czech Republic, 15-19 September 2003.
-
7th US National Congress on Computational Mechanics, New Mexico, USA, 28-30 July 2003.
-
XV Всероссийская конференция «Теоретические основы и конструирование численных алгоритмов для решения задач математической физики», г. Новороссийск, Абрау-Дюрсо, 8-11 сентября 2004 г.
-
7th Hellenic Hydrogeological Conference and 2nd Workshop on Fissured Rocks Hydrology, Athens, Greece, 5-6 October 2005.
-
Ежегодные научные школы-семинары (конференции) ИБРАЭ РАН (2002, 2003, 2004, 2005, 2006).
Публикации. По теме диссертации опубликовано 18 работ, среди которых 9 препринтов ИБРАЭ РАН, 2 статьи в журнале «Известия академии наук», 7 докладов в сборниках трудов всероссийских и международных конференций. Приняты к печати 2 статьи в журналы «Дифференциальные уравнения» и "Transport in Porous Media".
Структура и объём диссертации. Диссертация состоит из введения, четырех глав, заключения, списка литературы и приложения.
Определение Грюнвальда-Летникова
Математический аппарат дробного интегрирования и дифференцирования имеет давнюю историю [89, с. 9-16]. В настоящее время теория дробного интегро-дифференциального исчисления представляет собой одну из ветвей теории функций комплексных переменных и приобретает важное практическое значение в различных областях знаний [90-92]. Математические модели, базирующиеся на дробных производных, находят применение в математической биологии [93], гидрогеологии при моделировании процессов тепло- и массопереноса в сильно неоднородных средах [94-100], задачах уп-ругопластичности, трансзвуковых течений, исследованиях в области полупроводников, в эпидемиологии, финансах и других направлениях.
Дробные интегралы и производные целого порядка — это обычные интегралы и производные. Однако в случае дробного порядка эти понятия имеют своеобразную специфику, которая проявляется, например, в том, что для них в разных ситуациях совершенно естественно возникают их различные модификации. Таким образом, существует несколько подходов к определению производной дробного порядка [89]. Среди них определения Вейля, Грюнвальда-Летникова, Римана-Лиувилля — все они применяются в разработанных и описанных в данной работе вычислительных алгоритмах.
Определение дробной производной в смысле Вейля (Weyt) [89] основано на обобщении оператора дифференцирования в Фурье-пространстве. Если функция ф(х) абсолютно интегрируема на отрезке [-ж, ж], то она представима в виде бесконечной суммы (ряда Фурье): #(х) фке«", фк = -)ф(х)-е-,1аск. — 2л- _л
Оператор первой производной в Фурье-пространстве — это множитель (-/), производная целой степени п — это соответственно {-iky. Приняв, что последнее утверждение справедливо и для нецелых значений п, получим определение дробной производной по Вейлю для абсолютно интегрируемой, 2л- -периодической функции: W) I &кТФкеікх. (1.1) к =-оо
Дробная производная Вейля 2ж -периодической функции — также 2ж -периодическая функция. Здесь (±ik)a =\к\а ехр-±— -signA L
Определение Грюнвальда-Лепгникоеа
Совершенно другой подход используется в определении дробной производной Грюнвальда-Летникова {Grunwald-Letnikov). Рассмотрим конечно-разностные определения целых производных различного порядка: F{x)-F{x-h) F(x) -2F(x-h) + F(x - 2K) h hZ (1.2) F(x)-3F(x-h) + 3F(x-2h)-F(x-3h)
Эти соотношения при /г-»0 дают соответственно первую, вторую и третью производные от F(x). Формулы (1.2) могут быть обобщены на производные произвольного порядка как: A- 0 ft F("\x) = Um F, A?F(x) = (-1)4 i=0 fc Kkj F(x-k-h), (1.3) где биномиальные коэффициенты определяются из выражения: а(а-\)...(а-к + 1) Г(а + 1) Г(к + 1)-Г(а-к + 1)
В выражениях (1.3) бесконечный ряд сходится абсолютно и равномерно при всех а 0 для любой ограниченной функции [89]. Таким образом, производная Грюнвальда-Летникова подчиняется равенству: Z#Kx) = lim . a 0, (1.4) и существует для интегрируемых по Лебегу функций [89].
Определение Римана-Лиувилля Идея еще одного определения дробной производной заключается в обобщении формулы Коши для п -кратного интеграла на дробный порядок п: X „-\ \ (" l)!aJ \ \ ...\ф{х,)ск,...с1х„_{ =__ J(x-/)-V(0« . а а а Отсюда получаем определение дробного интеграла по Риману-Лиувиллю {Riemann-Liouville):
Конечно-разностный метод
Необходимость повышения порядка аппроксимации по пространству вызвана тем, что разностная схема Грюнвальда-Летникова (первого порядка) не достаточно точно описывает поведение профиля концентрации для близких к единице значений параметра дробности a. Увеличение числа узлов (или уменьшение шага расчетной сетки h) частично решает проблему, но приводит к существенному увеличению вычислительных затрат.
В определение Римана-Лиувилля входит интеграл, который может быть численно вычислен с любой заданной точностью. Для получения второго порядка аппроксимации по пространственной переменной дробные потоки, относящиеся к серединам расчетных ячеек, следует аппроксимировать с третьим порядком методом трапеций.
При т = \ операции дробного дифференцирования на отрезке [a,b] по Риману и Лиувиллю имеют вид: дха Г(l-a) dx3a(x) даС -1 a Ь -1-.±)Ші, (2.20) -a) dxht-хУ к даС 1 d xrC(t)dt д(-х)а Г(І-ог) dxj(t-x) хє[а,Ь], 0 ог 1. Представим исходное уравнение дробной диффузии в разностном потоковом виде, подобно (2.10): /-гл+1 - и і і /г + г?" _ + FT" 1 _ /? FT" _ - FT" - _ 1Т А гі+\/2 Гі-\І2 - Х У ГМІ2 Гі-\І2 (2 ОН г 2 h 2 И К J
Выражение (2.21) дает первый порядок аппроксимации по г. Порядок аппроксимации по h будет зависеть от того, с какой точностью мы представим потоки F. Рассмотрим, например, поток F"+U2: +Г,: п= -, 1 « 2. (2.22) Тогда, согласно определению (2.20), после разбиения входящего в него интеграла на сумму интегралов по отрезкам, ограниченным узлами расчетной сетки, получим: +F" = Л+1/2 Г(2-а) dx, +1/2 h)t(xi+mri+ і іхмпrl (2.23) x0 = a, 1 a 2. Совершенно аналогично записывается поток F"+V2: Fn = 1 /+1/2 -1 Г(2-ог) dxMI2 xN_i = b, 1 a 2. ,2 J f C(t)dt Cjt)dt ым Xk (t - xi+V2) XMII (t - xMI2) (2.24) Для получения разностного метода второго порядка аппроксимации по И, достаточно представить подынтегральную функцию С(х) в виде непрерывной, кусочно-линейной функции /Зк-х+/к, где х є іхк,хк+1\, С -С п _ W+i W и Ук=Ск-рк-хк. (2.25)
Вообще говоря, погрешность в определении потоков по формулам (2.23) и (2.24) не превосходит /гт+2, где т — степень многочлена, которым проинтерполирована функция С(х) внутри отрезков [хк, хк+] ]. В нашем случае т = 1, потоки F имеют третий порядок аппроксимации по h, и, следовательно, вся схема (2.21) — второй. Далее, интегрирование в формулах (2.23) и (2.24) можно произвести аналитически. В результате получим довольно громоздкие выражения для обоих потоков, которые удобно представить в операторном виде: L,C =Т (+ +1/2 - Fi-\ll). (2.26) к(сп=-.(- +1/2-- .1/2).
Как показал расчет, вид матриц L и R совершенно аналогичен матрицам, полученным методом Грюнвальда-Летникова в предыдущей главе. Однако коэффициенты вк в операторе (2.14) будут иметь более сложный вид: с _ 1 0 (2-аг).Г(2-аг)-22-" 0,=0о(з2--2), (2.27) вк = 0О [(2k - 2f-a -2-(2k- \f-a + (2 к + \f-a ], к 2, \ а 2.
Смысл коэффициентов вк в методе Римана-Лиувилля, несмотря на более сложный по сравнению с коэффициентами Грюнвальда-Летникова вид, не изменился. Значения коэффициентов Qk при к Ъ достаточно быстро стремятся к нулю. Быстрое убывание коэффициентов Qk наводит на мысль о возможности замены значительной их части нулевыми значениями, что могло бы привести к существенным упрощениям вычислительных алгоритмов. Тестовые расчеты, однако, показывают, что такого рода упрощения могут приводить к качественно неверным результатам.
Прежде всего следует отметить, что механическое обнуление малых коэффициентов в матрицах L и R приведет к нарушению дивергентности разностных уравнений, которое будет выражаться в появлении аппроксима-ционных источников или стоков переносимой величины, полностью меняющих всю картину диффузионного переноса. И этот эффект будет весьма значительным. Чтобы сохранить свойство консервативности уравнений при уменьшении шаблона аппроксимации, следует вернуться к потоковой форме уравнений (2.21) и исследовать последствия отбрасывания малых коэффици 71 ентов в выражениях для соответствующих потоков. Для потока aF"+V2, например, мы имеем: , {(Ь-хЩ -aKv2=-rh- Z Xk-C(xt+k-h,t); h - (2.28) Х0 =1; X, =-(а-1) ; Х2 =(а-1)-(а-2)/2 ;...; Хк+, = - " Из (2.28) видно, что все коэффициенты Хк, начиная со второго, оказываются отрицательными. Можно показать, что коэффициенты Хк имеют следующее свойство: 00 оо А/ Qo=2 =0; Q2=I =a-2; Q2M=2X Q2=a-2; (2.29) 4=0 А=2 i=2 На рис. 9-10 представлено поведение частичных сумм Q"=±0k (2.30) i=0 для методов Грюнвальда-Летникова и Римана-Лиувилля.
Видно, что для достижения соответствующими суммами близких к нулю значений, требуется количество расчетных узлов, превосходящее многие сотни. Простое отбрасывание малых коэффициентов приводит к фиксации значительного дисбаланса в сумме коэффициентов, что не может не сказаться на качественном поведении рассчитываемых функций хотя бы в области «тяжелых» хвостов.
Постановка вычислительной задачи
Как и для задачи Коши, для краевой задачи также можно было бы записать выражение, характеризующее асимптотическое поведение решений при JC — оо, и получить соответствующие рисунки. Однако в настоящей работе Зададимся целью численно решить уравнение диффузии (3.1) с дробной производной по времени. Рассмотрим случаи как «быстрой» (1 / 2), так и «медленной» диффузии (0 / 1). Уравнение будем решать в квадрате x 2л-, 0 t t0. Для этого выберем равномерную расчетную сетку с шагом h по оси х и с шагом г по оси t. Значения концентраций u(x,t) отнесем к узлам расчетной сетки, которые занумеруем от 0 до N-1 вдоль оси х и от 0 до Г вдоль оси /. Проекцию функции u(x,t) на сетку обозначим и1„, где /,weZ, 0 / JV-l, 0 л 7\
Рассмотрим задачу Коши (3.12), выбрав начальное условие в виде сеточной дельта-функции, «пик» которой смещен на середину отрезка 0 х 2л: u(x,t)\=0. = Six-ж). При выборе второго начального условия u ,{x,t)\t=0. (для случая 1 / 2) можно пойти разными путями. Например, положить u t(x,t)\t_Q. =0» как это Делает Майнарди в своей статье [105]. Другой способ — приравнять м (х,0,_0 к некоторой функции, полученной в результате решения уравнения (3.7), которое выведено из квазиволнового представления и отвечает за волну, движущуюся в заданном направлении. Таким образом, на нулевой временной слой проецируется сеточная дельта-функция 8{х-ж), а на первый временной слой (для 1 / 2) — либо та же дельта-функция (если u t(x,t)\t . = 0), либо задается некоторое распределение и[, полученное численно.
Методы аппроксимации уравнения (3.1) могут быть различными. Во-первых, можно использовать разные определения дробной производной по времени: Римана-Лиувилля или Капуто. Не исключено применение других определений, например, Грюнвальда-Летникова (Griinwald-Letnikov) или Маршо (Marcho) [89], но в этой работе они использоваться не будут. Во-вторых, можно по-разному аппроксимировать целую производную по пространству. Методы, базирующиеся на том или ином определении дробной производной, объединены в данной работе в соответствующие параграфы.
Эти методы основаны на определении дробной производной Римана-Лиувилля, введенном в гл. 1. Напомним общий вид правостороннего оператора дифференцирования по Риману-Лиувиллю: +ВЖх) = \ L\ «№ (3.27) где x a, m = 1,2,..., 0 m-\ a m .
Адаптируем это определение к поставленной задаче, положив а = 0 и x = t. Тогда для производной степени у формула (3.27) примет вид: Если 0 у \, 1 d\ +DMt) = —і— \{$-?Гф{&йЬ t 0. (3.28) Если 1 7 2, + / (/) = —і—-iL J(/ - КЖ, t 0. (3.29)
Единственное отличие определений Римана-Лиувилля от определений в форме Капуто (см. уравнения (3.9), (3.10) и (3.11)) заключается в том, что в формулировках Капуто классический оператор дифференцирования целого порядка внесен под знак интеграла и непосредственно действует на рассматриваемую функцию. Если эта функция — константа, то производная Капуто, в отличие от производной Римана-Лиувилля, обращается в нуль. О других особенностях определений Римана-Лиувилля и Капуто написано, например, в статье [109].
Фундаментальные решения для уравнения диффузии с дробной производной по времени, как правило, выводятся на основе определения Капуто (работы [96, 105-107]). Вероятно, аналитическими методами вывести эти решения исходя из определения Римана-Лиувилля достаточно трудно. Так как численные решения уравнения диффузии с дробной производной в смысле Римана-Лиувилля сравнивать пока не с чем, ограничимся случаем «быстрой» диффузии (1 у 2 ).
Начальное условие на первом временном слое
Для того чтобы заставить волну двигаться в одну сторону, выше было предложено задать начальное распределение на первом временном слое и[ таким, чтобы оно удовлетворяло сеточной аппроксимации уравнения (3.7), записанного для функции u(x,t). Можно поступить другим образом: принять уравнение (3.7) в качестве исходного и численно решать его вместо уравнения дробной диффузии (3.1), при этом второе начальное условие и[ не понадобится, но это уже другая задача (она будет рассмотрена ниже в главе «Квазиволновое представление»).
Принимая во внимание определение дробной производной Римана-Лиувилля (3.28), численно аппроксимируем уравнение (3.7). Методы его аппроксимации также могут быть различными. Для простоты остановимся на неявной разностной схеме типа «уголок».
Алгоритм численного решения двумерного уравнения конвективного переноса
Постановка и дискретизация задачи
Рассмотрим систему уравнений конвективного переноса в двумерном случае с дивергентным полем скоростей: да ди р dva . ди dv п .. .. -2- + —— + —— = 0, — + — = 0, и = сЛх,у), v = c(x,y) dt дх ду дх ду Л " Л " (4.10) (x,y)eD:{[0 x L]x[0 y L]}, t є[0,Т] где p = p(x,y,t) характеризует концентрацию переносимой субстанции, сх и с — скорости конвективного переноса вдоль направлений х и у соответственно.
Для прямого численного моделирования процесса стохастической адвекции (конвекции) необходимо использовать алгоритмы численного решения уравнения (4.10), обладающие минимальной сеточной диффузией и улучшенными дисперсионными свойствами.
Покроем область D: {[0 х L] х [0 у L]} равномерной расчетной сеткой с шагами hx=hy = L/(M -1), где М — число узлов по каждому из направлений и будем относить концентрации примеси и конвективные скорости к серединам граней расчетных ячеек. Обозначим их как: Рм/2,; u+i/2 vf+i/2,; "/+1/2» где верхний индекс п указывает на номер временного слоя tn =tn_x+Tn, tQ = 0, т„ 0 — переменный шаг по времени. Кроме того, введем т.н. «консервативную» переменную концентрации, которую будем относить к серединам расчетных пространственно временных ячеек и обозначать как "V/2j+i/2- Величины p?+U2J, p"j+i,2 в дальнейшем будем называть «потоковыми» переменными.
Вычисление новых значений консервативных переменных
В качестве начальных условий в момент времени /0 = 0 зададим значения консервативных переменных i/s?+V2J+i/2 = +,/2,,+,/2 как первичные. Величины потоковых переменных pli,2J, # ;+1/2 при t0 = 0 будем считать вторичными и зададим их следующим образом: ", +1/2=1. о .. . J Vi+U2J= _0 _ K+./2J+./2 / ",,,+,/2 0. 0 к0-./2,,Ч./2 / "/,,+1/2 0 +1/2,,+1/2 {f VMI2J MI2J-\I2 If VMI2J (4.11) На первом шаге по времени аппроксимируем уравнение баланса (4.10) как rfimj.m-V Mnj.m , (ц- ,,+,/2 (ц )1+,/2 , (v tf+,/2,,+, (v )1,/2,/ hr h.. rj2 (4.12) На последующих временных шагах дискретное уравнение баланса будет иметь вид: Kxllj+V2-Kx!Lm , ("- +,,,+,/2-( - ./+,/2 , ( )/+,/2,,+,- )/+,/2., Q (4ЛЗ) И. h.. (г„+г„_,)/2
Нетрудно видеть, что схема (4.13) на гладких решениях аппроксимирует исходное уравнение со вторым порядком аппроксимации по пространственным переменным и первым порядком по времени. В случае постоянства временных шагов будет обеспечен и второй порядок аппроксимации по времени.
Коррекция консервативных переменных Уравнение адвекции (4.10) на соленоидальном поле скоростей можно представить в характеристическом виде: EV + u.?V+v.EV=0 (4.14) dt дх By откуда непосредственно следует, что при переносе по характеристикам значение концентрации остается постоянным. Это значит, что внутри области новых локальных минимумов или максимумов образовываться не может. Этот факт часто называют принципом максимума.
При вычислениях новых консервативных переменных в соответствии с (4.12), (4.13) их значения могут оказаться больше или меньше всех окружающих расчетную ячейку потоковых переменных, что приведет к нарушению принципа максимума и появлению немонотонности в численном решении. Для того чтобы численное решение оставалось монотонным, применим процедуру нелинейной коррекции, заключающуюся в следующем.
Вычислим максимальные и минимальные значения потоковых переменных по расчетной ячейке: (max p)"+ll2j+v2 = max { .+I/2, , ,+I/2, pnMll J, 1/2y+I}
Будем считать консервативные величины, вычисленные по формулам (4.12), (4.13) предварительными и обозначать их как wltvlj n- Рассчитаем величину дисбаланса D" HJJrm в ячейке, образующегося за счет нарушения принципа максимума:
