Моделирование тепловой конвекции в пористых средах с учетом объемного тепловыделения в коде GeRa

ЧЕБЫШЕВСКИЙ СБОРНИК Том 18 Выпуск 3

УДК 621.039 DOI 10.22405/2226-8383-2017-18-3-234-253

МОДЕЛИРОВАНИЕ ТЕПЛОВОЙ

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

GERA

Ф. В. Григорьев1, И. В. Капырин2, Ю. В. Василевский3 (г.

Москва)

Аннотация

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

В первом разделе представлен краткий обзор известных гидрогеологических расчетных кодов, обладающих функционалом для решения тепловых задач (FEFLOW, SUTRA, SEAWAT, TOUGH2). Отмечается отсутствие возможности учета объемного тепловыделения вследствие радиоактивного распада в данных программных комплексах. Дается описание расчетного кода GeRa собственной разработки.

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

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

1 Григорьев Федор Владимирович, инженер-исследователь ИБРАЭ РАН, grig-fedor@yandex.ru

2 Капырин Иван Викторович, заведующий лабораторией геомиграционного моделирования Института проблем безопасного развития атомной энергетики РАН, старший научный сотрудник Института вычислительной математики РАН, ivan.kapyrin@gmail.com

3 Василевский Юрий Викторович, член-корреспондент Российской академии наук, заместитель директора ИВМ РАН, yuri.vassilevski@gmail.com

стью неявного итерационного сопряжения, на каждой итерации которой последовательно решаются задачи фильтрации, тепло- и массопе-реноса.

В четвертом разделе приводится тестовая задача о конвекции тепловыделяющей жидкости в замкнутой двумерной полости с изотермическими стенками, заполненной пористой средой. Результаты решения кодом ОеЫа сравниваются с асимптотическим решением, полученным Наа^гаёеЬ. В пятом разделе представлены результаты моделирования в коде ОеЫа экспериментов Вигейа п Ветшал по изучению режимов свободной тепловой конвекции жидкости с объемным тепловыделением в пористой среде. Сравниваются зависимости числа Нуссельта от числа Рэлея, полученные в натурном эксперименте и численно.

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

Ключевые слова: численное моделирование, тепловая конвекция в пористых средах, объемное

Библиография: 24 названий.

MODELING OF THERMAL CONVECTION IN POROUS MEDIA WITH VOLUMETRIC HEAT

SOURCE USING THE GERA CODE

F. V. Grigoriev, I. V. Kapyrin, Y. V. Vassilevski

Abstract

This article is devoted to the problem of thermal convection in porous media with volumetric heat generation modelling, arising in practice of radioactive waste (RW) disposal safety assessment.

In the first section a brief overview of widespread hvdrogeological codes (FEFLOW, SUTRA, SEAWAT, TOUGH2) featuring'the ability to solve thermal problems is done. We point out the lack of heat generation caused by radioactive decay model in these programs. The GeRa numerical code developed by the authors is presented.

In the second section we consider the mathematical model of coupled groundwater flow, solute and heat transport, which is implemented in GeRa. The model describes these processes in saturated porous media and takes into account radioactive decay, sorption on the rock, the dependences of density and viscosity on temperature. The heat transport equation is written assuming thermal equilibrium between the fluid and the rock. The model includes heat transport by convection and conduction-thermal

dispersion. The heat source terms can be wells and volumetric heat generation due to radioactive decay.

The numerical scheme implemented in GeRa to solve the aforementioned coupled problem is introduced in the third section. The space discretization is done using finite volume methods (FVM). Sequential iterative coupling implicit scheme is used for temporal discretization. On each iteration of the scheme the flow, heat transport and solute transport problems are solved sequentially.

The fourth section is devoted to the test problem of heat generating fluid convection in a closed two-dimensional cavern filled by porous material with isothermal walls. The results obtained using GeRa code are compared to the asymptotical solution deduced by Haajizadeh.

In the fifth section we present the results of modelling with GeRa the experiments of Buretta and Berman in which they investigated the regimes of free thermal convection of fluid with volumetric heat generation in porous media. The dependences of Nusselt number on the Ralev number measured in the experiments and calculated numerically are compared.

In the sixth section we consider the test problem of continuous injection of high-level RW into an aquifer. Here the ability to model coupled flow, heat and solute transport processes is shown. The numerical solution obtained using GsRxt is compared to a known analytical one.

Keywords: numerical modeling, thermal convection in porous media, the GeRa code, volumetric heat source.

Bibliography: 24 titles.

1. Введение

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

