Development of heat and mass transfer modulus to Feflow code for calculation of brine loaded to permafrost ground

Texto integral

ВВЕДЕНИЕ

В настоящий момент времени ведется разработка программных продуктов, позволяющих решать задачи геофильтрации в многолетнемерзлых породах (ММП). Как правило, эти коды учитывают процессы теплопереноса без массопереноса. Обзор этих разработок дан в статье [11]. Расчеты фильтрации с учетом тепломассопереноса в ММП опубликованы в работах российских специалистов [7–9]. Эти работы рассматривают одномерные условия фильтрационного потока.

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

Предлагаемый модуль может быть использован при расчете многомерных задач фильтрации и тепломассопереноса. Его особенностью является подключение как внешней библиотеки DLL в среду DHI Feflow [13]. Также дорабатывается возможность подключения расчетов в области ММП к кодам MODFLOW.

ХАРАКТЕРИСТИКА ОБЪЕКТА ИССЛЕДОВАНИЯ

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

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

Для закачки дренажных рассолов в криогенные структуры используются преимущественно участки повышенной трещиноватости. Процесс закачки дренажных вод приводит к высвобождению дополнительного объема пустотности в ММП и образованию техногенного водоносного горизонта. Происходит увеличение емкостных показателей пород и их фильтрационной проницаемости [2, 3]. Дополнительные емкостные параметры ММП на участках закачки формируются главным образом в результате фазового перехода воды из твердого состояния в жидкое при плавлении льда рассолами. Естественная (открытая) пустотность мерзлого массива низкая и не превышает первых процентов от общей потенциальной емкости резервуара.

Как показывает практика и многолетний положительный опыт на разных месторождениях, использование таких участков для удаления дренажных вод – реализованное направление [1–6].

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

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

МЕТОДИКА

Влияние промерзания на проницаемость

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

Сложность задачи заключалась в том, что в программе FEFLOW отсутствовала возможность учета промерзания и изменения теплофизических свойств грунтов в результате этого. Однако такая возможность была реализована в отдельном программном коде на С++ с использованием API функций библиотеки IFM FEFLOW и включена в модель в среде FEFLOW в виде динамически подключаемого модуля FreezeThaw75. Это позволило учесть переход от талого состояния к мерзлому при переменной влагонасыщенности.

Основные уравнения

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

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

$Q_f=\nabla \left(\overrightarrow{q}\right)+n\frac{\partial \left(s_w\right)}{\partial t}+s_w{S}_O\frac{\partial h}{\partial t}.$ (1)

Теплоперенос, распределение температур и теплота плавления льда учитываются в подземной среде и вычисляются на основе (2):

$\nabla \left(\lambda \nabla Т\right)-{С}_w\nabla \overrightarrow{q}T+Q_h=C_a\frac{\partial T}{\partial t}+L_i\frac{\partial {\mathrm{\theta }}_{\mathrm{w}}}{\partial t}.$ (2)

Влагонасыщенность в зависимости от температуры (3) задается с помощью эмпирического уравнения [12]:

$s_w=\left(1-s_{wres}\right)e^{{\left(-T/W\right)}^2}+s_{wres}.$ (3)

Массоперенос оценивается следующим образом (4):

${\theta }_w\nabla \left(D\nabla C\right)-\nabla \overrightarrow{q}C+Q_t=\frac{\partial \left({\theta }_{w}C\right)}{\partial t}.$ (4)

Вариант вычисления относительной проницаемости в зависимости от влагонасыщенности (5):

$K_{rw}={10}^{-\Omega \left(1-s_w\right)n}.$ (5)

Коэффициенты теплопроводности и теплоемкости для твердой фазы, учитывающие льдистость, вычисляются по (6):

$\lambda =n\cdot s_w{\lambda }_w+n\left(1-s_w\right){\lambda }_i+\left(1-n\right){\lambda }_s,$

$C_a=n\cdot s_w{\rho }_w{C}_w+n\left(1-s_w\right){\rho }_i{C}_i+\left(1-n\right){\rho }_s{C}_s.$(6)

В уравнениях (1)–(6):

