автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Некоторые математические модели переноса радионуклидов в сильно неоднородных геологических формациях

российская академия наук

Институт Проблем Безопасного Развития Атомной Энергетики

На правах рукописи

Короткин Иван Александрович °'

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

05.13.18 - Математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ

диссертации на соискание ученой степени кандидата физико-математических наук

Москва - 2006

Работа выполнена в Институте проблем безопасного развития атомной энергетики Российской академии наук.

Научный руководитель: д.ф.-м.н., проф.

Головизнин Василий Михайлович

Официальные оппоненты: д.ф.-м.н., проф.

Гасилов Владимир Анатольевич,

д.ф.-м.н., проф.

Иванов Михаил Федорович.

Ведущая организация: Троицкий Институт Инновационных и

Термоядерных Исследований (ТРИНИТИ)

Защита состоится «_»_ 2006 г. в_на заседании диссертационного совета К 002.058.01 при Институте математического моделирования РАН по адресу: 125047, Москва, Миусская пл., 4а.

С диссертацией можно ознакомиться в библиотеке Института математического моделирования Российской академии наук.

Автореферат разослан« С » _ 2006 г.

Ученый секретарь диссертационного совета к.ф.-м.н.

Н.Г. Прончева

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

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

Решение о переработке или окончательном захоронении отработанного ядерного топлива зависит в основном от затрат на ядерный топливный цикл. Результаты экономического анализа замкнутого (с переработкой отработавшего топлива с последующим возвращением в цикл) и открытого (прямое удаление отработавшего топлива без переработки) ядерных топливных циклов в Европе показали некоторое превышение затрат в замкнутом цикле по сравнению с открытым циклом. Тем не менее, ряд стран (Россия, Великобритания, Франция, Япония), учитывая экологическое преимущество (снижение объема и радиоактивности отходов при захоронении после переработки) и возможность рециклирования извлеченных плутония и урана, приняли стратегию замкнутого цикла. США, Канада, Швеция и другие страны предпочли открытый цикл и имеют программы строительства хранилищ для удаления отработанного ядерного топлива в геологических формациях. Начало эксплуатации геологических хранилищ ожидается в США (Юкка-Маунтин) не ранее 2010 г., Швеции — в 2015 г., Финляндии (Олкилуото) — в 2020 г.

Накопление больших объемов радиоактивных отходов и отработанного ядерного топлива с тепловых реакторов поставило вопрос о выработке и реализации долгосрочной стратегии обращения с отходами и в России. Ближайшее время предполагается определить места сооружения крупных могильников для централизованного хранения долгоживущих высокоактивных отходов на территории РФ. Их создание позволило бы решить на самом современном уровне проблему обращения с отходами ядерной энергетики в стране.

Несмотря на интенсивные исследования и большие капиталовложения в разработку и строительство хранилищ, проблема высокоактивных отходов считается сейчас одной из атомрой

библиотека