Модели тепловых процессов существуют во многих гидрогеологических расчетных кодах. Один из наиболее известных - FEFLOW (Finite Element Subsurface Flow and Transport Simulation System). Разработка данной программы началась в 1979 году Диршем (Hans-Jorg G. Dierseh) и продолжается до сих пор подразделением MIKE организации DHI. Код способен моделировать трехмерную фильтрацию в пористых средах с учетом плотноетных эффектов, тепло-, массопереноса, химических взаимодействий. Для дискретизации уравнений используются конечно-элементные методы. Подробная документация возможностей кода в области физико-химических моделей и численных методов приведена в книге [1].

SUTRA - конечно-элементный код, выпущенный геологической службой США (USGS) в 1984 году [2] для двумерных задач и в 2002 году [3] для трехмерных задач. Последнее обновление программы вышло в 2010 году. С помощью этого кода было исследовано множество задач, связанных со свободной конвекцией, ставших затем бенчмаркинговыми [4, 5]. Моделирование может проводиться только на гексаэдральных сетках. Одновременное моделирование тепло- и массоперепоса невозможно.

SEAWAT - еще один код, выпущенный USGS. Дискретизация уравнений основана на методе конечных разностей. Важной особенностью программы можно назвать то, что она является результатом синтеза двух популярных гидрогеологических пакетов: MODFLOW (моделирование фильтрации) и MT3DMS (моделирование массопереноса) [6]. Последнее обновление SEAWAT v I вышло в 2012 году.

TOUGH2 (Transport of Unsaturated Groundwater and Heat) [7] - численный код, разрабатываемый национальной лабораторией им. Лоуренса в Беркли для моделирования геотермальных процессов и задач, связанных с изоляцией высокоактивных РАО. Пространственная дискретизация в программе основана напрямую на законах сохранения в интегральной форме без перевода их в уравнения в частных производных (т. п. "integral finite difference method" - IFDM). Используя различные модули программы, можно учесть двойную пористость и проницаемость среды, многофазность течения (в том числе фазовые переходы).

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

В данной работе излагается модель теплоперепоса и тепловой конвекции, реализованная в коде GeRa. Расчетный код GeRa разрабатывается специалистами ИБРАЭ РАН и IIBM РАН. Он предназначен для гидрогеологического моделирования в задачах по обоснованию безопасности захоронений РАО (см. [8]). Моделирование геофильтрации, геомиграции и теплоперепоса в GeRa осуществляется на неструктурированных многогранных сетках с помощью метода конечных объемов. Для выполнения операций с сетками и их элементами, для сборки и решения линейных систем, а также для хранения данных в последовательном и параллельном режимах в коде используется программная платформа INMOST [9, 10]. Способность кода GeRa моделировать процессы теплопереноса и тепловой конвекции с учетом объемного тепловыделения вследствие радиоактивного распада в данной работе проверяется сравнением численного решения верификационных задач с референтными результатами, основанными на экспериментальных данных и на аналитических расчетах.

2. Модель теплопереноса в пористой среде

Система уравнений, описывающая сопряженную модель фильтрации, массо-, теплопереноса в насыщенной пористой среде с учетом радиоактивного распада и сорбции, выглядит следующим образом [11]:

ма

дк дТ Ксотр

рв— - РРо3~^ + к'ио1'1 ~+ У(р^) = РзЯз

дЪ