$q=-K_{rw}K\nabla \left(\frac{{\rho }_w}{\mu }\left(\overrightarrow{h}+\frac{{\rho }_w-{\rho }_0}{{\rho }_w}\right)\right)g$ −

скорость Дарси, [м/сут]; n − пористость, [д.е.]; h − напор, [м]; T − температура, [K]; ρ_o,ρ_w − плотности для пресной воды и воды в зависимости от концентрации солей, [кг/м³]; ρ_i,ρ_s − плотности для льда и породы, [кг/м³]; g − ускорение силы тяжести, [м/с²]; S_{w =} θ_w/n − влагонасыщенность, [д.е.]; θ_w – объемная влажность, [д.е.]; s_wres − остаточная влагонасыщенность, [д.е.]; S_o − удельная упругоемкость, [1/м]; µ − динамическая вязкость жидкой фазы, [Па с]; L_i − удельная теплота плавления льда, [Дж/кг]; C − концентрация солей в жидкой фазе, [мг/л]; D − тензор гидродинамической дисперсии, [м²/сут]; Ω − показатель сопротивления, д.е.]; W − эмпирический параметр, задающий интервал фазового перехода между твердым и талым состояниями воды, [К];  K_rw, K− относительная проницаемость и собственная проницаемость породы, соответственно [д.е.], [м²]; λ_w, λ_i, λ_s − коэффициент теплопроводности воды, льда и породы, [Вт/м/К]; C_w, C_a, C_i, C_s − теплоемкость воды, твердой фазы, льда, породы, [Дж/м/К]; Для проверки достоверности результатов, полученных с помощью разработанного программного средства, были выполнены тестовые расчеты, результаты которых приведены ниже.

РЕЗУЛЬТАТЫ ИССЛЕДОВАНИЙ И ИХ ОБСУЖДЕНИЕ

Промерзание при отсутствии фильтрационного потока

Для постановки данных условий учитываются фазовые переходы в насыщенной пористой среде при отсутствии гидродинамического потока. Для таких условий имеется аналитическое решение Лунардини [11].

На рис. 1 приведена постановка этой задачи. Изначально в горизонтальной модельной области задана температура +4°C. На левой границе задана постоянная температура –6°C, приводящая к постепенному промерзанию. В аналитической модели оцениваются термодинамическое влияние левой границы и распространение промерзания во времени. На рис. 2 приведено сопоставление расчетов на основе модуля FreezeThaw75 с аналитическим решением на момент времени 1 сут.

Рис. 1. Постановка задачи Лунардини.

Рис. 2. Сопоставление аналитического расчета и расчета с помощью модуля FreezeThaw75.

Промерзание при наличии фильтрационного потока

Такая модельная задача была предложена при разработке программы SUTRA [10, 11]. Двумерная напорная фильтрация решается в прямоугольной области размером 10 м по вертикали и 50 м по горизонтали при наличии слабопроницаемой стены в грунте с постоянной температурой ниже точки замерзания воды и являющейся преградой фильтрационному потоку. Течение жидкости вызвано перепадом давления между левой и правой границей, на которых заданы ГУ1 с постоянным напором и перепадом 0.05 м между ними. Постоянная температура на левой границе принимается равной +5°C. Стена в грунте внутри модельной области поддерживается при температуре –5°C. Начальная температура всей модельной области +5°C. Дискретизация – равномерная прямоугольная с шагом 0.5 м. Продолжительность моделирования устанавливается равной 800 сут, в течение которых достигается равновесное состояние. Результаты сопоставления расчетов температур и насыщенностей, выполненных с помощью разработанного модуля и программой SUTRA, показаны на рис. 3. Расхождения можно объяснить особенностями сеточной аппроксимации каждой из программ.

Рис. 3. Сопоставление расчетов на модуле FreezeThaw75 (слева) и SUTRA (справа). Верхний ряд – сопоставление по температуре, нижний – по влагонасыщенности.

Тепломассоперенос, фазовые переходы и фильтрационный поток

