Трехмерное моделирование процессов переноса примесей в пористых средах сложной структуры

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

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

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

Вышеизложенные свойства коэффициентов и расчетных сеток сужают класс возможных дискретизаций. Так, применение хорошо разработанных методов конечно-разностных дискретизаций [16, 19] для неструктурированных тетраэдральных сеток затруднительно. Наиболее перспективными аппроксимациями для рассматриваемого класса задач являются, па взгляд автора, методы конечных элементов (МКЭ)[15],[18] и методы конечных объемов (МКО)[35], в рамках которых будут рассмотрены некоторые подходы к построению вычислительных схем.

Разработка численных методов для уравнения конвекции-диффузии традиционно отталкивается от методов пространственной аппроксимации диффузионного оператора. Самым простым и популярным является МКЭ с пространством непрерывных кусочно-линейных функций [18], Pi-МКЭ. Недостатками этого метода являются немонотонность на произвольных неструктурированных сетках и отсутствие аппроксимации дискретных потоков на границах ячеек сетки. Несмотря на то, что возможно восстановление дискретных потоков на границах ячеек вспомогательной сетки (например, сетки Вороного), эти ячейки могут быть невыпуклы, содержать много граней и включать разрывы коэффициента диффузии. Стремление построить метод, приближающий диффузионные потоки на границах ячеек заданной сетки, привело к появлению других аппроксимаций. Метод смешанных конечных элементов (МСКЭ) [30] одновременно аппроксимирует как само решение, так и диффузионные потоки. Решение может иметь разрывы на границах элементов, диффузионные потоки находятся в пространстве Равьяра-Тома [56], благодаря чему обеспечивается непрерывность их нормальной составляющей на общих гранях двух элементов и, следовательно, локальная консервативность на ячейках сетки. К недостаткам МСКЭ в приложении к диффузионным задачам относится немонотонность схемы на неструктурированных и анизотропных сетках [28], [41]. Метод Галеркина с разрывными элементами (БС, 018соп{-тюи8 Са1ег1ап) [32], [25] является современной альтернативой вышеперечисленным методам в тех случаях, когда необходима гибкость при задании степени аппроксимирующего полинома для каждой ячейки сетки или используемые сетки являются неконформпыми. По своим свойствам метод близок к МСКЭ [27].

Метод конечных объемов (МКО) широко применяется в инженерных вычислениях. Существует множество вариаций МКО: методы для разных типов сеточных ячеекметоды со степенями свободы в вершинах сетки или в ячейкахметоды с двухточечным или многоточечным шаблоном, аппроксимирующим поток через грань ячейкиметоды с линейной или нелинейной аппроксимацией потока, и т. д. В применении к неструктурированным триан-гуляциям для МКО с двухточечной аппроксимацией диффузионного потока накладываются существенные ограничения на сетки (ортогональность сеток или сетки Делоне и Вороного [55],[37]), а для МКО с многоточечными шаблонами для потоков требуются сетки с регулярными ячейками [21], [34]. Кроме того, наличие полного тензора диффузии существенно усложняет задачу построения линейной схемы МКО. В диссертационной работе предложен новый монотонный нелинейный МКО для произвольных сеток и тензоров диффузии (или проницаемости для фильтрационных задач) [11]. Метод использует идеи, опубликованные в работах [52],[53] для двумерных задач диффузии. Под монотонностью метода здесь и далее подразумевается свойство неотрицательности получаемого численного решения при соответствующих граничных условиях и правой части, обеспечиваемое монотонностью матрицы ([3], стр.269) аппроксимации диффузионного оператора, аналогично работам [52], [53]. В случае однородных краевых условий типа Дирихле такая трактовка понятия монотонности соответствует принципу максимума с диагональным преобладанием по столбцам в матрице аппроксимации, сформулированному А. А. Самарским и П. Н. Вабищевичем в [16](стр.63). Неотрицательность численного решения приобретает исключительную важность в задачах, где требуется учет химических взаимодействий, что отмечено в работе [36]. Наличие отрицательной концентрации в расчетной области может приводить к неверному расчету химических реакций и неконсервативности модели.

Переход от задачи диффузии к задаче конвекции-диффузии сопряжен с дополнительными трудностями при построении устойчивых схем в случае доминирующей конвекции. Во-первых, для стабилизации схем могут применяться методы регуляризации [16]. Классическим примером является Pi-МКЭ в сочетании с алгоритмом SUPG (Streamline Upwinding Petrov Galerkin) [42]. В случае МСКЭ и МКО используется противопотоковая аппроксимация конвективного члена. Во-вторых, для аппроксимации конвективного оператора возможно использование схем сквозного счета, основанных на методах Годунова [4], Лакса и Фридрихса [49] и др. Примерами таких схем являются методы, использующие разрывные конечные элементы (РКЭ) [58], [33], а также метод коррекции потока [44], обеспечивающий монотонность и минимальную численную диффузию, который, однако, не применим к случаю полных и неоднородных тензоров диффузии. Эти методы позволяют аппроксимировать негладкие решения с малой численной диффузией и без осцилляций.

Для аппроксимации по времени нестационарных задач конвекции-диффузии могут использоваться явные схемы, неявные схемы и методы расщепления [14]. При этом явные и неявные схемы используют одну и ту же пространственную аппроксимацию конвективных и диффузионных членов (МКЭ, МСКЭ, МКОсравнение методов на двумерной тестовой задаче проведено в [28]), а методы расщепления допускают построение разных пространственных аппроксимаций отдельно для каждого из физических операторов. Так, например, МСКЭ может использоваться для аппроксимации диффузионных потоков, РКЭ — для аппроксимации решения [58]. Доусон [33], Аккерер и др. [23] предлагают схемы расщепления с явной аппроксимацией конвекции и неявной — диффузии. Вохралик [61] использует МКО для конвективного члена и иеконформные КЭ или МСКЭ для диффузионного оператора в сочетании с полностью неявной схемой. Схемы расщепления позволяют для каждого из пространственных операторов выбирать наиболее адекватный метод аппроксимации, поэтому именно они рассматриваются в работе как наиболее приспособленные для моделирования процессов переноса примесей в пористых средах.

Рассматриваемый класс трехмерных проблем относится к разряду больших задач, требующих больших вычислительных мощностей, машинной памяти, расчеты занимают дни и недели. Для повышения эффективности вычислений применяются различные технологии: в работах [33], [23] предлагаемые схемы расщепления допускают применение разных шагов по времени для конвекции и диффузии, а также локальных шагов по времени в подобластяхВ.И.Агошковым и др. [24] предложен метод декомпозиции области для процессов конвекции-диффузии в неоднородных средах, который позволяет разбить задачу на несколько подзадач в подобластях, решая только задачу конвекции в подобластях с сильным преобладанием коивекции над диффузией. В настоящей работе предложен ряд методов ускорения вычислений. Па-раллелизация алгоритма [2],[10] дает возможность решать задачи на сетках, содержащих несколько миллионов тетраэдров. Технология автоматического динамического перестроения (локального сгущения и разгрубления) сеток сокращает время работы программы благодаря меньшему числу неизвестных и обеспечивает высокую точность в областях быстрого изменения решения [9]. Двухуровневый переобуславливатель для задач диффузии в пластовых средах [46],[47], разработанный и исследованный в сотрудничестве с Ю. А. Кузнецовым, убирает зависимость скорости сходимости метода сопряженных градиентов от коэффициентов диффузии и шага сетки по вертикальной оси.

Работы с результатами численного моделирования реальных объектов захоронения загрязнений редко встречаются в печати. В тематическом выпуске журнала Computational Geosciences (т.8, № 2, 2004) опубликованы результаты расчетов задания Couplex, предложенного французским агентством по ядерным отходам. Задание состояло из трех задач: трехмерного моделирования мелкомасштабных процессов конвекции-диффузии с учетом химических реакций в малой окрестности захоронения ядерных отходовдвумерного моделирования крупномасштабного переноса загрязнения от источникасовмещения мелкомасштабной модели источника загрязнений с крупномасштабной моделью двумерного переноса. П. И. Чаловым и др. в статье [20] представлены результаты расчета переноса загрязнений от хвостохрапилища карабал-тинского гидрометаллургического завода по переработке урапосодержащих руд (Кыргызстан). В этой задаче поле конвективных потоков было двумерным, для дискретизации использовались гексаэдральпые сетки. В настоящей работе проведено численное моделирование новой полностью трехмерной задачи крупномасштабного переноса загрязнения па неструктурированных тетраэдральных сетках. Задание представляет большой интерес с методической точки зрения, т.к. расчет затрудняют следующие факторы: анизотропия (2 порядка) и неоднородность (до 10 порядков) тензора проницаемостисильная анизотропия области, и, как следствие, расчетной сетки (диаметр области — 55 км, глубина — 500м) — неоднородность и анизотропия полного тензора диффузиивыклинивания геологических пластовнеобходимость прогноза на большой период времени (сотни тысяч лет).

Основные результаты диссертационной работы:

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

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

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

4. Проведено численное моделирование трехмерной задачи с реальными параметрами.

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

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

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

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

Автор выражает глубокую признательность научному руководителю Ю. В. Василевскому за постановку задачи и полезные обсужденияЮ.А.Кузнецову за предоставленную возможность участвовать в разработке и исследовании нового переобуславливателяД.А. Святскому, В. Н. Чугунову, и К. Н. Липникову за плодотворные дискуссии по возникавшим в ходе написания диссертации вопросамколлективу ИВМ РАН за создание благоприятной рабочей атмосферы. Такэюе автор благодарен первому научному руководителю И. В. Новоэюилову, пробудившему интерес к научной работе.

Выводы по главе 4.

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

Решение задачи переноса также осложнялось сжатием сетки по вертикали, анизотропией и неоднородностью тензора диффузии, разными временными масштабами процессов в верхних и нижних геологических пластах. В результате расчетов получен прогноз времени и места выхода на земную поверхность загрязненных вод с концентрациями радиоактивного Йода 129, равными Ю-14, Ю-12, Ю~10кг/м3. Показана динамика распространения загрязнения в трехмерной области с течением времени, вплоть до 350 000 лет.

Заключение

.

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

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

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