С.-Пет(р§ург

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

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

примеси на больших временах растет как а убывание кон-

центрации в области «хвостов» (то есть на больших расстояниях, при г » R(t) ) происходит по гауссову закону. Вместе с тем, результаты лабораторных и полевых экспериментов показывают, что для сильно неоднородных неупорядоченных сред, какими являются трещиноватые горные породы, неклассический характер переноса представляется скорее правилом, чем исключением. При этом рост облака частиц примеси со временем (R(t) ) может идти быстрее, чем убывание концентрации при г s> R(t) — медленнее, чем по гауссову закону. Если «хвосты» убывают по степенному закону, то их называют «тяжелыми». Отсюда следует, что предельно допустимые концентрации радионуклидов могут распространяться на расстояния, которые на много порядков превышают найденные из классических представлений. Все это означает, что для получения адекватных оценок надежности радиоактивных захоронений, требуется разработка новых подходов.

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

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

Если дисперсия приращений координат становится бесконечной, то в одномерном случае средняя плотность описывается уравнением дробной диффузии с оператором дробного пространственного дифференцирования, определяемым двумя параметрами, характеризующими функцию распределения. Это параметр а, задающий степенное убывание «хвостов», и параметр /?, задающий степень асимметрии распределения. При а = 2 распределение стремится к нормальному, при а = 1 — переходит в распределение Коши. При а < 1 функция распределения приращений координат блуждающей частицы характеризуется не только бесконечной дисперсией, но и бесконечно большой средней величиной. И хотя, в этом случае, мы не можем получить устойчивых средних для каждой из частиц в отдельности, средние значения плотности частиц в заданной точке пространства остаются устойчивыми переменными и имеют привычный физический смысл.

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

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

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

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

Модель стохастической адвекции представляет собой яркий пример того, как во фрактальной, сильно неупорядоченной среде может реализоваться режим дробной диффузии (супердиффузии) на умеренных расстояниях с одновременным отсутствием «тяжелых» степенных «хвостов» в распределении концентрации на больших расстояниях.

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

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

Для модели дробной диффузии разработаны новые вычислительные алгоритмы решения уравнений диффузии с дробной производной по пространству и с дробной производной по времени в одномерном случае. Предложены способы обобщения алгоритмов на многомерный случай.

Для модели стохастической адвекции разработан численный алгоритм генерации случайных соленоидальных двумерных полей скоростей и создан новый численный метод решения задачи конвективного переноса примеси в прямоугольной области. Данный метод обладает минимальной сеточной диффузией, улучшенными дисперсионными свойствами и является модернизацией и обобщением одномерного алгоритма «КАБАРЕ». В ходе совершенствования двумерной схемы «КАБАРЕ» разработаны и апробированы на тестовых задачах новые алгоритмы коррекции консервативных и потоковых величин.

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

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

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

Личный вклад автора. Соискателем лично:

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

2. Разработаны и программно реализованы вычислительные алгоритмы решения уравнения диффузии с дробной производной по времени в одномерном случае. Проведен сравнительный анализ полученных численных решений с аналитическими фундаментальными решениями и асимптотиками. Исследован характер убывания решения в области далеких «хвостов».

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

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

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

Защищаемые положения. На защиту выносятся: 1. Вычислительные алгоритмы решения уравнения диффузии с дробной производной по пространству в одномерном случае с

применением различных подходов к заданию производной дробного порядка.

2. Вычислительные алгоритмы решения уравнения диффузии с дробной производной по времени в одномерном случае.

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

Апробация работы. Основные результаты были представлены на следующих конференциях и семинарах:

1. XIV Всероссийская конференция «Теоретические основы и конструирование численных алгоритмов решения задач математической физики», г. Новороссийск, Абрау-Дюрсо, 19-21 сентября 2002 г.

2. International Conference on "Groundwater in Fractured Rocks", Prague, Czech Republic, 15-19 September 2003.

3. 7th US National Congress on Computational Mechanics, New Mexico, USA, 28-30 July 2003.

4. XV Всероссийская конференция «Теоретические основы и конструирование численных алгоритмов для решения задач математической физики», г. Новороссийск, Абрау-Дюрсо, 8-11 сентября 2004 г.

5. 7th Hellenic Hydrogeological Conference and 2nd Workshop on Fissured Rocks Hydrology, Athens, Greece, 5-6 October 2005.

6. Ежегодные научные школы-семинары (конференции) ИБРАЭ РАН (2002, 2003, 2004, 2005, 2006).

Публикации. По теме диссертации опубликовано 18 работ, среди которых 9 препринтов ИБРАЭ РАН, 2 статьи в журнале «Известия академии наук», 7 докладов в сборниках трудов всероссийских и международных конференций. Приняты к печати 2 статьи в журналы «Дифференциальные уравнения» и "Transport in Porous Media".

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

СОДЕРЖАНИЕ РАБОТЫ

Во введении приведен краткий обзор результатов лабораторных и полевых экспериментов, описанных в зарубежной литературе.

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

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

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

Стандартное уравнение конвекции-диффузии в одномерном случае выглядит следующим образом:

дС дС пд2С ч

где С — концентрация переносимой субстанции, V — скорость конвективного переноса вдоль оси х, £> — коэффициент диффузии (дисперсии).

Замена в уравнении (1) второй пространственной производной на дробную (степени а, где а е [1,2]) приводит к одномерному уравнению диффузии с дробной производной по пространству. Обобщение первой временной производной в уравнении (1) на производную

дробной степени у, где у е (0,2], приводит к уравнению диффузии с дробной производной по времени в одномерном случае.

Глава состоит из 3 параграфов.

В первом параграфе первой главы введено понятие дробной производной, приведены определения Вейля, Римана-Лиувилля и Грюн-

вальда-Летникова, показана их эквивалентность на достаточно «хороших» функциях. Приведены примеры производных дробного порядка от некоторых элементарных функций.

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

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

Глава 2 диссертации посвящена описанию разработанных и протестированных численных методов решения уравнения дробной диффузии с дробной производной по пространству:

дС дС /)[,, -V даС „ _ч даС С = С(х,/),

где а е [1,2] — степень дробной производной, характеризующая неоднородность среды, в которой осуществляется перенос примеси, Р е [—1,1] — параметр «скошенности», характеризующий пространственную асимметрию среды. При а - 2 уравнение дробной диффузии переходит в классическое диффузионное уравнение.

Роль показателя а в случае 1 < а < 2 при оценке концентрации на больших расстояниях и, следовательно, безопасности захоронений состоит в следующем. Пусть Т — время, за которое описываемая уравнением дробной диффузии концентрация от мгновенного точечного источника на заданном единичном расстоянии (используем безразмерные переменные) достигает уровня Ю-8 от максимального. Отметим, что даже столь малые концентрации фактически могут оказаться опасными из-за превышения санитарных норм. Время Т сильно зависит от порядка а при прочих равных условиях (считаем дисперсию равной 1, это означает, что основное ядро распределения за единичное время достигает единичной ширины). Даже небольшое отклонение параметра а от значения а = 2 (классическая диффузия) приводит к уменьшению времени Т на 2-3 порядка величины. Чтобы достичь заданной концентрации за то же время при классиче-

ской диффузии, потребовалось бы соответственно увеличить коэффициент диффузии. Зависимость увеличенного таким образом коэффициента диффузии от параметра а показана на рис. 1. Наиболее сильна эта зависимость при а —> 2.

а

Рис. 1. Зависимость эффективного коэффициента классической диффузии, необходимого для обеспечения того же времени достижения концентрации 10-8 в заданной точке, что и при дробной диффузии, от порядка дробной производной а

Существенным отличием переноса, описываемого уравнением дробной диффузии, от классического является степенной характер убывания концентрации C{x,t) диффундирующего вещества на бесконечности. В классической диффузии концентрация убывает по экспоненте, поэтому при анализе безопасности подземных захоронений радионуклидов принципиально важным становится вопрос о возможности неклассического (супердиффузионного) переноса.

На рис. 2 в логарифмическом масштабе, приведены графики фундаментального решения при различных величинах а. Хорошо видно резкое увеличение концентрации на «хвостах» при аФ 2 по сравнению со случаем а = 2.

ТаМв: 1 = 0.1; Ьега = 0.0

-а!рЬа = 2 0 || - а!рЬа = 17!

Рис. 2. Вид фундаментального решения классической диффузии (а = 2) и уравнения дробной диффузии при а = 1.7

Глава состоит из 7 параграфов.

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

В третьем параграфе описан вычислительный алгоритм решения уравнения дробной диффузии, получивший название «метод Фурье». Данный алгоритм использует определение дробной производной по Вейлю, которое основано на обобщении оператора дифференцирования в Фурье-пространстве и для абсолютно интегрируемой, 2ж -периодической функции выглядит следующим образом:

(3)

к~—

Методом Фурье можно получить только периодическое решение С(х, Г) = С(х + 2л, г) уравнения (2).

В четвертом параграфе рассмотрен конечно-разностный метод решения уравнения дробной диффузии, основанный на определении производной дробного порядка Грюнвальда-Летникова (а > 0 ):

(

а

F(x-^t•Л)• (4)

V"-/

Данное определение является обобщением конечно-разностных определений целых производных.

Метод Грюнвальда-Летникова свободен от условия периодичности, которое накладывается на решение, если его искать методом Фурье. Предложенная схема аппроксимирует исходное уравнение (2) с первым порядком по пространству и времени и является устойчивой при ат I И* < 1, где г и А — шаги равномерной расчетной сетки по времени и пространству соответственно. Разработан и описан безусловно устойчивый, частично неявный вариант данной схемы.

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

где т = 1,2,... и 0<т-\<а<т.

Метод Римана-Лиувилля обладает вторым порядком аппроксимации по пространству и первым по времени. Предложены две схемы: устойчивая при ат ! И* < 1 и безусловно устойчивая.

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

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

Глава 3 посвящена описанию методики численного решения уравнения диффузии с дробной производной по времени:

где у — степень оператора дробного дифференцирования по времени.

(5)

(6)

В отличие от уравнения диффузии с дробной производной по пространству, в данной модели оператор дробного дифференцирования по времени является правосторонним, производная вычисляется на отрезке [0, /]. При у —>1 уравнение (6) переходит в классическое уравнение диффузии с экспоненциальным затуханием решения на бесконечности. Для его вывода необходимо одно начальное условие: С(х,0|,=0- При у2 уравнение (6) стремится к в волновому

уравнению, которому требуется два начальных условия: С(л:,/)|, и

Глава состоит из 5 параграфов.

В промежуточном случае, когда 0 < у < 2, возникает вопрос о количестве начальных данных, необходимых для корректной постановки задачи Коши. В первом параграфе вводится так называемое квазиволновое представление уравнения (6), данное представление позволяет не только решить вопрос о количестве начальных данных, но и ввести коэффициент «кососимметричности» решения.

Фундаментальные решения уравнения диффузии с дробной производной по времени приведены во втором параграфе. В третьем параграфе приводится постановка вычислительной задачи.

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

В пятом параграфе за основу взято распространенное определение дробной производной в смысле Капуто:

=-1-, 0 (7)

Г{т-у)1К 4 к ¿Г

где т = 1,2,... и 0<т-1<у<т.

Дробная производная по Капуто является простой адаптацией определения Римана-Лиувилля к рассматриваемой задаче, что позволяет аналитически вывести фундаментальные решения уравнения (6).

Одна из основных особенностей определения Капуто — равенство нулю производной любого порядка от константы.

Продемонстрированы сходства и различия численных фундаментальных решений уравнения (6) при использовании в расчетах определений дробной производной Римана-Лиувилля и Капуто. На рис. 3 показаны численные решения уравнения (6) в моменты времени / = 0.1 и (= 0.5 в логарифмическом масштабе в сравнении с соответствующими классическими гауссовыми распределениями. Для уравнения диффузии с дробной производной по времени (как Римана-Лиувилля, так и Капуто) характерно отсутствие «тяжелых» степенных «хвостов».

ватта = 1.2

Рис. 3. Графики решений уравнения дробной диффузии различными методами при у = 1.2 в логарифмическом масштабе в сравнении с классическими распределениями

Глава 4 посвящена прямому численному моделированию задачи о стохастической адвекции примеси в сильно неоднородных геологических средах. Математически такой процесс в случае двух пространственных измерений можно записать как:

^- + &\(с<р) = 0, (Ну С = 0, (8)

где (р — концентрация переносимой консервативной примеси, с — вектор случайного поля адвекционных скоростей, являющийся

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

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

Глава состоит из 3 параграфов.

В первом параграфе предлагается вычислительный алгоритм генерации случайных соленоидальных полей скоростей с заданной корреляцией в двумерном случае. На данном этапе принципиальным является вопрос об установлении связи между полученными случайными полями и распределениями сетей трещин в горных породах. Эту проблему можно назвать задачей конструирования «электронного» представления (по аналогии с понятием «электронного документа») изучаемой неоднородной геологической среды.

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

Представленный в диссертации двумерный алгоритм решения уравнения конвективного переноса является обобщением одномерного алгоритма «КАБАРЕ». Данная схема включает процедуры нелинейной коррекции консервативных и потоковых переменных. Эти процедуры делают схему «КАБАРЕ» монотонной в областях сильных градиентов и не нарушают второго порядка аппроксимации в

областях гладкости. Схема является явной и условно устойчивой при |гг|<1, <1. Можно показать, что в схеме «КАБАРЕ» аппрокси-мационная диффузия является наименьшей из всех возможных.

На рис. 4 представлен пример численного решения задачи стохастической адвекции на одной из реализаций случайного поля скоростей в некоторый момент времени. В начальный момент в центре расчетной области задано единичное возмущение. На графике значения концентраций в закрашенных серым цветом областях в плоскости г - О составляют 1 (Г6 -И О-8. Ядро распределения в данный момент еще сосредоточено в центральных ячейках, но значимые значения концентрации уже достигли границ расчетной области.

Pite. 4. Пример решения задачи стохастической адвекции

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

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

размерности 256x256 и вычислительных алгоритмов с минимально возможной аппроксимационной диффузией прямое стохастическое моделирование позволяет выявлять основные качественные особенности, присущие переносу в сильно неоднородных средах. При переходе к более подробным сеткам 512x512 отклонение расчетных результатов от теоретических уменьшается примерно в два раза.

Приложение диссертации состоит из двух параграфов, которые посвящены краткому описанию специальных функций Миттаг-Леффлера (Mittag-Leffler) и Райта (Wright). Некоторые свойства данных функций применялись при выводе фундаментальных решений уравнения диффузии с дробной производной по времени.

ОСНОВНЫЕ РЕЗУЛЬТАТЫ ДИССЕРТАЦИИ

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

2. Разработаны вычислительные методы и создан комплекс программ для решения уравнения диффузии с дробной производной по времени в одномерном случае. За основу предложенных алгоритмов приняты определения Римана-Лиувилля и Капуто. Проведено сравнение результатов численных расчетов в случае Капуто с фундаментальными решениями. Показано отсутствие «тяжелых» степенных «хвостов» в полученных решениях.

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

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

СПИСОК ПУБЛИКАЦИЙ ПО ТЕМЕ ДИССЕРТАЦИИ

1. В.М. Головизнин, В.П. Киселев, И.А. Короткин, Ю.И. Юрков, «Некоторые особенности вычислительных алгоритмов для уравнения дробной диффузии». Препринт №ГВЯАЕ-2002-01, Москва: ИБРАЭ РАН, 2002. 57 с.

2. В.М. Головизнин, В.П. Киселев, И.А. Короткин, «Численные методы решения уравнения дробной диффузии в одномерном случае». Препринт №1В11АЕ-2002-10, Москва: ИБРАЭ РАН, 2002. 35 с.

3. В.М. Головизнин, В.П. Киселев, В.Н. Семенов, И.А. Короткин, А.Г. Хромов, Ю.И. Юрков, «Исследование возможностей практического использования искусственных нейронных сетей для идентификации основных параметров дробной диффузии». Препринт ЖВКАЕ-2002-12, Москва: ИБРАЭ РАН, 2002. 37 с.

4. В.М. Головизнин, В.П. Киселев, И.А. Короткин, Ю.И. Юрков, «Решение обратной задачи по идентификации основных параметров дробной диффузии». Препринт №1ВЯАЕ-2002-20, Москва: ИБРАЭ РАН, 2002.20 с.

5. В.М. Головизнин, В.П. Киселев, И.А. Короткин, «Численные методы решения уравнения диффузии с дробной производной по времени в одномерном случае». Препринт №ГОЯАЕ-2003-12, Москва: ИБРАЭ РАН, 2003. 35 с.

6. В.М. Головизнин, В.П. Киселев, В.И. Питербарг, И.А. Короткин, Ю.И. Юрков, «Генерация случайных полей трещиноватости с заданной корреляционной функцией и заданным средним». Препринт №1В11АЕ-2003-17, Москва: ИБРАЭ РАН, 2003. 23 с.

7. В.М. Головизнин, В.П. Киселев, И.А. Короткин, Ю.И. Юрков, «Решение задач диффузионного переноса на случайных полях трещиноватости и дробная диффузия». Препринт №1ВЯАЕ-2004-01, Москва: ИБРАЭ РАН, 2004. 28 с.

8. В.М. Головизнин, В.П. Киселев, И.А. Короткий, «Численные методы решения уравнения дробной диффузии в одномерном случае». Сборник тезисов докладов XIV Всероссийской конференции «Теоретические основы и конструирование численных алгоритмов решения задач математической физики», г.Новороссийск, Абрау-Дюрсо, 19-21 сентября 2002 г. С. 101-104.

9. В.М. Головизнин, В.П. Киселев, В.Н. Семенов, И.А. Короткин, Ю.И. Юрков, «Решение обратных задач для уравнения дробной диффузии с помощью искусственных нейронных сетей». Сборник тезисов докладов XIV Всероссийской конференции «Теоретические основы и конструирование численных алгоритмов решения задач математической физики», г. Новороссийск, Абрау-Дюрсо, 19-21 сентября 2002 г. С. 104-1 Об.

10. V.M. Goloviznin, V.P. Kiselev, V.N. Semenov, I.A. Korotkin, Yu.I. Yurkov, "Solution methods for direct and inverse problems of fractional dispersion". Proceedings of 7th US National Congress on Computational Mechanics, New Mexico, USA, 28-30 July 2003.

11. Ivan A. Korotkin, Yuriy I. Yurkov, Vladimir P. Kiselev, Va-siliy M. Goloviznin, Vladimir N. Semenov, "Solution of the inverse problem of identification of fractional diffusion parameters". Proceedings of the International Conference on "Groundwater in Fractured Rocks", Prague, Czech Republic, 15-19 September 2003. Extended Abstracts. Edited by J. Krasny, Z. Hrkal and J. Bruthans. IHP-VI, Series on Groundwater No. 7.

12. JI.A. Большое, В.М. Головизнин, A.M. Дыхне, В.П. Киселев, П.С. Кондратенко, В.Н. Семенов, «Новые подходы к оценке безопасности захоронений радиоактивных отходов». Журнал «Известия РАН. Энергетика». 2004. №4. С. 99-108.

13. В.М. Головизнин, В.П. Киселев, И.А. Короткин, Ю.И. Юрков, «Прямые задачи неклассического переноса радионуклидов в геологических формациях». Журнал «Известия РАН. Энергетика». 2004. №4. С. 121-130.

14. В.М. Головизнин, В.П.Киселев, В.Н.Семенов, И.А. Короткин, Ю.И. Юрков, «Решение обратной задачи неклассического переноса радионуклидов в геологических формациях с использованием нейронных сетей». Журнал «Известия РАН. Энергетика». 2004. №5. С. 81-87.

15. В.М. Головизнин, В.П. Киселев, И.А. Короткин, «Решение задач диффузионного переноса на двумерных случайных полях тре-щиноватости и дробная диффузия». Сборник тезисов докладов XV Всероссийской конференции «Теоретические основы и конструирование численных алгоритмов для решения задач математической физики», г. Новороссийск, Абрау-Дюрсо, 8-11 сентября 2004 г. С. 29-30.

16. В.М. Головизнин, М.К. Мелькина, И.А. Короткин, «Решение уравнения дробной диффузии методом Монте-Карло». Сборник тезисов докладов XV Всероссийской конференции «Теоретические основы и конструирование численных алгоритмов для решения задач математической физики», г. Новороссийск, Абрау-Дюрсо, 8-11 сентября 2004 г. С. 30-31.

17. V.M. Goloviznin, P.S. Kondratenko, L.V. Matweev, V.N. Semenov, I.A. Korotkin, "Direct Numerical Modeling of Stochastic Radionuclide Advection in Highly Non-Uniform Media". Preprint №IBRAE-2005-01. Moscow: NSI RAS, 2005. 37 p.

18. В.М. Головизнин, O.C. Сороковикова, М.К. Мелькина, И.А. Короткин, «Стохастический подход к моделированию распространения примеси в трещиноватых и сильно неоднородных средах. Связь стохастического моделирования и подходов, использующих уравнения дробной диффузии». Препринт №IBRAE-2005-02, Москва: ИБРАЭ РАН, 2005. 34 с.

19. V.M. Goloviznin, V.N. Semenov, I.A. Korotkin, "Numerical modeling of stochastic radionuclide advection in fissured media" (full paper). 7th Hellenic Hydrogeological Conference and 2nd Workshop on Fissured Rocks Hydrology, key lectures and workshop proceedings. Eds.: G. Stournaras, K. Pavlopoulos, Th. Bellos. Vol. 2, P. 127-131. Athens, 2005.

3.Р06 fi

»- 95 61

Введение.

Глава 1. Элементы дробного интегро-дифференциального исчисления. Дробная диффузия.

1.1. Понятие дробной производной.

Определение Вейля.

Определение Грюнвальда-Летникова.

Определение Римана-Лиувилля.

Примеры производных дробного порядка.

Дробный Лапласиан.

1.2. Модель дробной диффузии.

1.3. Вывод уравнения дробной диффузии.

Способ Эйлера.

Распределение Леви.

Глава 2. Методы численного решения уравнения диффузии с дробной производной по пространству.

2.1. Постановка задачи.

2.2. Фундаментальные решения.

2.3. Метод Фурье.

2.4. Конечно-разностный метод.

2.5. Метод сплайн-аппроксимаций.

2.6. Анализ численных методов.

Устойчивость явной разностной схемы.

Стационарная задача. Аппроксимация.

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

Глава 3. Методы численного решения уравнения диффузии с дробной производной по времени.

3.1. Квазиволновое представление уравнения диффузии с дробной производной по времени.

3.2. Фундаментальные решения.

Задача Коши (Cauchy Problem).

Краевая задача (Signalling Problem).

3.3. Постановка вычислительной задачи.

3.4. Методы, основанные на определении Римана-Лиувилля.

Начальное условие на первом временном слое.

Вычислительные особенности.

Квазиволновое представление.

3.5. Методы, основанные на определении Капуто.

Начальное условие на первом временном слое.

Вычислительные особенности.

Квазиволновое представление.

Глава 4. Численное моделирование стохастической адвекции радионуклидов в сильно неоднородных средах.

4.1. Генерация случайных соленоидальных полей скорости с заданной корреляцией в случае двух пространственных измерений.

Результаты работы алгоритма.

4.2. Алгоритм численного решения двумерного уравнения конвективного переноса.

Постановка и дискретизация задачи.

Вычисление новых значений консервативных переменных.

Коррекция консервативных переменных.

Вычисление новых значений потоковых переменных.

Тестирование алгоритма.

4.3. Задача о стохастической адвекции примеси.

Случайные поля скоростей.

Распространение примеси с течением времени.

Тяжелые хвосты».

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

Долгосрочное наращивание производства энергии за счет увеличения потребления традиционных источников органического топлива возможно, в первую очередь, за счет увеличения относительной и абсолютной добычи угля, поскольку мировых запасов нефти и газа при нынешнем темпе потребления хватит лишь на несколько десятилетий. Этот путь потребует изменения соотношения потребления природных ресурсов в большинстве развитых стран в пользу угля. Но при этом возникают три серьезные проблемы: транспортная (запасы угля распределены территориально крайне неравномерно), экологическая (обусловленная накоплением отходов — золы, выбросами в атмосферу вредных газов и возрастающими темпами сжигания кислорода), техническая (необходимость приемлемого решения проблемы получения жидкого топлива и газа из угля).

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

Решение энергетических проблем человечества возможно на основе развития ядерной энергетики. Как известно, имеются две принципиальные возможности получения ядерной энергии: во-первых, в реакции деления, во-вторых, в реакции синтеза. Ядерная энергетика сегодняшнего дня развивается только на базе АЭС с реакторами деления, главным образом, с реакторами на тепловых нейтронах, использующими в качестве ядерного горючего только один изотоп урана — уран-235. Другой изотоп природного урана, уран-238, составляющий основную его часть (более 99,3%), используется незначительно. С учетом не только надежно известных, но и предполагаемых ресурсах урана, ядерная энергетика с реакторами на тепловых нейтронах обеспечена энергоресурсами на 50-100 лет, то есть примерно по запасам сопоставима с запасами энергетики на органическом топливе и не может обеспечить создание мировой энергетики на длительный (в несколько столетий) период. Однако в настоящее время существует несколько возможных стратегий развития топливных ресурсов ядерной энергетики [1].

В основе первой группы стратегий лежит усовершенствование реализованных на сегодняшний день так называемых открытых топливных циклов без регенерации ядерного горючего, основанных на использовании урана-235 (совместно с ураном-238 и торием-232). Топливным циклом обычно называют всю последовательность повторяющихся производственных процессов в топливно-энергетическом комплексе, начиная с добычи топлива, включая производство энергии и кончая удалением отходов. Основное направление научно-технического поиска в развитии первой группы стратегий связано с повышением глубины выгорания делящегося материала (как урана-235, так и образующегося в процессе облучения плутония и урана-233). Основной недостаток первой группы стратегий — увеличение затрат и масштабов добычи сырья, необходимость строительства обогатительных заводов и специальных хранилищ для длительного хранения отработанного топлива.

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

Полномасштабное развитие ядерной энергетики может решить проблемы покрытия мирового топливно-энергетического баланса в течение многих столетий без каких-либо ограничений, естественно при условии обеспечения безопасной работы АЭС и успешного решения возникающих экологических проблем, прежде всего — проблемы захоронения радиоактивных отходов.

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

Главное возражение противников ядерной энергетики основано на том, что при работе реакторов АЭС в ядерном топливе образуется значительное количество продуктов деления урана и плутония (осколков). Продукты деления классифицируются по периоду полураспада Г//2 и виду излучения. Существуют короткоживущие и долгоживущие продукты деления. Короткоживущие имеют значения Тщ в диапазоне от нескольких секунд до десятков часов. Примерно через 10 значений Тщ они практически полностью распадаются. Нас больше интересуют долгоживущие продукты деления, представляющие значительную радиационную опасность. Кроме того, при эксплуатации АЭС, как и в любом промышленном производстве, помимо продуктов деления образуется некоторое количество отходов. Значительная часть этих отходов имеет наведенную активность, возникшую в результате активации нейтронами конструкционных материалов, теплоносителя и воздуха, заполняющего шахту реактора. Эти отходы возникают в связи с необходимостью замены различных внутрикорпусных устройств датчиков контроля с наведенной активностью, замены фильтров водяного теплоносителя I контура, содержащих радиоактивные продукты коррозии, дезактивации загрязненного оборудования и производственных помещений и т.д. Все это ведет к появлению твердых, жидких и газообразных радиоактивных отходов. Таким образом, в результате работы АЭС имеются отработавшее ядерное топливо с высокой радиоактивностью и радиоактивные отходы [5-6].

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

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

Решение о переработке или окончательном захоронении отработанного ядерного топлива зависит в основном от затрат на ядерный топливный цикл. Результаты экономического анализа замкнутого (с переработкой отработавшего топлива с последующим возвращением в цикл) и открытого (прямое удаление отработавшего топлива без переработки) ядерных топливных циклов в Европе показали некоторое превышение затрат в замкнутом цикле по сравнению с открытым циклом. Тем не менее, ряд стран (Россия, Великобритания, Франция, Япония), учитывая экологическое преимущество (снижение объема и радиоактивности отходов при захоронении после переработки) и возможность рециклирования извлеченных плутония и урана, приняли стратегию замкнутого цикла. США, Канада, Швеция и другие страны предпочли открытый цикл и имеют программы строительства хранилищ для удаления отработанного ядерного топлива в геологических формациях. Начало эксплуатации геологических хранилищ ожидается в США (Юкка-Маунтин) не ранее 2010 г., Швеции — в 2015 г., Финляндии (Олкилуото) — в 2020 г.

Накопление больших объемов радиоактивных отходов и отработанного ядерного топлива с тепловых реакторов поставило вопрос о выработке и реализации долгосрочной стратегии обращения с отходами и в России. Ближайшее время предполагается определить места сооружения крупных могильников для централизованного хранения долгоживущих высокоактивных отходов на территории РФ. Их создание позволило бы решить на самом современном уровне проблему обращения с отходами ядерной энергетики в стране.

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

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

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

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

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

Лабораторные эксперименты. Хорошо известно, что процесс проникновения влаги в пористые материалы за счет капиллярных сил носит диффузионный характер и описывается уравнением диффузии, в котором коэффициент диффузии, вообще говоря, есть функция насыщенности. В однородной проницаемой среде этот коэффициент не зависит от координат. В этом случае в отсутствие влияния гравитационных сил (речь идет об образцах небольшого размера и/или горизонтальном направлении проникновения жидкости) и при соответствующих граничных условиях на поверхности образца (5 = 0 при г = О и 5 = 1 при г >0) насыщенность 5 в образце может зависеть от координаты и времени при одномерном распространении только как s(x//1/2). В частности, полное количество жидкости, впитавшейся в образец за время t, будет пропорционально /1/2.

В неоднородной пористой среде, какой фактически является любой реальный материал, характер распространения жидкости может существенно измениться. Имеются данные специальных лабораторных экспериментов по пропитке водой таких материалов, как кирпич и известняк [8]. Предварительно высушенные образцы приводились в контакт с водой одной из поверхностей, и пространственный профиль насыщенности внутри образца прецизионно измерялся с помощью ЯМР-метода. Результаты измерений показали, что в упомянутых реальных материалах диффузия воды имеет аномальный характер. Экспериментальная зависимость насыщенности от координаты и времени может быть представлена как функция от x/ta, где а ф 0.5. Оцененное по экспериментальным данным значение показателя оказалось а »0.6. Таким образом имеет место «супердиффузия» при впитывании влаги в реальные материалы.

В работе [9] дана численная интерпретация этих экспериментальных результатов. Моделирование диффузионного впитывания осуществлялось с помощью двумерной модели сеточных автоматов (lattice gas automaton). Принципиальным для использованной модели является учет неоднородности среды — в модели однородной среды аномальной диффузии не наблюдалось. Наряду с неоднородностью, по утверждению авторов, принципиальным является также сильная зависимость проницаемости среды от насыщенности. При этом «супердиффузия» наблюдается в случае, если проницаемость растет с увеличением насыщенности (3D/8s > 0). Этот случай реализовывался в описанных эксперментах и в расчетах [8, 9]. В противоположном случае (dD/ds < 0) должна иметь место «субдиффузия», то есть а <0.5. В качестве иллюстрации такого вывода приводятся результаты диффузионного эксперимента в условиях dD/ds < 0 [10].

Природа наблюдаемой аномальной диффузии воды в пористые материалы детально не ясна. Можно отметить, что модельная неоднородность среды, которая вводилась в численные расчеты, не имела фрактального характера. К сожалению, длительность наблюдения в экспериментах и в расчетах не позволяет утверждать, что «супердиффузионный» режим с а > 0.5 является предельным режимом при больших хи/.Ив экспериментах с обоими материалами, и в численных расчетах величина а, вообще говоря, не оставалась постоянной, и наблюдалась тенденция к уменьшению а и приближению к 1/2. Но даже в том случае, если наблюдавшийся в экспериментах режим является переходным, учет существования такого режима принципиально важен для оценки скорости проникновения влаги в реальные материалы.

В [11] исследовалась адвекция-диффузия в искусственной пористой среде, представляющей насыпку мелких стеклянных пористых шариков. Вода прокачивалась сквозь трубку, заполненную такой средой, с постоянной скоростью, и в течение некоторого времени прослеживалась эволюция распределения концентрации протонов, поляризованных в момент / магнитным импульсом. Профиль концентрации поляризованных протонов измерялся методом ЯМР через определенные промежутки времени в течение интервала, меньшего времени релаксации поляризации. В определенных режимах течения, при промежуточных числах Пекле, обнаружено негауссовское поведение профиля концентрации, что связано с сильной неоднородностью скорости течения между шариками и в порах внутри шариков.

Можно отметить еще ряд экспериментов, в которых наблюдалось не-фиковский (поп^Ыап) характер диффузии растворенных в жидкости трассеров при фильтрации ее через сильно неоднородные пористые среды. К ним можно отнести эксперимент [12] по изучению фильтрации через искусственную неоднородную среду в виде слоя стеклянных шариков разных размеров, неоднородно перемешанных.

Большое количество лабораторных экспериментов проведено с целью исследования протекания через отдельные насыщенные водой натуральные трещины [13-18]. В работах [19-21] исследовалось протекание воды через ненасыщенные образцы натуральных трещин. В [21] использовалась реплика натуральной трещины из прозрачного материала, что позволило визуализировать течение и исследовать его в естественном режиме неполного насыщения, возникающем при подаче воды на входное сечение трещины.

Течение воды в реальных трещинах сильно отличается от течения между двумя гладкими параллельными плоскостями. Сильно нерегулярная поверхность стенок трещины может приводить к анизотропии потока, а также к нарушению известного «кубического закона». При изменении расстояния между двумя стенками общая проницаемость образца меняется быстрее, чем пропорционально кубу расстояния, как это должно быть для гладких плоскостей [15]. Нерегулярная шероховатость стенок трещин приводит к тому, что даже в малых масштабах (5-40 см) поток проявляет стохастические свойства, в нем наблюдаются предпочтительные пути, когда до 90% потока проходит через 5-20% сечения на выходе [21, 22]. В этих условиях можно ожидать аномальную дисперсию переносимых примесей [23].

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

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

Важное место в исследованиях свойств водоносных слоев (aquifer) занимает изучение характеристик самой вмещающей пористой среды. Для трещиноватых сред такие исследования проводятся как в лабораторных условиях, так и непосредственно на местах.

Исследования структуры системы трещин. В последнее время появилось много сообщений об исследовании геометрических характеристик сети трещин непосредственно, наблюдая следы выхода трещин на поверхность трещиноватой породы, либо изучая пересечения трещин с внутренней поверхностью скважин с помощью специальной техники [25, 26].

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

В обзоре [26] показано, что наблюдаемые сети трещин в широком диапазоне масштабов усреднения обладают фрактальными свойствами. Приведенные в [26] гистограммы частоты появления фрактального индекса показывают, что для двумерных сечений трехмерных сетей трещин значения индекса могут пробегать все значения от 1 до 2. Наиболее часто встречающиеся значения близки к 1.5, либо к 2. Первое значение соответствует либо однородным сеткам трещин, либо перколяционным кластерам в двумерной геометрии. Значение из второго диапазона соответствует двумерным кластерам, являющимся проекцией трехмерных перколяционных кластеров [27].

Теория перколяции [27-30] является наиболее подходящим инструментом для исследования сетей трещин: благодаря фрактальной структуре пер-коляционных кластеров они являются хорошей моделью сети трещин.

Важнейшим с точки зрения водопроницаемости свойством сети трещин в непроницаемых (или почти непроницаемых) породах является их связность, взаимные соединения. Даже густая сеть трещин может быть непроницаемой для влаги, если отдельные трещины не соединены между собой.

В теории перколяции свойства проницаемости определяются порогом перколяции. Особый интерес в этой связи представляют состояния перколя-ционной среды вблизи порога перколяции. Во многих случаях, по крайней мере для трещин, возникших от разгрузки напряжений в блоках горных пород, есть основания считать, что эта сеть находятся в состоянии, близком к порогу (с проводящей стороны от порога). В механике разрушения предполагается, что разгрузка напряжений происходит в момент, когда сеть трещин, имеющаяся в сплошном блоке, достигает порога перколяции [30].

В трещиноватых породах, как правило, имеет место очень большой разброс локальной проницаемости в пределах системы трещин. Например, измерения, проведенные в гранит-урановом руднике Рапау-Аи§егез (Франция) более чем в 100 точках в пределах 100 м, показали, что разброс локальной проницаемости в разных точках составляет 4 порядка величины [31]. При этом, несмотря на то, что вся сеть геометрически хорошо соединена, только 0.1% всех трещин пропускают большую часть потока между разными точками. Аналогичный результат описан в [24]. Поэтому даже хорошо связная сеть при условии широкого распределения проницаемости отдельных трещин фактически оказывается вблизи порога перколяции [32].

В последнее время достигнут значительный прогресс в исследовании свойств поверхностей, ограничивающих реальные трещины. Они существенно отличаются от плоских параллельных поверхностей, какие в первом приближении используются при описании фильтрации жидкости в трещине. Профилометрические измерения, проведенные на стенках натурных трещин, показали, что их рельеф носит крайне нерегулярный характер и обнаруживает самоаффинные (self-affine) фрактальные свойства [33-35], причем эти свойства наблюдаются независимо от типа пород и происхождения тещин [32]. То же самое, по-видимому, можно сказать и об апертуре натурных трещин. Эти свойства сами по себе могут приводить к аномальным эффектам при переносе в потоке жидкости в отдельной трещине [23].

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

Гидравлические тесты. Один из наиболее распространенных методов испытаний — гидравлические испытания (pumping tests, pressure transient test, borehole hydraulic tests, pressure response tests), в которых наблюдается изменение напора со временем на разных расстояниях при отборе или инжекции воды в скважину.

Для интерпретации результатов испытаний необходимо знать геометрические характеристики исследуемой среды. В частности, важнейшей характеристикой является зависимость сечения, через которое проходит поток, исходящий из точечного источника, от расстояния от источника. Так, для трехмерной геометрии (неограниченная однородная пористая среда) Acer2, для двумерной геометрии (тонкий пласт, ограниченный параллельными плоскостями или система параллельных трещин) Acer1, для одномерной геометрии (канал или система каналов постоянного сечения) Аосг0.

Эксперименты в трещиноватых горных породах (fractured rocks) часто дают результаты, которые не могут быть интерпретированы на основе этих соотношений ни для какой размерности. Это радикально отличает трещиноватые среды от неконсолидированных пористых, нетрещиноватых сред [36-38].

Все указанные выше соотношения годятся по существу только для однородных (или упорядоченных неоднородных) сред. Для неоднородных неупорядоченных сред, например для фрактальной системы трещин, имеющей фрактальную размерность, отличную от 1, 2 или 3, эти соотношения не имеют места.

Yortsos [37] предложил формализм расчета переходного режима, основанный на теории фракталов, а в [39] эта теория была проверена на численной модели фрактальной сети трещин. В работе [40] было проведено формальное обобщение геометрических соотношений уравнения переноса влаги на пространство произвольной размерности (не обязательно целой). При этом предполагается, что Accrß, где 0 < ß < 2 — некоторое число, вообще говоря, не целое. Показатель ßсвязан с размерностью пространства п: ß = n-\. О величине п в [40, 38] говорится как о размерности потока (flow dimension). Такое обобщение можно ассоциировать с представлением о поровом объеме среды как о фрактальном объекте с дробной фрактальной размерностью.

В работе [41] эти идеи были применены к исследованию свойств водоносного слоя больших размеров. Исследования проводились на площадке Campus Test Site, University of the Free State, Bloemfontein, South Africa. Ha этой площадке размерами приблизительно 200x200 м имеется более 20 скважин на глубину 20-30 м. Характерной особенностью этой площадки является наличие нескольких систем трещин почти горизонтальных и наклонных. Эти трещины являются основными каналами фильтрации подземных вод.

Наиболее интересные для нас результаты были получены в гидравлических испытаниях, когда производилась откачка воды из одной из скважин и наблюдался временной ход падения напора (drawdown) на самой откачиваемой скважине и на удаленных скважинах на расстоянии 15-50 м. Ход изменения напора со временем зависит от геометрии задачи: в цилиндрическом случае (п = 2) падение происходит по логарифмическому закону со временем, в плоском (п = 1) — по корневому. Наблюдаемые в эксперименте [41] кривые не относились ни к одному из этих типов. Полученные данные были использован для оценки величины "flow dimension" п. Результаты следующие.

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

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

Аналогичные методы и идеи в последнее время широко используются разными исследователями для характеристики свойств проницаемых сред на реальных площадках [43-47]. Так, в работе [44] описываются результаты гидравлических испытаний и тестов с трассерами, причем для интерпретации этих тестов использовались понятия размерности потока. Анализ результатов полевых гидравлических тестов, проведенных на Waste Isolation Pilot Plant site, Carlsbad, привели к выводу о дробной размерности потока я, характеризующей исследуемую фильтрующую среду. Это можно считать указанием на фрактальную структуру сети проводящих трещин и каналов.

В [44] обращено внимание на то, что подобные гидравлические тесты не всегда дают определенное значение размерности потока п, характеризующее проницаемую среду. Результаты экспериментов могут и не привести к определенному значению п из-за того, что это значение меняется на разных стадиях эксперимента. В [44] приведены примеры сред, которые на основании полевых испытаний могут быть характеризованы определенной размерностью потока п, и тех, которые не могут. В последнем случае тесты разных типов, например, гидравлические испытания (pumping tests) и тесты с трассером в сходящемся потоке (radial convergent tracer test) могут давать разные результаты. Это связано с тем, что тесты разных типов охватывают разные объемы среды вокруг испытательной скважины.

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

В экспериментах, проведенных в crystalline aquifer in Britanny (France) [48], помимо оценки размерности потока в гидравлических испытаниях (которая оказалась п «1.6), исследована зависимость характерного времени отклика давления от расстояния от откачиваемой скважины. Эта зависимость оказалась степенной с неклассическим дробным показателем ~2.8.

В [46, 47] показано, что в переходные режимы отклика напора в гидравлических тестах могут проявлять свойства дробной размерности потока не обязательно только в том случае, если система каналов, проводящих воду в окрестности испытательной скважины, является фрактальной, то есть автомодельной (selfsimilar) на всех масштабах. Отсутствие свойства фрактально-сти системы каналов на всех масштабах может приводить к смене типов переходных режимов. В любом случае обнаружение переходных режимов, характеризуемых дробной размерностью служит указанием на существенно неоднородный характер сети трещин и каналов, и стохастическая адвекция в такой сети может проявлять аномальные свойства.

В [45] отмечается, что извлекаемая их данных опыта размерность потока может зависеть еще от степени нестационарности течения и граничных условий проведения эксперимента.

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

Трассерные тесты (Borehole Tracer Tests), в частности, в трещиноватых породах (fractured rocks) широко используются в последнее время для исследования переноса загрязнений в поземных водах с учетом их физико-химического взаимодействия с твердой матрицей [49-51]. Такие эксперименты дают также важную информацию о дисперсионных свойствах проницаемых сред.

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

Основная информация таких экспериментах содержится во временном ходе концентрации трассера в откачиваемой скважине (breakthrough curves) при залповом выбросе примеси через инжекционную скважину. Эти кривые содержат информацию о дисперсионных и сорбционных свойствах исследуемой среды.

Типичной особенностью трассерных экспериментов (tracer tests) в трещиноватых породах (по сравнению с неконсолидированными средами (unconsolidated porous media)) является наличие пологого и продолжительного заднего фронта. Обычно считается, что задний фронт кривой обязан диффузии трассера в стенки проводящих каналов и процессам адсорбции-десорбции на стенках. На пути к откачиваемой скважине (pumping well) часть трассера задерживается на стенках за счет упомянутых механизмов, а затем, когда концентрация примеси в основном потоке становится мала, эта часть поступает снова в поток и принимается в откачиваемой скважине в течение длительного времени.

Однако есть ряд экспериментальных данных, показывающих, что происхождение хвоста может не иметь отношения к диффузии в стенки матрицы. К ним относится работа [49], где описываются результаты трассерных тестов, проведенных в Hubbard Brook Experimental Forest, the southern part of the White Mountains, Grafton Country, New Hampshire. В этих экспериментах расстояние между инжекционной и приемной скважинами было 36 м.

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

Форма хвоста в широком диапазоне изменения концентрации (около 2 порядков величины) оказывается с хорошей точностью степенной, то есть концентрация меняется со временем по степенному закону C(t) ос гг. Показатель у, оцененный на основе опубликованных в [49] данных, оказался равным 1.7-1.8.

Эксперимент проводился в конфигурации «слабого диполя» (5% откачиваемой воды возвращалось в систегиу через инжектирующую скважину). Это течение практически эквивалентно просто сходящемуся потоку, без инжекции воды. В таких условиях если задний фронт временного профиля концентрации определяется целиком диффузией трассера в матрицу, то он должен быть степенным, но показатель у в этом случае должен быть у = 3/2. Если же ход концентрации в приемной скважине определяется только «классической» дисперсией при фильтрации, приводящей к гауссовому пространственному профилю, то оба фронта, передний и задний, должны быть экспоненциальными.

Таким образом, обнаруженный в экспериментах [49] показательный задний фронт с отличным от 3/2 показателем может служить указанием на проявление в данном эксперименте «аномальной» природы дисперсии, связанной со стохастической адвекцией во фрактальной системе трещин.

Что касается переднего фронта профиля концентрации, то он оказывается экспоненциальным. При малых (по сравнению с максимумом) концентрациях его можно представить в виде С(/)осехр{-(//г)Л|, где время t отсчитывается от максимума концентрации, а т — некоторая константа. Важно, что оценка значения Я по данным [49] дает Л * 2, как должно было быть при классической диффузии трассера.

В эксперименте [50] представлен пример другого поведения трассеров в отношении взаимодействия с твердой матрицей, когда это взаимодействие значительно влияет на ход концентрации в приемном колодце. Эксперимент проведен в Aspo Hard Rock Laborarory, Oskarshamn, Sweden. В эксперименте использовался радиально сходящийся поток, создаваемый откачкой воды из приемной скважины. Расстояние между инжектирующей и приемной скважинами было около 5 м. Основным каналом перетока между скважинами, по-видимому, являлась большая трещина, простирающаяся далеко от места расположения скважин. В [23] показано, что аномальные режимы дисперсии трассеров могут возникать и при фильтрации через одиночную трещину, еели ее апертура является сильно флуктуирующей величиной с дальнодейст-вующими корреляциями.

Для различных трассеров наблюдается различное время прихода примеси к приемнику, различные положения максимума концентрации и разные значения концентрации в максимуме. Измерение концентраций так же, как и в [49], проведено в широком динамическом диапазоне, перекрывающем 2-4 порядка величины.

Из представленных в [50] данных можно заключить, что для всех трассеров форма заднего фронта, как и в [49], является показательной. Что касается показателя у,; то он оказывается разным для разных трассеров. Для трассеров, которые значительной степени подвержены влиянию взаимодействия с матрицей, этот показатель оказывается близок к 1.5.

Среди использованных трассеров были имеющие низкую сорбционную способность и низкую диффузию в матрицу. Профили концентрации, относящиеся к этим трассерам, близко совпадают друг с другом, и можно предположить, что влияние взаимодействия с матрицей на них незначительно (что соответствует условиям эксперимента [50]). Для этих трассеров показатель степенного закона оказывается значительно больше — порядка 3.1. Таким образом, данные, опубликованные в [50], также можно интерпретировать как проявление аномальной дисперсии при адвекции в трещиноватой горной породе.

В работе [51] описан трассерный эксперимент, проведенный на ground water system surrounding a high-volume wastewater disposal well in the Florida Keys. Это эксперимент с инжекцией трассера в скважину, которая использовалась для сброса отходов, причем концентрация трассера измерялась одновременно на нескольких наблюдательных скважинах, расположенных на разных расстояниях (от 5 до 180 м) от инжекционной. Водоносный слой, в который инжектировался трассер, представляет собой известняк. Пористая среда является бимодальной. Исследование образцов матрицы показало, что наряду с микромасштабными порами имеются трещины и каналы макроскопических размеров, которые, по-видимому, обеспечивают проницаемость среды.

В эксперименте использовался «консервативный» трассер SF6, который практически не взаимодействует с твердой матрицей. Как и во многих других трассерных экспериментах в трещиноватых породах, "breakthrough curve" имеет пологий задний фронт (tailing) на всех наблюдательных скважинах (monitoring well). К сожалению, представленные в [51] данные о концентрации перекрывают сравнительно узкий диапазон изменения, что не дает возможности определить форму хвоста профиля концентрации. Однако можно сделать некоторые оценки на основании данных о времени прихода максимума концентрации и величине концентрации в максимуме на разных расстояниях.

В однородном плоском слое (размерность п = 2) время t прихода максимума концентрации на расстоянии г от источника должно зависеть от г как /осг2. На основании экспериментальных данных, приведенных в [51], можно сказать, что если эта зависимость выполняется, то только на малых расстояниях от источника — до 5 м. На больших расстояниях t растет с г медленнее. В духе представления о дробной размерности потока при фильтрации в системе трещин фрактальной структуры можно предположить, что tazrs. Показатель S трудно оценить численно на основании экспериментальных данных [51] в силу значительного разброса значений t для разных скважин, расположенных на близких расстояниях от источника, но можно определенно утверждать, что S заметно меньше 2.

Для консервативного трассера (то есть не поглощаемого матрицей) максимальная концентрация примеси Ст на расстоянии г в однородном плоском слое при «классической» дисперсии должна быть Cmocl//15. На основании экспериментальных данных [51] показатель степенной зависимости между Ст и t существенно меньше классического значения 1.5.

Все это дает основания предполагать, что здесь мы имеем дело с проявлением аномальной дисперсии, связанной со стохастической адвекцией в трещиноватых горных породах.

Вообще, неклассическое, нефиковское поведение передних и задних фронтов (early and late time tailing) профилей концентрации примеси (breakthrough curves) является типичным для трещиноватых сред. Негауссовы хвосты концентрационных профилей наблюдались в экспериментах [52-54]. Для формального описания таких результатов в [55, 56] выдвинута гипотеза о масштабозависимой дисперсии (scale-dependant dispersion). Предполагается, что движение облака примеси (plume) описывается перемещением центра тяжести и диффузионным расширением. При этом коэффициент дисперсии, по аналогии с турбулентной диффузией, является функцией размеров облака примеси. Рассеяние примеси и увеличение размеров облака происходят благодаря флуктуациям скорости фильтрации масштаба порядка или меньше размеров облака. Поэтому для осуществления супердиффузионного режима необходимо, чтобы имелись дальнодействующие корреляции поля скоростей.

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

Трассерные тесты в естественных условиях без искусственных воздействий (natural gradient tracer tests). В экспериментах этого типа исследуется распространение трассера, инжектированного за короткое время в водоносный слой. Течение в слое после инжекции примеси происходит под действием естественного градиента и не подвергается никаким искусственным воздействиям (откачка, инжекция). Концентрация трассера измеряется с помощью сети наблюдательных скважин. Эксперименты этого типа в принципе дают возможность прямого наблюдения изменений размеров облака примеси со временем, а также пространственного профиля концентрации.

Имеются подробные описания двух серий данных трассерных тестов, проведенных с достаточно хорошим пространственным и временным разрешением, что позволяет детально проследить профиль облака (plume) и его эволюцию. Один из них проведен на площадке Cape Cod, Massachusetts [57]. Другой — macrodispersion experiment (MADE) — был проведен на базе Columbus air force, Mississippi [58, 59].

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

В экспериментах в водоносный слой инжектировались около 10 м3 воды, содержащей одновременно разные трассеры, включая «консервативные», такие как бромид, тритий (тритированная вода). Наблюдения продолжались от 15 до 20 месяцев в разных экспериментах.

Наблюдения распределений концентрации примесей, проводившиеся в течение почти года, дали следующие результаты. Расширение облака происходит в горизонтальном направлении, то есть двумерным образом, вертикальное расширение мало. Размеры облака увеличиваются со временем быстрее, чем по диффузионному закону (r2}oct2/a. При этом продольный (вдоль направления среднего течения) и поперечный размеры характеризуются разными значениями показателя or. а «1.2 для продольного размера, и а «1.5 для поперечного (or = 2 соответствует обычному диффузионному закону). С этим согласуются данные об изменении со временем максимальной и средней концентрации в плюме (plume).

Кроме того, обнаружено, что профиль концентрации трития (эксперимент MADE-2) существенно несимметричен в направлении среднего распространения: он имеет сильно вытянутый в сторону течения хвост, причем этот хвост является не экспоненциальным, а степенным.

Качественно похожие в отношении «супердиффузионного» характера расширения облака результаты были получены и на площадке Cape Cod. Водоносный слой в этом случае также был песчано-гравийного состава, однако степень его неоднородности, по-видимому, была несколько ниже, чем на площадке MADE. Так же, как и в MADE, вертикальная дисперсия плюма оказалась мала, так что диффузия была по существу двумерной.

Было установлено, что расширение облака в продольном к среднему потоку и поперечном направлениях происходит по-разному. Продольное расширение имеет супер диффузионный характер с показателем а «1.65. Поперечное рассеивание плюма описывается обычным диффузионным законом (а «2). Подробный анализ данных этих экспериментов можно найти в [68, 69].

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

Зона неполного насыщения (Vadose zone). Основной вопрос, касающийся оценки проводящих свойств ненасыщенных трещиноватых пород — это вопрос о механизме распространения влаги в таких средах. Согласно классической капиллярной модели [60,61], вода за счет капиллярного эффекта впитывается в твердую матрицу и распространяется за счет фильтрации по ней. При этом трещины являются препятствиями для движения воды на большие расстояния. Для сред, наиболее интересных с точки зрения задачи захоронения отходов (например, трещиноватый туф Yucca Mountain), по-видимому, более реалистична другая модель, согласно которой основным каналом распространения воды являются именно трещины, а капиллярное впитывание представляет сравнительно слабый эффект [62, 63]. В этом случае фильтрационный поток оказывается крайне неоднородным и нестационарным, наблюдаются наличие предпочтительных путей распространения [64]. Режим распространения воды и переноса примесей определяется статистическими свойствами сети трещин.

К настоящему времени проведено немного экспериментальных исследований фильтрации воды и транспорта примесей в трещиноватых ненасыщенных породах. В [65, 66] описана серия полевых экспериментов с измерениями структуры фильтрационного потока и транспорта примесей через трещиноватую породу, а также внутренней структуры трещин. Наблюдения [65,66] показали пространственную и временную нестабильность потока, сильное каналирование, когда большая часть потока (70-100%) проходит через малую часть доступного сечения трещин (15-20%). При этом «активные» пути движения воды постоянно меняются в зависимости от циклов смачивания/осушения, химического взаимодействия потока со стенками, отложений на стенках растворенных в воде материалов. Все это значительно усложняет картину фильтрации в ненасыщенной зоне по сравнению с насыщенной.

В [67] описаны результаты пневматических испытаний, проведенных в туннеле на Yucca Mountain, Nevada, USA, который является предположительным местом захоронения отходов. В блоке с размерами 15x20x15 м было просверлено в разных направлениях около 30 скважин длиной 5-10 м. В каждой скважине имелся уплотненный участок, через который подавался воздух с постоянным расходом. Одновременно с подачей воздуха измерялось давление в самой нагнетаемой скважине и во всех остальных. Процедура нагнетания и измерения давления повторялась последовательно со всеми скважинами. Отклик давления быстро, в течение нескольких минут устанавливался на стационарном значении. Эти данные позволили оценить проницаемость вблизи каждой скважины (усредненную по длине уплотненной области), а также с помощью моделирования тестов (программа TOUGH2) путем подбора проницаемости во всей расчетной области. В результате получена трехмерная карта проницаемости на сетке 34x26x24.

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

Таким образом, приведенные результаты лабораторных и полевых экспериментов показывают, что для сильно неоднородных неупорядоченных сред, какими являются трещиноватые горные породы, неклассический, нефи-ковский характер переноса представляется скорее правилом, чем исключением. При этом рост облака частиц примеси со временем (/?(/)) может идти быстрее, чем R~yft, а убывание концентрации при г » R(t) — медленнее, чем по гауссову закону. Если «хвосты» убывают по степенному закону, то их называют «тяжелыми». Причиной наличия аномальных режимов переноса и «тяжелых хвостов» в данном случае является сложная неоднородная структура натуральных сетей трещин, которые, как показывают непосредственные наблюдения, проявляют фрактальные свойства и попадают под определение перколяционных сред. Отсюда следует, что предельно допустимые концентрации радионуклидов могут распространяться на расстояния, которые на много порядков превышают найденные из классических представлений. Все это означает, что для получения адекватных оценок надежности радиоактивных захоронений, требуется разработка новых подходов.

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

Если дисперсия приращений координат становится бесконечной, то в одномерном случае средняя плотность описывается уравнением дробной диффузии с оператором дробного пространственного дифференцирования, определяемым двумя параметрами, характеризующими функцию распределения. Это параметр а, задающий степенное убывание «хвостов», и параметр р, задающий степень асимметрии распределения. При а = 2 распределение стремится к нормальному, при а = 1 — переходит в распределение Коши. При а < 1 функция распределения приращений координат блуждающей частицы характеризуется не только бесконечной дисперсией, но и бесконечно большой средней величиной. И хотя, в этом случае, мы не можем получить устойчивых средних для каждой из частиц в отдельности, средние значения плотности частиц в заданной точке пространства остаются устойчивыми переменными и имеют привычный физический смысл.

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

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

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

Модель стохастической адвекции представляет собой яркий пример того, как во фрактальной, сильно неупорядоченной среде может реализоваться режим дробной диффузии (супердиффузии) на умеренных расстояниях с одновременным отсутствием «тяжелых» степенных «хвостов» в распределении концентрации на больших расстояниях.

Целью настоящей работы является разработка вычислительных методов решения уравнений диффузии с дробной производной по пространству и по времени в одномерном случае, разработка алгоритма для проведения прямого численного моделирования стохастической адвекции примеси в сильно неоднородных геологических средах, создание и отладка соответствующих программ. Результаты проделанной работы используются [72-73, 78-81, 83, 121-122] для отработки методик решения обратных задач распространения радионуклидов в сильно неоднородных средах для оценки безопасности подземных захоронений и хранилищ высокоактивных ядерных отходов.

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

Заключение

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

Разработаны вычислительные методы и создан комплекс программ для решения уравнения диффузии с дробной производной по времени в одномерном случае. За основу предложенных алгоритмов приняты определения Римана-Лиувилля и Капуто. Показана аппроксимация и устойчивость полученных схем. Проведено сравнение результатов численных расчетов в случае Капуто с фундаментальными решениями. Показано отсутствие «тяжелых» степенных «хвостов» в полученных решениях.

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

1. Ю.А. Коровин, В.М. Мурогов, «Экологически приемлемый и безопасный топливный цикл ядерной энергетики». Обнинск, 1991.

2. Ю.А. Коровин и др., «Современные проблемы ядерной энергетики».

3. Г.В. Киселев, «Проблемы развития ядерной энергетики». Москва, 1990.

4. A.M. Бахметьев и др., «Ядерная энергетика будущего: атомная энергетика России на пороге XXI века». Москва, 1995.

5. В.А. Красноселов, «Введение в атомную энергетику». Ульяновск, 2004.

6. А.В. Кацай, «Атомная энергетика в вопросах и ответах». Москва, 2003.

7. National Research Council. Rock fractures and fluid flow: contemporary understanding and applications. Washington, DC: National Academy Press, 1996.

8. M. Kuntz, P. Lavalee. Experimental Evidence and Theoretical Analysis of Anomalous Diffusion during Water Infiltration in Porous Building Materials. J. Phys. D: Applied Phys. 34 (2001), 2547-2554.

9. M. Kuntz, P. Lavalee. Anomalous Spreading of a Density Front from an Infinite Continuous Source in a Concentration-Dependent Lattice Gas Automaton Diffusion Model. J. Phys. D: Applied Phys. 36 (2003), 1135-1142.

10. M. Kuntz, P. Lavalee. Anomalous diffusion is a rule in concentrated-dependent diffusion process. J. Phys. D: Applied Phys. 37 (2004), L5-L8.

11. U. M. Scheven and P.N. Sen. Spatial and temporal coarse graining for dispersion in randomly packed spheres. Phys. Rev. Lett. V. 89, No 25, 2002.

12. S.P.K. Sternberg. Dispersion Measurements in Highly Heterogeneous Laboratory Scale Porous Media. Transport in Porous Media. 54, 107-124, 2004.

13. Raven K.G., Gale J.E. Water flow in natural rock fracture as a function of stress in and sample size. Int. J. Rock. Mech. Min. Sci. Geomech. Abstr. 1985,22(4), 251-261.

14. Pyrac-Notle L.J., Meyer L.R., Cook N.G.W., Witherspoon P.A. Hydraulic and mechanic properties of natural fractures in low permeable rock. Proc. 6th Int. Cong. Mech. 1987, 225-231.

15. Durham W.B., Bonner B.P. Self-propping andfluidflow in slightly offset joint at high effective pressures. J. Geophys. Res. 1994, 99(B5), 9391-9399.

16. Keller A. A., Roberts P.V., Kitanidis P.K. Prediction of single-phase transport parameters in a variable aperture fracture. Geophys. Res. Lett. 1995, 22(11), 1425-1428.

17. Vandergraaf T.T. Radionuclide migration experiments under laboratory conditions. Geophys. Res. Lett. 1995, 22(11), 1409-1412.

18. Nolte D.D., Cook N.G.W., Pyrac-Notle LJ. The fractal geometry of flow paths in natural fractures in rock and the approach to percolation. Pure Appl. Geophys. 1985, 131, 111-138.

19. Nicholl M.J., Glass R.J., Nguen H.A. Small-scale behavior of single gravity-driven fingers in an initially dry fracture. In: Proc. High Level Radioactive Waste Manag. Conf. Las Vegas, NV, 1993.

20. Rasmussen T.C. Steady fluid flow and travel times in partially saturated fractures using a discrete air-water interface. Water Resour. Res. 1991, 27, 67-76.

21. Rasmuson A., Neretnieks I. Radionuclide transport in fast channels in crystalline rock. Water. Resour. Res. 1986, 22, 1247-1256.

22. G.W., Su J.T. Geller, K. Prues, F. Wen. Experimental studies of water seepage and intermittent flow in unsaturated, rough-walled fractures. Water Re-sour. Res. V. 35, No 4. p. 1019-1037, 1999.

23. S. Roux, F. Plouraboue, J.P. Hulin. Tracer Dispersion in Rough Open Cracks. Transport in Porous Media. 32, 87-116, 1998.

24. Karasaki K., Friefeld B., Cohen A., Grossenbacher K., Cook P., Vasco D.A. Multidisciplinary fractured rock characterization study at Raymond field site, Raymond, SA. Journal of Hydrology, 236 (2000), 17-34.

25. Odling N.E. Scaling and connectivity ofjoint systems in sandstone from western Norway. J. Struct. Geol. 1997, 19, 563-571.

26. Bonnet E., Bour O., Odling N.E., Davy P., Main I., Cowie P., Bercowitz B. Scaling of fracture systems in geological media. Reviews of Geophys., 2001, 39(3), 347-383.

27. Isichenko M.B. Percolation, statistical topography and transport in random media. Rev. Mod. Phys., 1992, V. 64, No 4, p. 961-1043.

28. Stauffer D., Aharony A. Introduction to percolation Theory. London: Taylor & Francis, 1992.

29. Bercowitz B. Analysis of fracture network connectivity using percolation theory. Math. Geol., 1995, V. 27, 467-484.

30. Gueguen Y., David C., Gavrilenko P. Percolation and fluid transport in the crust. Geophis. Res. Lett., 1999, 26, 2001-2004.

31. K. Ando, A. Kostner, S. Neuman. Stochastic continuum modeling of flow and transport in crystalline rock mass: Fanay-Augeres, France, revisited. Hydro-geology Journal, (2003) 11, 525-531.

32. Berkowitz B. Characterizing flow and transport in fractured geological media: A review. Advances in Water Resources, 25 (2002), 861-884.

33. Brown S.R., Scholz C.H. Broad bandwidth study of the topology of natural rock surfaces. J. Geophys. Res. 1985, 90(B24), 1275-1282.

34. Poon C.Y., Sayles R.S., Jones. T.A. Surface measurements and fractal characterization of naturally fractured rocks. J. Phys. D, 1992, 25(8), 1269-1275.

35. Schmittbuhl J., Schmitt F., Sholz C.H. Scaling invariance of crack surfaces. J. Geophys. Res. 1995, 100(B4), 5953-5973.

36. Novakowski K.S., Bikerton G.S. Borehole measurement of the hydraulic properties of low-permeability rock. Water Resour. Res. 1997, 33(11), 25092517.

37. Chang J., Yortsos Y.S. Pressure Transient Analyses of Fractal Reservoirs. In: Paper SPE 18170. Presented on 63rd Annual SPE Technical Conference and Exhibition, Soc. of Pet. Eng., Houston, Tex., October 2-5, 1988.

38. K. Riemann, G. Van Tonder. Interpretation of single-well tracer tests using fractional-flow dimensions. Part 1: Theory and mathematical models. Hydro-geology Journal (2002) 10, 351-356.

39. Acuna J.A., Yortsos Y.S. Applications of Fractal Geometry to the Study of Network of Fractures and their Pressure Transient. Water Resour. Res., 1995,31(3), 527-540.

40. Barker J.A. A generalized radial flow model for hydraulic tests in fractured rock. Water Resources Research, Vol. 24, No 10, 1988, p. 1796-1804.

41. K. Riemann, G. Van Tonder. Interpretation of single-well tracer tests using fractional-flow dimensions. Part 2: A case study. Hydrogeology Journal (2002) 10, 357-367.

42. J.M. Filho. A fractal dimension model for hydraulic parameter evaluation from single well tests in fractured media. Proceedings of Int. Conf. On Groundwater in Fractured Rocks (2003), Prague, p. 251-252.

43. R.L. Beauheim, R.M. Roberts, J.D. Avis. Well testing in fractured media: Flow dimensions and REVs. Proceedings of the International Groundwater Symposium, 2002, p. 277-280. IAHR, Paseo Bajo Virgen del Puero 3, 28005, Madrid, Spain.

44. D.D. Walker, R.M. Roberts. Summary of flow Dimensions Corresponding to Hydrologic Conditions. Proceedings of the International Groundwater Symposium, 2002, p. 295. IAHR, Paseo Bajo Virgen del Puero 3, 28005, Madrid, Spain.

45. H. Jourde, S. Piste, P. Perrochet, C. Drogue. Origin of Fractional Flow Dimension to a Partially Penetrating Well in Stratified Fractured Reservoir. New Results Based on Synthetic Fracture Network. Advances in Water Resources 25 (2002), 371-387.

46. Doe T.W. Fractional dimension analyses of constant pressure well tests. In Proceedings, Annual SPE Technical Conference and Exhibition, Formation Evaluation and Reservoir Geology, Soc. of Pet. Eng, Dallas, Tex, October 6-9, 1991, p. 461-467.

47. Le Borgne T., Bour O., De Dreuzy J-R, Davy P. Flow model relevant to crystalline aquifers: insights from a scaling interpretation of pumping.

48. M.W. Becker, A.M. Shapiro. Tracer transport in fractured crystalline rock: Evidence of non-diffusion breakthrough tailing. Water Resources Research, Vol. 36, No. 7, p. 1677-1686.

49. H. Widestrand, P. Anderson, J. Byegard, G. Skarnemark, M. Skalberg, E. Wass. In-situ Migration Experiments at Aspo Hard Rock Laboratory, Sweden: Results of Radioactive Tracer Migration Studies in a Single Fracture.

50. Journal of Radioanalytical and Nuclear Chemistry, V. 250, No. 3 (2001), 501-517.

51. K. Dillon, W. Burnett, G. Kim, J. Chanton, D.R. Corbett, K. Elliot, L. Kump. Groundwater Flow and Phosphate Dynamics Surrounding a High Discharge Waste Water Disposal Well in the Florida Keys. Journal of Hydrology, 284 (2003), 193-210.

52. Silde R.C., Nilsson B., Fredericia J. Spatially varying hydraulic and solute transport characteristics offractured till determined by field tracer tests, Fu-nen, Denmark. Water. Resour. Res. 1998, 34, 2515-2527.

53. McKay L.D., Sanford W.E., Strong J.M. Field-scale migration of colloidal tracers in a fractured shale saprolite. Ground Water, 2000, 38(1), 139-147.

54. Kosakowski G., Bercowitz B., Scher H. Analysis offield observation of tracer transport in fractured till. J. Contam. Hydrol. 2001, 47(1), 29-51.

55. Gelhar L.W., Wetly C., Rehfeldt K.R. A critical review of data on field-scale dispersion in aquifers. Water. Resour. Res. 1992, 28(7), 1955-1974.

56. H. Rajaram, L.W. Gelhar. Plume-scale dependent dispersion in aquifers with a wide range of scales of heterogeneity. Water Resour. Res. V. 31, No 10, 1995, p. 2469-2482.

57. Adams E.E., Gellar L.W. Field study of dispersion in heterogeneous aquifer, 2. Spatial moments analysis. Water Resour. Res. 1992, 28(12), 3293-3307.

58. Boggs J.M., Beard L.M., Long S.E., McGee M.P. Database of the second macrodispersion experiment (MADE-2). EPRI report TR-102072, electric Power Res. Inst., Palo Alto, Ca., 1993.

59. Wang J.S.Y., Narasimhan T.N. Hydrologie mechanisms governing fluid flow in a partially saturated, fractured porous medium. Water Resour. Res. 1985, 21, 1861-1874.

60. Peters R.R., Klavetter E.A. A continuum model for water movement in unsaturatedfracture rock. Water Resour. Res. 1988, 26, 415-424.

61. Pruess K., Faybichenko B., Bodwarsson G. Alternative concepts and approaches for modeling flow and transport in thick unsaturated zones of fractured rocks. J. Contam. Hydrol. 1999, 38, 281-322.

62. Pruess K. A mechanistic model for water seepage through thick unsaturated zones of fractured rocks. Water Resour. Res. 1999, 35(4), 1039-1051.

63. Foster S.S.D., Carington S.A. The infiltration of tritium in the chalk unsaturated zone. J. Hydrol., 1980, 46, 343-364.

64. Dahan O., Nativ R., Adar E.M., Berkowitz B., Ronen Z. Field observation of flow in a fracture intersecting unsaturated chalk. Water Resour. Res. 1999, 35(11), 3315-3326.

65. Dahan O., Nativ R., Adar E.M., Berkowitz B., Weisbrod N. On fracture structure and preferential flow in unsaturated chalk. Ground Water, 2000, 3893,444-451.

66. Huang K., Tsang Y.W., Bodvarsson G.S. Simultaneous inversion of air-injection tests in fractured unsaturated tuff at Yucca Mountain. Water Resour. Res. 1999, 35(8), 2375-2386.

67. Meerschaert M., Benson D.A., Baumer B. Operator Levy motion and multis-caling anomalous diffusion. Phys. Rev. V. 63, 2001, 021112.

68. Benson D.A., Shumer R., Meerschaert M., Wheatcraft S.V. Fractional dispersion, Levy motion, and the MADE tracer tests. Transport in porous media, 42, 2001,211-240.

69. В.М. Головизнин, В.П. Киселев, И.А. Короткин, Ю.И. Юрков, «Некоторые особенности вычислительных алгоритмов для уравнения дробной диффузии». Препринт №1ВЯАЕ-2002-01, Москва: ИБРАЭ РАН, 2002. 57 с.

70. В.М. Головизнин, В.П.Киселев, И.А. Короткин, «Численные методы решения уравнения дробной диффузии в одномерном случае». Препринт №1ВЯАЕ-2002-10, Москва: ИБРАЭ РАН, 2002. 35 с.

71. В.М. Головизнин, В.П. Киселев, И.А. Короткин, Ю.И. Юрков, «Решение обратной задачи по идентификации основных параметров дробной диффузии». Препринт №1ВЯАЕ-2002-20, Москва: ИБРАЭ РАН, 2002. 20 с.

72. В.М. Головизнин, В.П. Киселев, И.А. Короткин, «Численные методы решения уравнения диффузии с дробной производной по времени в одномерном случае». Препринт №1В11АЕ-2003-12, Москва: ИБРАЭ РАН, 2003. 35 с.

73. В.М. Головизнин, В.П. Киселев, В.И. Питербарг, И.А. Короткин, Ю.И. Юрков, «Генерация случайных полей трещиноватости с заданной корреляционной функцией и заданным средним». Препринт №1В11АЕ-2003-17, Москва: ИБРАЭ РАН, 2003. 23 с.

74. В.М. Головизнин, В.П. Киселев, И.А. Короткин, Ю.И. Юрков, «Решение задач диффузионного переноса на случайных полях трещиноватости идробная диффузия». Препринт №IBRAE-2004-01, Москва: ИБРАЭ РАН, 2004. 28 с.

75. JI.A. Большов, В.М. Головизнин, A.M. Дыхне, В.П.Киселев, П.С. Кондратенко, В.Н. Семенов, «Новые подходы к оценке безопасности захоронений радиоактивных отходов». Журнал «Известия РАН. Энергетика». 2004. №4. С. 99-108.

76. В.М. Головизнин, В.П. Киселев, И.А. Короткин, Ю.И. Юрков, «Прямые задачи неклассического переноса радионуклидов в геологических формациях». Журнал «Известия РАН. Энергетика». 2004. №4. С. 121-130.

77. В.М. Головизнин, В.П. Киселев, В.Н. Семенов, И.А. Короткин, Ю.И. Юрков, «Решение обратной задачи неклассического переноса радионуклидов в геологических формациях с использованием нейронных сетей». Журнал «Известия РАН. Энергетика». 2004. №5. С. 81-87.

78. V.M. Goloviznin, P.S. Kondratenko, L.V. Matweev, V.N. Semenov, I.A. Korotkin, "Direct Numerical Modeling of Stochastic Radionuclide Ad-vection in Highly Non-Uniform Media". Preprint №IBRAE-2005-01. Moscow: Nuclear Safety Institute RAS, 2005. 37 p.

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

80. A.M. Нахушев. «Элементы дробного исчисления и их применение». Нальчик, издательство КБНЦ РАН, 2000. 299 с.

81. Ross В. (Ed.). Fractional Calculus and its Applications, lecture notes in mathematics. Springer-Verlag, Berlin, 1975.

82. Saichev A.I., Zaslavsky G.M. Fractional kinetic equations: Solutions and applications. J. Chaos 7, p. 753-764, 1997.

83. A.M. Нахушев. «Уравнения математической биологии». М.: Высш. шк., 1995.301 с.

84. Динариев О.Ю. «Фильтрация в трещиноватой среде с фрактальной геометрией трещин». Механика жидкости и газа. 1990. №5. С. 66-70.

85. Кочубей А.Н. «Диффузия дробного порядка». Дифференциальные уравнения. 1990. Т. 26, №4. С. 660-670.

86. Mainardi F. Fractals and Fractional Calculus in Continuum Mechanics. Eds. Carpinteri A. and Mainardi F., CISM Courses and Lectures, Springer-Verlag, Wie,p. 291-348, 1997.

87. Montroll E.W., Weiss G.H. Random walks on lattices. J. Mathematical Phys. 6, p. 167-181, 1965.

88. Shlesinger M., Klafter J., Wong Y.M. Random walks with infinite spatial and temporal moments. J. Statist. Phys. 27, p. 499-512, 1982.

89. Benson D. The Fractional Advection-Dispersion Equation: Development and Application. A dissertation submitted in partial fulfillment of the Doctor of Philosophy in Hydrogeology, University of Nevada, Reno, 1998.

90. Schumer R., Benson D., Meerschaert M., Wheatcraft S. Eulerian derivation of the fractional advection-dispersion equation. J. Contaminant Hydrology, 38(1/2), p. 69-88, 2001.

91. Lu S., Molz F.J., Liu H.H. An efficient, three-dimensional, anisotropic, fractional Brownian motion and truncated fractional Levy motion simulation algorithm based on successive random additions. Computers & Geosciences, 29, p. 15-25,2003.

92. Lutz E. Fractional transport equations for Levy stable processes. Physical Review Letters, vol. 86, Issue 11,2001, pp. 2208-2211.

93. F.R. Gantmacher. Theory of Matrix, 1973, Moscow (In Russian).

94. W. Chen. A new definition of the fractional Laplacian. Simula Research Laboratory, P. O. Box. 134, NO-1325 Lysaker, Norway, 2002. Eprint arXiv:cs/0209020.

95. F. Mainardi, P. Paradisi and R. Gorenflo. Probability Distributions Generated by Fractional Diffusion Equations. TeX Pre-print, Bologna, January 1998, edited for J. Kertesz and I. Kondor (Editors).

96. R. Gorenflo and F. Mainardi. Fractional Calculus and Stable Probability Distributions. Arch. Mechanics, 50 (3), 377-388, 1998.

97. L. Vazques. Fractional Diffusion Equations with Internal Degrees of Freedom. Journal of Computational Mathematics (2002, in press).

98. M. Caputo. Linear Models of Dissipation whose Q is Almost Frequency Independent, Part II, Geophys. J. R. Astr. Soc., 529-539, 1967.

99. W. Chen. Positive Time-Fractional Derivative. Simula Research Laboratory, P. O. Box. 134, 1325 Lysaker, Norway, 30 September 2002. Eprint arXiv:cs/0210005.

100. Dranikov I.L., Kondratenko P.S., Matweev L.V. Anomalous Contaminant Transport in Stochastic Advection Model. Preprint IBRAE-2003-20. Moscow: Nuclear Safety Institute RAS, December 2003. 23 p.

101. A.M. Дыхне, П.С. Кондратенко, Л.В. Матвеев. «Перенос примеси в пер-коляционных средах». Письма в ЖЭТФ, том 80, вып. 6, с. 464-467.

102. Hundsdorfer W., Koren В. et al. A Positive Finite-Difference Advection Scheme. Journal of Computational Physics, 1995, 117, pp. 35-46.

103. Карабасов С.А. «Разностные схемы с пространственным расщеплением временной производной для задач двухфазной фильтрации». Диссертация на соискание ученой степени кандидата физ.-мат. наук. ИБРАЭ РАН. Москва. 1999.

104. Головизнин В.М., Самарский А.А. «Разностная аппроксимация конвективного переноса с пространственным расщеплением временной производной». Мат. моделирование. 1998. Т. 10. №1. С. 86-100.

105. Головизнин В.М., Самарский А.А. «Некоторые свойства разностной схемы «КАБАРЕ». Мат. моделирование. 1998. Т.10. №1. С. 101-116.

106. Головизнин В.М., Карабасов С.А. «Нелинейная коррекция схемы «КАБАРЕ». Препринт ИБРАЭ №IBRAE-98-08. Москва: ИБРАЭ РАН, 1998. 17 с.

107. Головизнин В.М., Карабасов С.А., Кобринский И.М. «Балансно-характеристические схемы с разделенными консервативными и потоковыми переменными». Мат. моделирование. 2003. Т. 15. №9. С. 29-48.

108. V.M. Goloviznin, V.N. Semenov, I.A. Korotkin, S.A. Karabasov, "A Novel Computational Method for Modelling Stochastic Advection in Heterogeneous Media ", Transport in Porous Media, 2006 (accepted).