Для верификации работы модуля на задачах, где учитывается тепломассоперенос и фазовые переходы, было произведено сопоставление с расчетами, проведенными авторами статьи [12]. В статье [12] рассматривалась профильная модель, разработанная в программе FlexPDE. В ней учтена зона ММП небольшой мощности и характерные для нее фильтрационные и температурные условия. Задана толща сезонно-талого слоя и талики. Один из таликов образовался за счeт утечек муниципальных сточных вод (содержащих в своем составе хлориды), которые разгружаются в другой талик, образованный урезом реки. За счет минерализации снижается температура замерзания, которая изначально равна 0°C для пресной воды. Повышенные концентрации растворенных веществ снижают температуру замерзания поровой воды. Процесс оттаивания, как под действием тепла, так и из-за термоосмотического сдвига за счeт растворенных в воде веществ уменьшает содержание льда в ММП и увеличивает их проницаемость.

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

Рис. 4. Начальные и граничные условия для задачи тепломассопереноса. Сине-голубой заливкой показана начальная льдистость, изолиниями – начальные напоры, красной – начальная концентрация промстоков 100 мг/л.

На рис. 5 приведено сопоставление расчетов тепломассопереноса с помощью модуля FreezeThaw75 и результатов расчета из статьи [12].

Рис. 5. Сопоставление расчетов тепломассопереноса в FreezeThaw75 (слева) и расчeтов из [12] (справа). Сверху льдистость, внизу концентрация хлоридов.

Расчеты выполнены при одинаковых условиях по фильтрации, тепло- и массопереносу. В качестве мигрирующего компонента рассматривается сценарий миграции хлоридов с максимальной концентрацией 100 мг/л. Если вода имеет повышенную минерализацию, то происходит депрессия температуры замерзания, которая может быть оценена следующим образом. Зависимость температуры замерзания T_f (°C) от минерализации C (мг/л) предлагается рассмотреть по следующей зависимости [12]:

$T_f=\left\{\begin{array}{c}C\ge 50\text{\hspace{0.33em}мг}/\text{}\mathrm{л},-0.08C+4\\ С<50\mathrm{мг}\text{}/\text{}\mathrm{л},0\end{array}\right..$ (7)

Различие в расчетах объясняется тем, что в статье [12] используется динамическое адаптивное измельчение сетки. Размеры ячеек и дискретность по вертикали во время моделирования могут варьироваться на несколько порядков.

Результаты моделирования режима закачки крепких рассолов на опытном участке

После верификационных расчeтов модуль FreezeThaw75 опробован на одном из участков отведения рассолов в многолетнемерзлые породы в Далдыно-Алакитском районе Якутии.

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

Методы построения фильтрационных моделей

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

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

Опыт построения фильтрационных моделей одного и того же участка на основе различного программного обеспечения показал, что при едином подходе получены расхождения в оценке конечных объемов закачки. Для расчетов использовали ПО Modflow/GMS и FEFLOW. Закачка рассолов осуществлялась в ММП. Калибровка выполнялась в нестационарных условиях.

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

Модель тепломассопереноса для участка закачки

Второй подход при моделировании применялся авторами и заключался в применении расчетов с учетом тепломассопереноса (ТМП) и фазовых переходов лед–вода и вода–лед. При расчете использовался модуль FreezeThaw75.

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

Моделирование закачки рассолов

В соответствии с зональностью фильтрационных свойств ММП в разрезе и необходимостью учета напоров верхнекембрийского водоносного комплекса (ВВК) в подмерзлотной толще пород, модель разделили на две части: верхняя – ММП, нижняя – ВВК (рис. 6, 7). Начальный уровень свободной поверхности +165 м. По границам заданы непроницаемые граничные условия. Влияние ручьев и инфильтрации на образующийся техногенный водоносный горизонт не учитываются, так как на участке закачки отсутствуют сквозные таликовые зоны.

Рис. 6. Схема расположения объекта моделирования с системой разломов (красные линии). Скважины: закачки − красные кружки, наблюдательные – зеленые; цифры − номера скважин. Линия разреза показана черным цветом.

Рис. 7. Схема разреза. Слева – южная часть модели, справа – восточная. Черные линии − скважины, белая линия – кривая репрессии, заливка – коэффициент фильтрации через 3 мес. после начала закачки.

Объем закачанных рассолов за весь период составил почти 22 млн м³. Для закачки были пробурены 18 скважин: 12 – западного ряда, и 6 – восточного. Заданы расходы в соответствии с режимом закачки на 10-летний период (см. рис. 6). Температура рассолов задавалась постоянной –5°C. Концентрация закачиваемого рассола до 350 г/л. Фоновая концентрация рассолов ВВК рассолов 100 г/л. Фоновая температура массива ММП –2°C. Для минимизации численной дисперсии модельные слои были раздроблены по вертикали.

В процессе моделирования выполнялась калибровка параметров. Проводилось сопоставление фактических и расчетных уровней. Сопоставление фактических и расчетных уровней по моделям фильтрации и ТМП показано на рис. 8. Отмечается более качественное совпадение фактических и расчетных уровней для модели с применением ТМП. Далее этот подход применялся для оценки оставшейся емкости после закачки, и был проведен расчет эксплуатации участка закачки. Результаты расчета показаны на рис. 9.

Рис. 8. Сопоставление фактических и расчетных уровней по ближайшим к закачным скважинам пьезометрам №№ 1–4.

Рис. 9. Результаты расчета закачки в течение 2013–2022 гг. Желтым показаны разломы, зеленым − ореол минерализации 200 г/л, белым − не размороженная часть.

ЗАКЛЮЧЕНИЕ

Разработанный модуль FreezeThaw75 позволяет решать задачи тепломассопереноса и фазовых переходов (лед–вода и вода–лед).

Для проверки достоверности результатов было выполнено сопоставление с:

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

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

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

Sobre autores

I. Rastorguev

Autor responsável pela correspondência
Email: stiar@mail.ru
Rússia, pl. Kurchatova, 1, Moscow, 123098; ul. Sovetskaya, 5, Novosibirsk, 630007

I. Litvinova

Email: LitvinovaIV@alrosa.ru
Rússia, ul. Sovetskaya, 5, Novosibirsk, 630007

N. Iost

Email: LitvinovaIV@alrosa.ru
Rússia, ul. Sovetskaya, 5, Novosibirsk, 630007

A. Ilin

Email: LitvinovaIV@alrosa.ru
Rússia, ul. Sovetskaya, 5, Novosibirsk, 630007

Bibliografia

Arquivos suplementares

Arquivos suplementares

Ação

1. JATS XML

Baixar

2. Fig. 1. Statement of the Lunardini problem.

Baixar (128KB)

Metadados

3. Fig. 2. Comparison of analytical calculation and calculation using the FreezeThaw75 module.

Baixar (176KB)

Metadados

4. Fig. 3. Comparison of calculations on the FreezeThaw75 module (left) and SUTRA (right). The top row is a comparison by temperature, the bottom row is a comparison by moisture saturation.

Baixar (374KB)

Metadados

5. Fig. 4. Initial and boundary conditions for the heat and mass transfer problem. The blue fill shows the initial ice content, the isolines show the initial pressures, and the red fill shows the initial concentration of industrial wastewater of 100 mg/l.

Baixar (209KB)

Metadados

6. Fig. 5. Comparison of heat and mass transfer calculations in FreezeThaw75 (left) and calculations from [12] (right). Ice content at the top, chloride concentration at the bottom.

Baixar (192KB)

Metadados

7. Fig. 6. Layout of the modeling object with a fault system (red lines). Wells: injections - red circles, observation - green; numbers - well numbers. The section line is shown in black.

Baixar (256KB)

Metadados

8. Fig. 7. Section diagram. On the left is the southern part of the model, on the right is the eastern part. Black lines are wells, white line is the repression curve, filling is the filtration coefficient 3 months after the start of injection.

Baixar (302KB)

Metadados

9. Fig. 8. Comparison of actual and calculated levels for piezometers No. 1–4 closest to injection wells.

Baixar (445KB)

Metadados

10. Fig. 9. Results of injection calculations for 2013–2022. Faults are shown in yellow, the 200 g/l mineralization halo is shown in green, and the unfrozen part is shown in white.

Baixar (135KB)

Metadados