и = -К[ У к + У г

Ро

г=1

0 У^

дъ

,дт

[<ррос^ +

д

т

+ лрИгСг + У

Выпишем для нее замыкающие соотношения:

К

р = Ро(1 -Р(Т - То)) +

сотр £

=1

к = —

К

+ г, о

крод

р

X = + (1 - р)Хг + ррос^Б^р, О = Ddiff + Ddisр, ^^р = \\и\\((а1 - аг)Е(и) + агI).

Здесь использованы следующие обозначения: р - плотность жидкости, Б -коэффициент упругой емкости, к - гидравлический напор, <р - пористость среды, р0 - плотность жидкости без примесей при референтной температуре, Т - температура среды, 3 = -- коэффициент объемного теплового расширения жидкости, - концентрация г-й примеси, = - коэф-

фициент объемного расширения для г-й примеси, и = (и1,и2,и3) - скорость фильтрации, р3 - плотность жидкости в источнике, qs - интенсивность источников и стоков, К - тензор фильтрации, г - вертикальная координата, отсчитываемая относительно некоторого уровня, Cf - удельная теплоемкость жидкости, рг - плотность породы, сг - удельная теплоемкость породы, Л -теплопроводность пористой среды (жидкости и породы), Xf - теплопроводность жидкости, Аг - теплопроводность породы, Т3 - температура источника, Ш - мощность объемного тепловыделения, Я^ - коэффициент сорбциоппой задержки для г-го компонента, Лг - константа радиоактивного распада для г-го компонента, С3,г - концентрация в источнике, р - давление жидкости, д

- ускорение свободного падения, & - проницаемость среды, ^ - динамическая вязкость жидкости, D¿iff - изотропный тензор эффективной молекулярной диффузии, Б ¿пар - тензор дисперсии, од - продольная дисперси вность, од -поперечная дисперсивпость, Е(и) = щ^-, I ~ единичная матрица размером 3x3. "

Уравнение (3) записано в предположении теплового равновесия между раствором и породой. Сорбция в уравнении (4) также предполагается равновесной.

Формула для тепловыделения вследствие радиоактивного распада с учетом собпрованного вещества выглядит следующим образом:

Ш =

г г

Здесь Мл - число Авогадро, Мг - молярная масса г-го вещества, Е^ - средняя энергия, выделяющаяся при одном акте распада, - доля энергии, уносимая нейтрино (для случая бета-распада).

Формула (6) записана в предположении о том, что тепловыделение происходит в той же точке, где произошел распад. Это предположение обосновано тем фактом, что длина ионизационного пробега в воде альфа-частиц (с типичными для радионуклидов РАО энергиями в диапазоне от 3 до 8 МэВ) составляет 30-100 мкм. Для бета-частиц - 1.5-5 мм. Длина половинного ослабления гамма-излучения с энергиями 0.5-4 МэВ в воде составляет 28-35 см [12]. Такие расстояния гораздо меньше характерных пространственных масштабов фильтрационных процессов, протекающих в геологических средах.

3. Численный метод

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

р (Кп+1 ,Сп+1,тп+1) = 0,

Н

М (Кп+1 ,Сп+1) = 0.

Здесь Кп+\ тп+\ сп+1 - векторы значений напора, температуры, концентраций растворенных компонентов в ячейках расчетной сетки на следующем, (п + 1)-ом шаге по времени. Система (7) состоит из трех подсистем нелинейных уравнений. Первая, Р(Кп+1, Сп+1 ,Тп+1) = 0 соответствует уравнению фильтрации (1). Вторая, Н(Кп+1,Сп+1,Тп+1) = 0, состоит из дискретных уравнений теплопереноса, соответствующих (3) (зависимость от Сп+1 обусловлена источником радиогенного тепла). Третья подсистема соответствует уравнению массопереноса (4).

При решении транспортных подзадач системы (7) могут быть использованы явно-неявная и полностью неявная дискретизации по времени, В первом случае для дискретизации оператора конвекции применяется явная ТУБ-схема метода конечных объемов (МКО) с кусочно-линейным восполнением концентрации на ячейках расчетной сетки (аналогично [13]), Во втором случае дискретизация конвективного оператора осуществляется с помощью неявной противопотоковой схемы первого порядка. Для оператора диффузии в обоих случаях используется неявная схема МКО с линейной двухточечной аппроксимацией потоков. Эта же схема применяется при решении фильтрационной подзадачи.

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

1, В качестве начального приближения выбираются значения с предыдущего шага по времени: Кп+1'° = Кп, Тп+1'° = Тп, Сп+1'° = Сп. Устанавливается счетчик итераций к = 0,

2, Рассчитывается фильтрационная задача с плотностью, концентрацией и температурой с предыдущей итерации: Р(Кп+1'к+1, сп+1,к,г^п+1,к) = 0_ Находятся напор Кп+1'к+1 и соответствующий ему фильтрационный поток иа+1'к+1.

3, Рассчитывается задача массоперепоса с известным потоком г!а+1'к+1\ М(Кп+1, Сп+1) = 0, Отсюда находится концентрация на новой итерации

£п+1,к+1

4, Рассчитывается задача теплопереноса с известным потоком цп+1,к+1 и концентрацией Сп+1'к+1\ Н(Ка+1 ,Сп+1 ,Тп+1) = 0, Отсюда находится температура на новой итерации тп+1'к+1.

5, Обновляется счетчик итераций к = к + 1,

6, Проверяется критерий сходимости итерационного процесса. Если критерий выполнен, положим Кп+1 = Кп+1'к, Тп+1 = Тп+1'к, Сп+1 = Сп+1'к. В ином случае возвращаемя к шагу 2,

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

\\1гп+1'к _ кп+1

к-11

ш.&х(\\Ьп+1'к\\, \\кп+1,к-1\\\

\\£<га+1'& _ ^п+1'к-1\\

тзх(\\Сп+1'к\\, \\Сп+1,к-1\

\\^га+1'& _ рп+1'к-1\\

[тах(\Тп+1'к\\, \\Тп+1'к

1

< ,

< ,

< &ге1.

4. Сравнение с аналитическим решением

В работе Наа]1гас1е]1 и др. [14] было выведено асимптотическое (для относительно малых чисел Рэлея Яа < 100) решение для свободной конвекции тепловыделяющей жидкости в замкнутой двумерной полости с изотермическими вертикальными стенками, заполненной пористой средой. Это решение предлагается в качестве тестовой задачи для кода Се11а.

Задача формулируется в двумерной прямоугольной области размерами Ах1 (А -аспектный параметр: отношение высоты к ширине). На вертикальных стенках задана постоянная температура Т = 0. На верхней и нижней границах задан нулевой поток тепла. На всей границе задано условие непротекания.

Рис. 1: Область решения задачи о свободной конвекции тепловыделяющей жидкости в прямоугольной каверне

Начальное условие задачи теплопереноса - однородное распределение нулевой температуры. Начальное условие фильтрационной задачи - однородное распределение напора к = 0. Расчет проводится до установления стационарного состояния.

В работе [14] показано, что основным параметром, определяющим режим течения в данной системе, является число Рэлея следующего вида (форма записи несколько отличается от [14] ввиду использования в данной работе других обозначений):

На =

где А - ширина прямоугольной каверны.

Приводится асимптотическое решение стационарной системы уравнений фильтрации-теплопереноса:

У(ри) = 0,

рос/Ъ

р = ро(1 - Р(Т - То)).

Оно представлено в виде разложения температуры и функции тока в ряд по степеням Яа:

те

Ф = ^2, RаПфп(х,У),

Г

Т = ^ ЯапТп(х,у).

п=0

Вычислены коэффициенты от п = 0 до п = 2,

Задача решается кодом (1еНа на равномерной гексаэдральной сетке с разбиением 25 х 1 х 50, Шаг по времени равен 10 сек. При решении подзадачи теплопереноса в (7) используется явно-неявная схема (см, раздел 3), Расчет производится до установления стационарного состояния. Числовые параметры задачи приведены в таблице 1,

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

Параметр Значение Система единиц

К 1 м / сек

в 10-6 м

V 0.2 -

Ро 1000 кг / м3

А/ 0.65 Вт/

А, 2.7125 Вт/

Р 10-4 1/°С

То 0 °с

с/ 4200 Дж/

1.259 Вт / м3

На рисунке 2 представлены профили обезразмеренной горизонтальной скорости вдоль прямой х = 0.2 для аспектпых отношений А = 5 А = 2 и числа Рэлея Яа = 50, полученные численно в коде СеКа в сравнении с аналитическим решением. Безразмерная компонента скорости и и размерная компонента и' связаны следующим об разом: и = u/dp0Cf /X. Рисунок 2 демонстрирует хорошее соответствие результатов аналитических расчетов [14] и численного решения кодом СеКа,

Рис, 2: Профиль горизонтальной скорости и вдоль прямой х = 0.2: численное решение кодом СеИа и аналитика |14| для различных аепектпых отношений А и числа Рэлея Яа = 50

5. Сравнение с экспериментом

В работе |15| приведены результаты экспериментов Buretta и Berman по изучению режимов свободной тепловой конвекции жидкости с объемным тепловыделением за счет джоулева тепла в пористой среде в цилиндрическом сосуде с изотермической верхней крышкой и теплоизолированной нижней крышкой и боковыми стенками. Дана теоретическая оценка критического числа Рэлея (такого числа Рэлея, при котором происходит переход из коп-дуктивного режима теплопереноса в конвективный): Rac = 32.8, Эта оценка была сверена с результатами эксперимента: Ra = 31.8, Также была построена зависимость определенного в эксперименте числа Нуссельта от числа Рэлея. Данная зависимость имела интересный вид: она разделялась па две ветки с разными наклонами в логарифмическом масштабе, причем при некотором числе Рэлея происходи,:: перескок с одной ветки па другую, В последующих экспериментах на аналогичных установках Hardee и Xilson |16|, Rhee и др. |17|, Knlacki и Freeman |18| так и не удалось воспроизвести этот скачок. Причем если в первых двух работах результаты лучше коррелировали с верхней ветвью зависимости, то в работе Knlacki и Freeman, наоборот, результаты лучше соотносились с нижней ветвью. Данные расхождения результатов до сих нор не получили окончательного убедительного объяснения |19|. Полученные в эксперименте критические числа Рэлея также значительно разнятся: Rac = 32 у Hardee и Nilson [16], Rac = 46 у Rhee и др. [17], Rac = 36 у Knlacki и Freeman |18|.

Приведем результаты моделирования в коде GeRa эксперимента Buretta

и Вегтап |15|, Рас четная область имеет форму цилиндра диаметром 30 см и высотой 4,64 см. Используется треугольно-призматическая сетка с числом слоев по вертикали 10, длина ребра ячейки у границы с1 ~ 1 см. Всего сетка содержит 28600 ячеек. Область моделирования изображена па рисунке 3,

Рис, 3: Область моделирования эксперимента по конвекции тепловыделяющей жидкости

Моделирование проводилось при числах Рэлея Ra = 15-245, Число Рэлея в данной задаче определяется следующим образом:

Ra = тЁШ±,

2А

Теплоперепое характеризуется числом Нуееельта, записываемым в следующей форме:

Wh2

N" = 2АДГ <13>

где AT = Тропот — Тир. Температура на верхней границе поддерживается постоянной Тир = 21°С, Температура на теплоизолированной нижней границе может быть неоднородной. Поэтому под величиной Tbottom подразумевается среднее значение:

Tbottam = 1 J J т

bottom

Числовые параметры задачи приведены в таблице 2,

Величина шага по времени для Ra = 15 — 105: At = 100 с, для Ra = 140 — 245: At = 25 с. При решении подзадачи теплопереноса в (7) используется полностью неявная схема (см, раздел 3), Результаты вычислительных экспериментов при различных числах Рэлея представлены па рисунке 4,

Таблица 2: Параметры моделирования эксперимента но конвекции тепловыделяющей жидкости

Параметр Значение Система единиц

К 30 см/с

V 0.4 -

р 0.0005 1/°С

с/ 4200 Дж/

0 Дж/

А/ 0.0065 Дж/

А, 0.013 Дж/

Ро 0.001 кг/см3

в 10-6 1/см

Т ± ир 21 °С

То 21 °С

W (На = 15 — 245) 0.00052 — 0.00842 Вт/см3

Рис. 4: Зависимость числа Нуссельта от числа Рэлея: результаты численного моделирования, аппроксимирующая прямая

В районе На = 30 конвективные потоки значительно изменяют картину ноля температуры, и число Нуссельта становится больше единицы. При На = 35 — 100 точки зависимости числа Рэлея от числа Нуссельта хорошо легли в логарифмическом масштабе па аппроксимирующую прямую (см. рисунок 4), соответствующую уравнению:

^ Ми = 0.5799 ^ На — 0.8617,

что неплохо соотносится с верхней веткой зависимости, полученной ВнгсТТа и Вегтап |15|:

lg Nu = 0.553 lg Ra — 0.871.

При Ra > 105 не удалоеь получить стационарное решение задачи. Соответствующие точки на рисунке 4 определены с использованием усредненной по времени температуры Tbottam и отклоняются вниз относительно зависимости (15), На пересечении прямой (15) с осью находим критическое число Рэлея Ra™m = 30.6.

6. Задача о непрерывной закачке высокоактивных РАО в пласт-коллектор

В статье [20] приведено аналитическое решение задачи о температурном поле при закачке тепловыделяющих РАО в пласт-коллектор с учетом отведения тепла кондукцией во вмещающие породы. В работе [21] это решение было использовано для верификации программы ТОГ(П 12-КОЯ7Н. модифицированной для учета радиогенного тепла.

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

При решении задачи кодом (1сНа область разбивается на ячейки треугольно призматической сеткой (см. рисунок 5). Пласт-коллектор содержит один слой ячеек, кровля - 14 слоев. В центре области находится скважина радиусом 0.1 м. От центра области ячейки увеличиваются в размере, и у границы длина ребра достигает <1 ~ 90 м. Всего сетка состоит из 47010 призм.

Расчет проводится на период 20 лет. Шаг по времени - 50 еут. Задачи маеео-, теплопереноеа считаются с помощью полностью неявной схемы (см. раздел 3). Результаты расчетов представлены на рисунке 6 в сравнении с аналитическим решением [20]. Видно хорошее соответствие результатов. Небольшое размазывание фронта можно объяснить численной диффузией и тем, что при выводе аналитического решения пренебрегли теплопроводностью в радиальном направлении.

Таблица 3: Параметры задачи о непрерывной закачке высокоактивных РАО в пласт-коллектор

Параметр Значение Система единиц

К (пласт-коллектор) 0.8 м/еут

К (вмещающие породы) 10-6 м/еут

V 0.2 -

Ро 1000 кг / м3

РЪи1к 1920 кг / м3

в 10-6 1/м

А/ 56160 Дж/

А, 234360 Дж/

с/ 4200 Дж/

1000 Дж/

То 0 °С

ТУд (темп, закачиваемой воды) 12 °С

Период полураспада 10585 сут

Ш 2.01 Вт / м3

Е (тепловыделение на 1 акт распада 1.126 МэВ/распад

Мощность пласта-коллектора 12.5 м

Мощность кровли 87.5 м

Радиус области 700 м

Дебит скважины 164.16 м3/сут

Концентрация РАО 5.58  1013 Бк/м3

в закачиваемой жидкости

Рис. 5: Область моделирования задачи о непрерывной закачке РАО

О 100 200 300 400 500 600 700

Рис, 6: Зависимость температуры в пласте-коллекторе от расстояния до скважины: (а) - через 5 .нет, (б) - через 20 .нет. Приведены результаты численного расчета кодом СеИа в сравнении с аналитическим решением |20|

7. Заключение

В данной работе изложена численная модель теилоиереиоеа и тепловой конвекции с учетом тепловыделения при радиоактивном распаде, реализованная в расчетном коде СеИа. Собрана верификационная база, включающая аналитические тесты и экспериментальные результаты. Численное решение верификационных задач показало хорошее совпадение с референтными результатами, свидетельствуя о том, что код СеИа адекватно моделирует сопряженные процессы геофильтрации, транспорта радионуклидов, теи-лоиереиоса и тепловыделения при радиоактивном распаде в геологических средах. Реализованная модель позволит использовать код СеИа в задачах оценки безопасности таких сложных объектов, как полигоны глубинной закачки жидких радиоактивных отходов |22|, для которых характерен сильный разогрев пласта в ирискважинной зоне, обусловленный объемным тепловыделением вследствие радиоактивного распада. Также возможным объектом моделирования является проектируемый пункт глубинного захоронения высокоактивных отходов в Нижиекаиском массиве (ПГЗРО) |23|. При исс.ле-

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

СПИСОК ЦИТИРОВАННОЙ ЛИТЕРАТУРЫ

1. Diersh H.-J. G. FEFLOW - Finite element modeling of flow, mass and heat transport in porous and fractured media. Springer, 2014, Berlin Heidelberg, XXXV, 996 pp.

2. Voss C.I. A finite-element simulation model for saturated-unsaturated, fluid-density-dependent ground-water flow with energy transport or chemically-reactive single-species solute transport: U.S. Geological Survey Water-Resources Investigations Report 84-4369. 1984, 409 pp.

3. Voss C.I., Provost A.M. SUTRA, A Model for Saturated-Unsaturated Variable-Density Ground-Water Flow with Solute or Energy Transport: Geological Survey Water-Resources Investigations Report 02-4231. 2002, 250 pp.

4. Weatherhill D,, Simmons C.T., Voss C.E., Robinson N.I. Testing density-dependent ground-water models: twodimensional steady state unstable convection in infinite, finite and inclined porous layers // Advances in Water Resources. 2004. Vol. 27, P. 547-562.

5. Ranganathan V., Hanor J.S. Basin density-driven groundwater flow near salt domes // Chem. Geol. 1988. №74. P. 173-188.

6. Langevin C.D., Thorne D.T., Jr., Dausman A.M., Sukop M.C., Guo, Weixing SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport: U.S. Geological Survey Techniques and Methods Book 6, Chapter A22, 2007, 39 pp.

7. Pruess K,, Oldenburg C,, Moridis G. TOUGH2 User's Guide, ver, 2.0. Lawrence Berkeley National Laboratory 1999.

8. Капырин И. В., Иванов В. А., Копытов Г. В., Уткин С. С. Интегральный код GeRa для обоснования безопасности захоронения РАО // Горный журнал. 2015. №10. С. 44-50.

9. Василевский Ю.В., Конынин И.Н., Копытов Г. В., Терехов К.М. INMOST - программная платформа и графическая среда для разработки параллельных численных моделей на сетках общего вида. М,: Издательство Московского университета, 2013. 144 с.

10. Официальная страница программной платформы INMOST: www.inmost.org.

11. Diersh H.-J.G,, Kolditz О, Coupled groundwater flow and transport: 2, Thermohaline and 3D convection systems // Advances in Water resources, 1998. Vol. 21, P. 401-425.

12. Кпмель Л. P., Машковпч В. П. Защита от ионизирующих излучений. Справочник. Изд. 2. М,: Атомиздат, 1972. 312 с.

13. Nikitin К., Vassilevski Yu, A monotone nonlinear finite volume method for adveetion-diffusion equations on unstructured polyhedral meshes in 3D // Russian J. Numer. Anal. Math. Modelling. 2010. Vol. 25. №4. P.335-358.

14. Haajizadeh M,, Ozguc A. F., Tien C.L. Natural convection in a vertical porous enclosure with internal heat generation // Int. J. Heat Mass Transfer. 1985. Vol. 27, №10. P. 1893-1902.

15. Buretta R. J., Berman A. S. Conveetive heat transfer in a liquid saturated porous layer // ASME J. Appl. Mech. 1976. Vol. 43, P. 249-253.

16. Hardee H. C., Nilson R. H. Natural convection in porous media with heat generation // Nucl. Sci. Engng, 1977. Vol. 63, P. 119-132.

17. Rhee S. J., Dhir V. K,, Catton I. Natural convection heat transfer in beds of inductively heated particles // ASME J. Heat Transfer. 1978. Vol. 100, P. 78-85.

18. Kulacki F.A., Freeman R. G. A note on thermal convection in a saturated, heat generating porous layer // ASME J. Heat Transfer 1979. Vol. 101, P. 169-171.

19. Nield D.A., Bejan A. Convection in Porous Media (5th edition). Springer, 2017.

20. Окуньков Г. А., Рыбапьченко А. И., Куваев А. А. Тепловой режим геологической среды при захоронении жидких радиоактивных отходов // Геоэкология. Инженерная геология. Гидрогеология. Геокриология. 2003. №3. С.237-244.

21. Кирюхин А. В., Куваев А. А. Моделирование тестовой задачи по закачке высокоактивных РАО в водоносный горизонт с учетом радиогенного разогрева пласта // Материалы международной научной конференции «Гидрогеология сегодня и завтра», МГУ. 2013. С.204-209.

22. Рыбапьченко А. И., Пименов М.К., Костин П. П. и др. Глубинное захоронение жидких радиоактивных отходов. М,: ИздАт, 1994. С. 120.

23. Мапьковекий В. И., Пэк А. А., Кочкин Б. Т., Озерский А. Ю. Оценка потенциального загрязнения геологической среды при подземном захоронении радиоактивных отходов на участке «Енисейский» Нижнеканеко-го массива (Красноярский край) // Геоэкология. Инженерная геология. Гидрогеология. Геокриология. 2013. .V"(i. С. 483-490.

24. Левин В. А., Модели и методы. Образование и развитие дефектов. М,: ФИЗМАТЛИТ. 2015. 456 е. (Нелинейная вычислительная механика прочности / Под общ. ред. В. А. Левина: В 5 т. Т. 1).

REFERENCES

1. Diersh, H.-J. G., 2014, FEFLOW - Finite element modeling of flow, mass and heat transport in porous and fractured media., Springer, Berlin Heidelberg, XXXV, 996 pp.

2. Voss, C.I., 1984, A finite-element simulation model for saturated-unsaturated, fluid-density-dependent ground-water flow with energy transport or chemically-reactive single-species solute transport: U.S. Geological Survey Water-Resources Investigations Report 84-4369, 409 pp.

3. Voss, C.I., Provost, A.M. 2002, SUTRA, A Model for Saturated-Unsaturated Variable-Density Ground-Water Flow with Solute or Energy Transport: Geological Survey Water-Resources Investigations Report 02-4231, 250 pp.

4. Weatherhill, D,, Simmons, C.T., Voss, C.E., Robinson N.I. 2004, "Testing density-dependent ground-water models: twodimensional steady state unstable convection in infinite, finite and inclined porous layers", Advances in Water Resources, vol. 27, pp. 547-562.

5. Ranganathan, V., Hanor, J.S. 1988, "Basin density-driven groundwater flow near salt domes" Chem. Geol., No 74, pp. 173-188.

6. Langevin, C.D., Thorne D.T., Jr., Dausman A.M., Sukop M.C., Guo, Weixing, 2007, SEAWAT Version 4: A Computer Program for Simulation of Multi-Species Solute and Heat Transport: U.S. Geological Survey Techniques and Methods Book 6, Chapter A22, 39 pp.

7. Pruess, K,, Oldenburg, C,, Moridis, G,, 1999, TOUGH2 User's Guide, ver. 2.0, Lawrence Berkeley National Laboratory

8. Kapyrin, I. V., Ivanov, V. A., Kopvtov, G. V., Utkin, S. S. 2015, "Integral code GeRa for RAW disposal safety validation" , Gornyi zhurnal, no. 10, pp. 44-50.

9. Vassilevski, Y. V., Konshin I.N., Kopvtov, G.V., Terekhov K.M. 2013, INMOST - programmnaja platforma i graficheskaja sreda dlja razrabotki paralleVnyh chislennyh modelej na setkah obshhego vida [INMOST - a Software Platform and Graphical Environment for Development of Parallel Numerical Models on General Meshes], Izdatelstvo Moskovskogo Universiteta Moscow.

10. INMOST project page, Available at: www.inmost.org.

11. Diersh, H.-J, G,, Kolditz, O, 1998, "Coupled groundwater flow and transport: 2, Thermohaline and 3D convection systems", Advances in Water Resources, vol. 21, pp. 401-425.

12. Kimel, L. R., Mashkovich, V. P. 1972, Z as chit a ot ioniziruyuschih izlucheniy [Protection against ionizing radiation], 2-nd ed,, Atomizdat, Moscow.

13. Nikitin, K,, Vassilevski, Yu. 2010, "A monotone nonlinear finite volume method for adveetion-diffusion equations on unstructured polyhedral meshes in 3D", Russian J. Numer. Anal Math. Modelling, vol. 25, no. 4, pp. 335-358.

14. Haajizadeh, M,, Ozguc, A. F., Tien, C.L. 1985, "Natural convection in a vertical porous enclosure with internal heat generation", Int. J. Heat Mass Transfer, vol. 27, no. 10, pp. 1893-1902.

15. Buretta, R. J., Berman, A. S. 1976, "Conveetive heat transfer in a liquid saturated porous layer" ASME J. Appl. Mech., vol. 43, pp. 249-253.

16. Hardee, H.C., Nilson, R. H. 1977, "Natural convection in porous media with heat generation" Nucl. Sci. Engng., vol. 63, pp. 119-132.

17. Rhee, S. J., Dhir, V. K,, Catton, I. 1978, "Natural convection heat transfer in beds of inductively heated particles" ASME J. Heat Transfer., vol. 100, pp. 78-85.

18. Kulacki, F.A., Freeman, R. G. 1979, "A note on thermal convection in a saturated, heat generating porous layer" ASME J. Heat Transfer, vol. 101, pp. 169-171.

19. Nield, D.A., Bejan, A., Convection in Porous Media (5th edition), Springer, 2017.

20. Okounkov, G. A., Rvbalehenko, A. I., Kuvaev A. A. 2003, "Thermal behaviour of geologic environment used for nuclear waste disposal", Geojekologija. Inzhenernaja geologija. Gidrogeologija. Geokriologija, no. 3, pp. 237-244.

21. Kirjuhin, A. V., Kuvaev, A. A. "Modeling of the test case of high-level radioactive waste injection in aquifer with radiogenic heat", Materialy mezhdunarodnoj nauchnoj konferencii «Gidrogeologija segodnja i zavtra», MGU ["Groundwater geology today and tomorrow" international scientific conference proceedings, MSU], Moscow, 2013, pp. 204-209.

22. Rvbalehenko, A.I., Pimenov, M. K,, Kostin, P.P. et al,, 1994, Glubinnoe zahoronenie zhidkih radioaktivnyh othodov [Deep geological liquid radioactive waste disposal], IzdAt, Moscow.

23. Mal'kovskii, V.I., Рек, A. A., Koehkin, В. Т., Ozerskii, A.Yu. 2013, "Assessment of potential pollution of the geological environment upon underground disposal of radioactive waste at the "Eniseiskii" site in the Nizhnekanskii mountain range (Krasnoyarsk region)", Geojekologija. Inzhenernaja geologija. Gidrogeologija. Geokriologija, no. 6, pp. 483-490.

24. Levin, V. A., 2015, Modeli i metodv. Obrazovanie i razvitie defektov. Nelinejnaja vychisliteVnaja mehanika prochnosti. T.l. [Models and methods. Defect formation and growth. Nonlinear computational mechanics of strength. Vol. 1.], FIZMATLIT, Moscow.

получено 22.05.2017 принято в печать 14.09.2017