Рассматриваются аспекты гидрогеологического моделирования и разработки соответствующих расчетных кодов для обоснования безопасности захоронения радиоактивных отходов (РАО). Анализ текущих разработок позволяет выделить три основных тенденции: возрастание сложности моделей; коды «рожденные параллельными»; переход на неструктурированные сетки и соответствующие им методы дискретизации. Эти концепции демонстрируются на примере расчетного кода GeRa (Геомиграция Радионуклидов) с использованием ряда сложных моделей: переноса с геохимическими взаимодействиями, насыщенно-ненасыщенной фильтрации, плотностной конвекции.
Гидрогеологическое моделирование в задачах обоснования
безопасности захоронений радиоактивных отходов с
помощью расчетного кода GeRa*
Разработка кода GeRa также ориентирована на использование адаптивных многогранных
неструктурированных сеток, параллельность, что обеспечивается разработкой кода на
программной платформе MSPP (близка к платформе INMOST [
К настоящему времени в коде GeRa реализованы возможности геологического моделирования, численного моделирования процессов фильтрации и переноса, оценки и визуализации результатов. В совокупности с дополнительным модулем расчета дозовых нагрузок реализован расчет доз для населения, получаемых при различных сценариях водопользования. Код снабжен графическим интерфейсом, обеспечивающим препроцессинг (задание геометрии расчетной области, начальных и граничных условий, параметров моделей, скважин, поверхностных водоемов и т.д.), запуск расчетов, а также двумерную и трехмерную визуализацию полученных результатов. На рисунке 1 показан интерфейс программы с двумерным (слева) и трехмерным (справа) окнами редактирования и визуализации.
Рис. 1. Интерфейс расчетного кода GeRa Геологическое моделирование включает построение структурной модели объекта и полей фильтрационных и миграционных параметров геологических сред. Для построения поверхностей кровель и подошв слоев, либо полей параметров могут быть использованы методы интерполяции (триангуляция, обратные расстояния) и геостатистики (кригинг, дополненный вариограммным анализом).
Для моделирования геофильтрационных и геомиграционных процессов в настоящее время
реализованы следующие модели:
фильтрация в напорной, напорно-безнапорной и насыщенно-ненасыщенной постановке.
адвективно-диффузионно-дисперсионный массоперенос.
радиоактивный распад.
сорбция по линейной изотерме с постоянным коэффициентом распределения Kd, либо с
переменным Kd, который зависит от концентрации одной из переносимых компонент. В
последнем случае Kdзадается таблицей значений.
плотностная конвекция.
перенос с полным расчетом химических взаимодействий. Непосредственно расчет
химических реакций выполняется с помощью свободно распространяемого кода
PHREEQC [
При моделировании учитываются различные гидрогеологические объекты: поверхностные водоемы, скважины, дренажные системы. 3. Особенности дискретизации задач
Методы дискретизации задач ориентированы на использование произвольных конформных (в обобщенном смысле) многогранных трехмерных сеток. В коде GeRa реализовано два сеточных генератора: генератор треугольно-призматических сеток с исключением вырожденных ячеек и возможностью образования тетраэдральных и пирамидальных ячеек при вырождении одной или двух образующих призм; генератор гексаэдральных сеток на основе восьмидеревьев с методикой «скалывания» ячеек [8].
Оба генератора позволяют аппроксимировать плановые границы расчетной области и поверхности кровель и подошв слоев со вторым порядком точности. Треугольно-призматический генератор обеспечивает двумерную плановую адаптацию сетки, в то время как гексаэдральный – трехмерную адаптацию сетки к особенностям задач.
Для дискретизации оператора диффузии по времени может быть использован один из двух
линейных методов конечных объемов (МКО): традиционная двухточечная схема
аппроксимации потоков и более точная, но вычислительно трудоемкая, многоточечная О-схема
аппроксимации потоков [
Для моделирования конвективно-диффузионного переноса могут быть использованы следующие численные схемы: явно-неявная схема расщепления по физическим процессам (явная конвекция, неявная диффузия), для адвекции в ней может быть использована низкодиссипативная схема высокого разрешения (MUSCL-типа) с кусочно-линейным восстановлением концентрации на ячейках сетки, либо стандартная схема МКО с кусочно-постоянной аппроксимацией концентрации; полностью неявная схема нелинейного монотонного МКО (с MUSCLсхемой для конвективного члена); полностью неявная схема с кусочно-постоянной аппроксимацией концентрации, использующая один из линейных МКО для дискретизации диффузионного члена.
При моделировании нелинейных задач, в частности, процессов безнапорной и насыщенноненасыщенной фильтрации, используется метод простой итерации. Для дискретизации по времени задач плотностной конвекции, которые неразделимы на подзадачи фильтрации и переноса, используется итерационное сопряжение моделей этих процессов, при котором на каждом шаге по времени цикл расчета отдельных подшагов выполняется до достижения заданного критерия сходимости. В силу высокой вычислительной сложности решения задач переноса с химическими взаимодействиями, в них выполняется безытерационное сопряжение, при котором реакционный и миграционный подшаги делаются лишь один раз на каждом шаге по времени. 4. Примеры применения кода
В настоящем разделе приводится ряд примеров использования кода GeRa для выполнения
расчетов в рамках достаточно сложных моделей: плотностной конвекции, безнапорной и
насыщенно-ненасыщенной фильтрации, переноса с химическими взаимодействиями.
4.1 Моделирование напорно-безнапорной и ненасыщенной фильтрации, переноса
в зоне аэрации
Данная задача используется в качестве обучаеющей к расчетному коду VS2DT[11], ее
постановка и численное решение можно найти в работе [
Рис. 2. Геометрия расчетной области и граничные условия Таблица 1. Параметры к двумерной задаче фильтрации и массопереноса
Параметр Коэффициент фильтрации при полном насыщении (м/сут) Влагосодержание при полном насыщении s Остаточное влагосодержание r Параметр уравнения состояния n Параметр уравнения состояния (м-1) Продольная дисперсивность, L (м) Поперченная дисперсивность, T (м) константа распада, (сут-1) Песок 12,4 0,39 0,06 5,8 2,3 0,1 0,01 0,03 0,1 0,03 Рис. 3. Распределение давления в расчетной области при T=110 суток, посчитанное моделью насыщенноненасыщенной фильтрации (слева) и напорно-безнапорной фильтрации (справа) Рис. 4. Концентрация загрязнения в расчетной области при T=110 суток, посчитанная при моделировании фильтрации насыщенно-ненасыщенной моделью (слева) и напорно-безнапорной моделью (справа) Моделирование проводится с помощью двух моделей: напорно-безнапорной и насыщенноненасыщенной фильтрации, для демонстрации различий между ними. Как видно из рисунка 3, картина распределения давления в расчетной области кардинально отличается в этих двух моделях. Более простая модель напорно-безнапорной фильтрации не способна учесть сложного нелинейного поведения гидрофизических параметров породы, что ставит под сомнение ее использование при моделировании фильтрации и переноса в зоне аэрации, хотя такие модели являются компонентами расчетных кодов MODFLOW и FEFLOW.
При переносе загрязнения (см. рисунок 4) модели также показывают принципиально разные результаты. Во-первых, в напорно-безнапорной постановке загрязнение распространяется значительно быстрее. Во-вторых, линза слабопроницаемых пылеватых грунтов выполняет различные роли. В модели напорно-безнапорной фильтрации она служит водоупором, и поток обтекает ее. В насыщенно-ненасыщенной постановке, напротив, она служит коллектором потока, а на ее нижней границе образуется капиллярный барьер, отводящий поток вправо, вдоль нижней границы линзы. 4.2 Моделирование плотностной конвекции
Одной из наиболее известных задач, используемых при верификации численных моделей
плотностной конвекции, является задача Элдера, предложенная Элдером в [
Стоит отметить, что в коде GeRa возможен адекватный учет чисто примесных источников
(когда поступление в систему примеси не сопровождается поступлением жидкости), без
которого расчет напоров и потоков в ряде задач (в частности, в задаче Элдера) неверен, как это
показано в работе [
Рис. 5. Решение задачи Элдера на сетке из 44×25 ячеек Рис. 6. Решение задачи Элдера на сетке из 400×100 ячеек 4.3 Моделирование переноса с химическими взаимодействиями жины, фильтры которых находятся на высоте от -2 м до 2 м. Расход каждой из этих скважин составляет 10 м3/сутки. В центре расчетной области, точке (0;0) расположена разгрузочная скважина, фильтр которой также находится на высоте от -2 м до 2 м. Расход разгрузочной скважины компенсирует закачку и составляет 40 м3/сутки. Порода пласта и раствор, заполняющий ее в начальный момент времени, находятся в химическом равновесии. На породе присутствуют две ионообменных фазы.
В нагнетательные скважины осуществляется закачка растворов, содержащих нитрат натрия и стронций. Концентрация стронция в закачиваемых растворах одинаковая и составляет 3,3х10-9 моль/л. Концентрация нитрата натрия варьируется и составляет: 0 моль/л в скважине P1; 0,2 моль/л в скважине P2; 0,4 моль/л в скважине P3; 0,6 моль/л в скважине P4.
Моделирование проводится на период в 3000 дней. На рисунках 7 и 8 показаны концентрации нитрата натрия и стронция в сечении слоя плоскостью Z=0 м в моменты времени 1000 дней и 3000 дней. Отмечены изолинии концентрации стронция 10-10 моль/л и 10-9 моль/л, а также изолинии концентрации нитрата натрия от 0,1 до 0,6 моль/л с шагом 0,1 моль/л. Во-первых, данная модель позволяет учесть влияние различного содержания нитрата натрия на миграцию стронция: чем выше концентрация нитрата, тем меньше сорбируется стронций. Во-вторых, в отличие от других моделей, наблюдается замедление миграции самого нитрата, то есть, он уже не является нейтральным трассером, как традиционно принято считать. а) T=1000 суток б) T=3000 суток Рис. 7. Распределение нитрата натрия в расчетной области при учете химических взаимодействий а) T=1000 суток б) T=3000 суток Рис. 8. Распределение стронция в расчетной области при учете химических взаимодействий 5. Заключение
В работе кратко описаны основные возможности интегрального кода GeRa для гидрогеологического моделирования в задачах обоснования безопасности захоронений РАО. В соответствии с современными тенденциями разработки кодов подобной направленности, GeRa является кодом, «рожденным параллельным» благодаря использованию программной платформы MSPP, и работает на адаптивных трехмерных неструктурированных сетках.
В расчетном коде реализованы базовые и усложненные модели фильтрации и переноса, включая модели напорно-безнапорной фильтрации, насыщенно-ненасыщенной фильтрации, плотностной конвекции, переноса в зоне аэрации, переноса с химическими взаимодействиями. Код отличает наличие ряда уникальных моделей, не включенных в современные широко распространенные коды (семейство MODFLOW, FEFLOW и др.), в частности, модель плотностной конвекции позволяет адекватно учитывать чисто примесные источники, а при моделировании сорбции по изотерме возможен учет влияния концентрации одной из переносимых компонент на коэффициент распределения других веществ.
Расчетный код GeRa позволяет провести полный цикл моделирования при решении задачи оценки безопасности захоронения РАО: геологическое моделирование, геофильтрационное и геомиграционное моделирование, оценка и визуализация результатов. К настоящему времени первая версия кода прошла опытную эксплуатацию во ФГУП «Национальный оператор по обращению с радиоактивными отходами» и готовится к подаче на аттестацию в Ростехнадзор. Литература 1. Капырин И.В., Уткин С.С., Василевский Ю.В. Концепция разработки и использования расчетного комплекса GeRa для обоснования безопасности пунктов захоронения радиоактивных отходов // Вестник атомной науки и техники, серия «Математическое моделирование физических процессов», – 2014 – №4 – C.44-54. 11. Официальная страница расчетного кода VS2D
http://wwwbrr.cr.usgs.gov/projects/GW_Unsat/vs2di1.3/index.html . Hydrogeological modeling in radioactive waste disposal safety assessment using the GeRa code Ivan Kapyrin, Igor Konshin, Kopytov German, Kirill Nikitin and Yuri Vassilevski Keywords: groundwater flow, transport in porous media, reactive transport, density-driven flow, disposal safety assessment The challenges of hydrogeological modeling in the framework of radioactive waste disposal safety assessment are discussed. Analysis of the ongoing projects shows three major trends in code development: increasing models complexity; codes "born parallel"; adaptive polyhedral grids and suitable high-order numerical schemes. These concepts are demonstrated using the GeRa (Geomigration of Radionuclides) numerical code in terms of a couple of advanced models: reactive transport, unsaturated flow and transport in vadose zone, density-driven flow.