Теплофизика высоких температур, 2022, T. 60, № 3, стр. 434-442

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

М. А. Семин^1, *, Л. Ю. Левин¹, М. С. Желнин², О. А. Плехов²

¹ Горный институт УрО РАН
Пермь, Россия

² Институт механики сплошных сред УрО РАН
Пермь, Россия

^* E-mail: seminma@inbox.ru

Поступила в редакцию 14.12.2020
После доработки 13.04.2021
Принята к публикации 19.05.2021

EDN: BVVLYH
DOI: 10.31857/S0040364422020296

Полный текст (PDF)

Аннотация

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

ВВЕДЕНИЕ

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

Экспериментальные исследования минерализации подземных вод должны сопровождаться теоретическими и лабораторными исследованиями закономерностей замораживания пород с содержанием минерализованных вод. Известно, что c повышением минерализации подземных вод температура кристаллизации воды (или в данном случае рассола) уменьшается [2–4], что может негативно сказаться на несущей способности ледопородных ограждений строящихся шахтных стволов [5, 6]. Тем не менее в настоящее время в литературе представлено небольшое количество исследований влияния минерализации подземных вод на процесс искусственного замораживания пород [7–10]. Сравнительно большее количество работ посвящено исследованию естественных процессов замораживания и оттаивания пород в условиях криолитозоны [11–17], однако данные исследования не могут быть в полной мере распространены на случай искусственного замораживания пород. Это связано с отличительными особенностями процесса искусственного замораживания пород: более высокие температурные перепады в породах, более быстрое изменение температуры и льдистости пород с течением времени, более высокие горное и гидростатическое давления вследствие большей глубины.

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

АНАЛИЗ ЛИТЕРАТУРЫ

При замораживании породного массива, содержащего минерализованную воду (рассол), происходит ряд взаимосвязанных тепломассообменных процессов, впервые классифицированных в работе [7]:

1) рост концентрации рассола вблизи фронта фазового перехода и связанное с этим снижение температуры кристаллизации содержащейся в рассоле воды;

2) диффузионный массоперенос соли, вызванный градиентом концентрации соли в рассоле;

3) конвективный перенос (или выпирание) незамерзшего рассола от границы фазового перехода вследствие объемного расширения воды при замерзании (~9%);

4) естественная конвекция, обусловленная разницей плотностей рассола в зоне охлаждения вблизи фронта фазового перехода и на удалении от него;

5) термодиффузия рассола (эффекты Людвига–Соре и Дюфура [18, 19]);

6) изменение теплофизических свойств породного массива при изменении фазового и компонентного составов рассола в массиве.

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

Несмотря на большое количество идущих в замораживаемом породном массиве тепломассообменных процессов, в настоящее время не представлено исследований, в которых бы производился их комплексный учет. Чаще всего в литературе описываются исследования на математических моделях с учетом только процессов 1, 2 [7, 8, 15, 20]. В работах [9, 10] к рассмотрению также добавляется процесс 3. Процессы 4, 5 в существующих работах рассматриваются применительно к не минерализованной поровой воде переменной плотности из-за разности температур [21–25]. Процесс 6 на сегодняшний день изучен в наименьшей степени, а существующие исследования содержат преимущественно качественный анализ влияния минерализации воды на теплофизические свойства мерзлых пород в естественном залегании [26]. Следует отметить, что влияние температуры пород на их теплофизические свойства в литературе изучены достаточно подробно [6].

Представляется, что при анализе тепломассопереноса в искусственно замораживаемых породах с содержанием рассолов важно учитывать и такой физический процесс, как выпадение соли в твердый нерастворимый осадок при достижении заданной отрицательной температуры. Это связано с тем, что при строительстве шахтных стволов породный массив зачастую замораживается до достаточно низких температур (менее –25°С) [1, 27]. По данным [1], соль может присутствовать в породах как в растворенной форме в воде, так и в форме твердого осадка. Для случая насыщенного раствора NaCl (при минерализации 290 г/л) полное замерзание с выпадением всей соли в твердый осадок происходит при –21.2°C, а для насыщенного раствора MgCl₂ эта температура достигает –33.6°C. Известно, что в процессе замораживания рассола содержащаяся в нем соль выпадает в осадок постепенно в таком количестве, чтобы при охлаждении насыщенного раствора не происходило его перенасыщения.

В литературе описан ряд упрощенных подходов к учету данного процесса. В работах [9–11] принято, что вся соль растворяется и остается в растворенной виде в рассоле, если концентрация рассола равна или меньше концентрации насыщенного рассола; в противном случае в осадок выпадет избыток соли в таком количестве, чтобы концентрация рассола не смогла превысить своего предельного значения, соответствующего насыщенному рассолу. В [7, 15] предполагается, что концентрация соли всегда ниже концентрации насыщенного раствора, при этом вся растворенная соль при образовании льда отторгается в раствор и в осадок не выпадает – это допущение уместно при рассмотрении малых вариаций температур вблизи нуля.

Для количественного описания процесса 1 из представленной выше классификации в литературе представлено несколько различных методов и подходов. Глобально можно выделить два подхода. Первый подход не учитывает переноса соли в пористом породном массиве и основан на использовании эмпирической кривой замораживания – в данном случае концентрация незамерзшего рассола зависит напрямую от температуры [28]. Второй подход связан с явным рассмотрением переноса соли в замораживаемом породном массиве, а температура кристаллизации рассола является функцией концентрации растворенной в нем соли – линейной [15] или квадратичной [7]. Еще более сложная полиномиальная зависимость от двух переменных (концентрации соли и температуры пород) введена в работе [9] на основании данных лабораторного эксперимента. В рамках второго подхода требуется явное рассмотрение уравнений переноса растворенной соли.

При моделировании процесса 2 в литературе обычно используется закон Фика, выражающий линейную зависимость массового потока растворенной в рассоле соли от градиента концентрации соли [7, 15]. Коэффициент диффузии в законе Фика зависит от многих факторов (пористость, проницаемость, вид растворенного вещества и растворителя, температура и пр.), однако чаще всего он предполагается постоянным в диапазоне 10^–14–10^–9 м²/с [9, 14, 15].

Еще одним интересным физическим процессом, особенно важным при рассмотрении естественного замораживания склонных к пучению засоленных грунтов, является наличие капиллярных сил [9, 13, 29]. Для задания зависимости насыщенности от капиллярного давления в [9] использована модель Ван Генухтена [30], в то время как в [29] использовалась модель Леверетта [31]. В [13] использовался несколько иной подход, в рамках которого температура замерзания поровой воды связывалась с такими параметрами, как кривизна и поверхностное натяжение границы раздела твердое тело–жидкость. В условиях замораживания частично водонасыщенных породных массивов и грунтов капиллярное давление приводит к дополнительной миграции влаги к фронту фазового перехода (криогенные течения) [29].

На сегодня существует два подхода к моделированию тепломассопереноса в грунтах и породных массивах, насыщенных минерализованной водой – первый подход заключается в явном выделении фронта фазового перехода [15], а второй (сквозного счета) в неявном рассмотрении фронта с использованием дополнительных функций температуры, в которые “зашивается” информация о фазовом переходе воды: энтальпии [7], кажущейся плотности [9], объемной доли незамерзшей воды или льда [8] и пр. Наиболее общий подход к построению уравнения состояния для описания фазового перехода поровой воды во взаимосвязи с ее минерализацией и гидростатическим давлением на основе свободной энергии Гиббса описан в [9, 10].

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

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

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

− свойства сухого скелета породного массива считаются однородными и изотропными;

− в начальный момент времени распределение всех параметров среды в объеме породного массива также считается однородным;

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

− фазовый переход испытывает только вода, в то время как соль, содержащаяся в рассоле, остается в нем;

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

− фазовый переход воды происходит полностью в заданном конечном интервале температур;

− температуры начала замерзания и полного замерзания воды в рассоле (или температуры ликвидуса и солидуса) линейно зависят от концентрации соли;

− конвективный перенос рассола не рассматривается;

− диффузионный перенос растворенной соли в рассоле происходит в соответствии с законом Фика;

− предполагается локальное тепловое равновесие фаз и компонентов в каждой точке рассматриваемой среды.

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

(1)

$\begin{gathered} n{{S}_{b}}{{\rho }_{b}}\frac{{\partial c}}{{\partial t}} + nc{{\rho }_{b}}\frac{{\partial {{S}_{b}}}}{{\partial t}} = \\ = \,\,\frac{1}{r}\frac{\partial }{{\partial r}}\left( {r{{\rho }_{b}}D\frac{{\partial c}}{{\partial r}}} \right) - {{\rho }_{b}}k\left( {{{c}_{{{\text{sat}}}}} - c} \right), \\ \end{gathered} $

(2)

$n{{\rho }_{s}}\frac{{\partial {{S}_{s}}}}{{\partial t}} = {{\rho }_{b}}k\left( {{{c}_{{{\text{sat}}}}} - c} \right),$

(3)

${{S}_{i}} + {{S}_{b}} + {{S}_{s}} = 1,$

(4)

$\rho {{C}_{p}}\frac{{\partial T}}{{\partial t}} + {{\rho }_{w}}nL\frac{{\partial {{S}_{b}}}}{{\partial t}} = \frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\lambda \frac{{\partial T}}{{\partial r}}} \right).$

Здесь $n$ – пористость массива; ${{\rho }_{b}}$ – истинная плотность рассола, кг/м³; ${{S}_{b}}$ – объемная концентрация рассола в поровом пространстве породного массива, м³/м³; $t$ – время, с; $c$ – объемная концентрация растворенной соли, м³/м³; $~{{\rho }_{i}}$ – истинная плотность льда, кг/м³; ${{S}_{i}}$ – объемная концентрация льда в поровом пространстве породного массива, м³/м³; $r$ – радиальная координата, м; $D$ – коэффициент диффузии растворенной соли, м²/с; $k$ – параметр, характеризующий инерционность процесса выпадения соли в твердый осадок; ${{c}_{{{\text{sat}}}}}$ – объемная концентрация насыщенного раствора соли, м³/м³; ${{\rho }_{s}}$ – истинная плотность соли, выпавшей в твердый осадок, кг/м³; ${{S}_{s}}$ – объемная концентрация выпавшей в твердый осадок соли в поровом пространстве породного массива, м³/м³; $\rho $ – плотность породного массива, содержащего воду, лед и соль, кг/м³; ${{C}_{p}}$ – удельная теплоемкость породного массива, Дж/(кг °С); $T$ – температура, °С; ${{\rho }_{w}}$ – истинная плотность воды, кг/м³; $L$ – удельная теплота кристаллизации воды, Дж/кг; $\lambda $ – теплопроводность породного массива, Вт/(м °С).

Уравнение (1) – диффузионное уравнение для объемной концентрации соли $c$, растворенной в рассоле; (2) – уравнение прироста объемной концентрации выпавшей в осадок соли в поровом пространстве массива; (3) – уравнение баланса массы в поровом пространстве массива; (4) – уравнение баланса теплоты в системе в целом. Переменные, отвечающие за компонентный и фазовый состав в элементарном объеме породного массива, определяются для некоторого произвольного малого объема замораживаемого влажного массива (рис. 1) согласно следующим формулам:

(5)

${{S}_{b}} = \frac{{{{V}_{b}}}}{{{{V}_{p}}}} = \frac{{{{V}_{w}} + {{V}_{{ds}}}}}{{{{V}_{p}}}},$

(6)

${{S}_{i}} = \frac{{{{V}_{i}}}}{{{{V}_{p}}}},$

(7)

${{S}_{s}} = \frac{{{{V}_{{ps}}}}}{{{{V}_{p}}}},$

(8)

$c = \frac{{{{V}_{{ds}}}}}{{{{V}_{b}}}}.$

Здесь ${{V}_{p}}$ – объем порового пространства, м³; ${{V}_{b}}$ – объем жидкого рассола (вода + растворенная соль), м³; ${{V}_{w}}$ – объем чистой воды, м³; ${{V}_{{ds}}}$ – объем растворенной соли, м³; ${{V}_{i}}$ – объем льда, м³; ${{V}_{{ps}}}$ – объем выпавшей в осадок соли, м³.

Рис. 1.

Компонентный и фазовый состав элементарного объема замораживаемого породного массива.

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

(9)

$\begin{gathered} S_{b}^{*} = \frac{{{{S}_{b}}}}{{{{S}_{i}} + {{S}_{b}}}} = \\ = \,\,\left\{ {\begin{array}{*{20}{c}} {1,} \\ {\left( {T - {{T}_{{{\text{sol}}}}}} \right){\text{/}}\left( {{{T}_{{{\text{liq}}}}} - {{T}_{{{\text{sol}}}}}} \right),} \\ {0,} \end{array}} \right.\,\,\,\,\begin{array}{*{20}{c}} {T > {{T}_{{{\text{liq}}}}};} \\ {{{T}_{{{\text{liq}}}}} \geqslant T > {{T}_{{{\text{sol}}}}};} \\ {{{T}_{{{\text{sol}}}}} \geqslant T.} \end{array} \\ \end{gathered} $

Здесь ${{T}_{{{\text{liq}}}}}$ – температура ликвидуса, °С; ${{T}_{{{\text{sol}}}}}$ – температура солидуса, °С; $S_{i}^{*}$ – относительная объемная концентрация льда, рассчитанная для доли порового пространства, занятого только льдом и водой без примеси соли. Решение и последующий анализ системы уравнений проводятся в терминах переменных ${{S}_{i}}$ и ${{S}_{b}}$, а не $S_{b}^{*}$.

Температуры ликвидуса и солидуса зависят от объемной концентрации растворенной соли согласно следующим линейным функциям:

(10)

${{T}_{{{\text{liq}}}}} = T_{{{\text{liq}}}}^{{\left( 0 \right)}} - \beta c,$

(11)

${{T}_{{{\text{sol}}}}} = T_{{{\text{sol}}}}^{{\left( 0 \right)}} - \beta c,$

где $\beta $ – коэффициент понижения температуры фазового перехода вода–лед, °С. В общем случае рассолов, включающих в себя различные виды соли, коэффициенты понижения температуры фазового перехода для температур солидуса и ликвидуса будут, скорее всего, не равны. Однако здесь рассматривается соль однородного состава.

Теплофизические свойства породного массива, присутствующие в системе балансовых уравнений (1)–(3), задаются как функции пористости, объемных концентраций рассола, льда, растворенной и выпавшей в осадок соли следующими формулами:

(12)

$\begin{gathered} \rho {{C}_{p}} = \left( {1 - n} \right){{\rho }_{r}}{{C}_{r}} + \\ + \,\,n\left[ {{{S}_{i}}{{\rho }_{i}}{{C}_{i}} + {{S}_{b}}{{\rho }_{w}}{{C}_{w}}\left( {1 - c} \right) + {{\rho }_{s}}{{C}_{s}}\left( {{{S}_{s}} + c{{S}_{b}}} \right)} \right], \\ \end{gathered} $

(13)

$\lambda = \lambda _{r}^{{1 - n}}\lambda _{i}^{{n{{S}_{i}}}}\lambda _{w}^{{n{{S}_{b}}\left( {1 - c} \right)}}\lambda _{s}^{{n\left( {{{S}_{s}} + c{{S}_{b}}} \right)}},$

(14)

${{\rho }_{b}} = {{\rho }_{w}}\left( {1 - c} \right) + {{\rho }_{s}}c.$

Здесь $C$ – удельная теплоемкость рассматриваемой фазы, Дж/(кг °С); индексы r, i, b, w и s соответствуют сухому скелету породного массива, льду, незамерзшему рассолу, чистой воде и выпавшей в осадок соли соответственно. Данные формулы являются обобщением известных формул [32, 33] на случай наличия дополнительной фазы – соли, выпавшей в твердый осадок.

Коэффициент диффузии растворенной соли рассчитывается как линейная функция [15, 34] эффективной пористости породного массива $n{\kern 1pt} *$, представляющая собой часть объема, занятую незамерзшим рассолом:

(15)

$D = {{D}_{0}}n{\kern 1pt} * = {{D}_{0}}n{{S}_{b}}.$

Здесь ${{D}_{0}}$ – коэффициент диффузии растворенной соли в рассоле в ситуации, когда пористая среда отсутствует и рассол занимает весь элементарный объем, м²/с.

Вся информация о фазовом переходе воды в лед “зашита” во втором слагаемом слева в уравнении (4). В данном случае отсутствует явное выделение фронта фазового перехода, а сам подход к решению задачи теплопереноса аналогичен подходу эффективной теплоемкости [35] и энтальпийному подходу [27, 28].

На границе с замораживающей колонкой $r = {{r}_{c}}$ задается граничное условие III рода теплообмена между движущимся хладоносителем и примыкающим к колонке породным массивом:

(16)

$\lambda \frac{{\partial T}}{{\partial r}}\left( {t,{{r}_{c}}} \right) = \alpha \left( {{{T}_{c}} - T\left( {t,{{r}_{c}}} \right)} \right),$

а также условие нулевого градиента концентраций растворенной соли

(17)

$\frac{{\partial c}}{{\partial r}}\left( {t,{{r}_{c}}} \right) = 0.$

Здесь $\alpha $ – коэффициент теплоотдачи от породного массива к хладоносителю в колонке, зависящий от геометрических параметров колонки, теплофизических и гидравлических свойств рассола [27, 36], Вт/(м² °С); ${{T}_{c}}$ – температура хладоносителя в колонке, °С.

На внешней границе расчетной области $r = {{r}_{{{\text{out}}}}}$ задается условие, соответствующее непотревоженному состоянию породного массива:

(18)

$T\left( {t,{{r}_{{{\text{out}}}}}} \right) = {{T}_{0}} > {{T}_{{{\text{liq}}}}},$

(19)

$c\left( {t,{{r}_{{{\text{out}}}}}} \right) = {{c}_{0}}.$

Здесь ${{T}_{0}}$ – начальная температура породного массива, °С; ${{c}_{0}}$ – начальная объемная концентрация растворенной соли в поровом пространстве породного массива, м³/м³.

Задача дополняется соответствующими начальными условиями:

(20)

$T\left( {0,r} \right) = {{T}_{0}},$

(21)

$c\left( {0,t} \right) = {{c}_{0}}.$

Получена замкнутая система уравнений (1)–(21) со следующими неизвестными функциями: $c,~{{S}_{b}},~{{S}_{i}},{{S}_{s}},~T$. Точное решение данной задачи может быть получено только численно.

ЧИСЛЕННЫЙ АЛГОРИТМ

Для численной реализации задачи (1)–(21) использовался метод конечных разностей. Расчетная область разбивалась на ячейки одинакового размера ${{\Delta }}r$. Шаг по времени рассчитывался на основании известного условия Куранта–Фридрихса–Леви для уравнения диффузии [35]:

(22)

$\Delta t = {\text{CFL}}\frac{{\Delta {{r}^{2}}}}{{2a}},$

где а – минимальное значение из коэффициентов тепловой диффузии в среде и молекулярной диффузии соли, растворенной в рассоле; ${\text{CFL}}$ – число Куранта–Фридрихса–Леви [37].

На каждой временной итерации k новых значений неизвестных функций $c,~{{S}_{b}},~{{S}_{i}},{{S}_{s}},~T$ определялись с использованием следующего алгоритма, включающего в себя четыре последовательных этапа.

Этап 1. Расчет прироста удельной энтальпии $\Delta {{h}_{j}}$ в результате действия только кондуктивного теплопереноса:

(23)

$\begin{gathered} \Delta {{h}_{j}} = \frac{{\Delta t}}{{\Delta r}}\left( {{{\lambda }_{{j + 1/2}}}\frac{{T_{{j + 1}}^{{\left( k \right)}} - T_{j}^{{\left( k \right)}}}}{{\Delta r}} - {{\lambda }_{{j - 1/2}}}\frac{{T_{j}^{{\left( k \right)}} - T_{{j - 1}}^{{\left( k \right)}}}}{{\Delta r}}} \right) + \\ + \,\,\Delta t\frac{{{{\lambda }_{j}}}}{r}\frac{{T_{{j + 1}}^{{\left( k \right)}} - T_{{j - 1}}^{{\left( k \right)}}}}{{2\Delta r}}~, \\ \end{gathered} $

(24)

${{\lambda }_{{j \pm 1/2}}} = \frac{{{{\lambda }_{{j \pm 1}}} + {{\lambda }_{j}}}}{2}.$

Этап 2. Расчет изменений температуры и объемной концентрации рассола, вызванных изменением удельной энтальпии $\Delta {{h}_{j}}$:

(25)

$\begin{gathered} \left( {\rho {{C}_{p}}} \right)_{j}^{{\left( k \right)}}\left( {T_{j}^{{\left( {k + 1} \right)}} - T_{j}^{{\left( k \right)}}} \right) + \\ + \,\,{{\rho }_{w}}nL\left( {S_{{bj}}^{{\left( {k + 1} \right)}} - S_{{bj}}^{{\left( k \right)}}} \right) = \Delta {{h}_{j}}, \\ \end{gathered} $

(26)

$S_{{bj}}^{{\left( {k + 1} \right)}} = S_{b}^{*}\left( {T_{j}^{{\left( {k + 1} \right)}},c_{j}^{{\left( k \right)}}} \right)\left( {1 - S_{{sj}}^{{\left( k \right)}}} \right).$

После решения данной нелинейной системы уравнений относительно $T_{j}^{{\left( {k + 1} \right)}}$ и ${{S}_{b}}_{j}^{{\left( {k + 1} \right)}}$ определяются узловые значения концентрации растворенной соли на следующем временном шаге:

(27)

$n{{S}_{b}}_{j}^{{\left( k \right)}}{{\rho }_{b}}\left( {c_{j}^{*} - c_{j}^{{\left( k \right)}}} \right) + nc_{j}^{{\left( k \right)}}{{\rho }_{b}}\left( {{{S}_{b}}_{j}^{{\left( {k + 1} \right)}} - {{S}_{b}}_{j}^{{\left( k \right)}}} \right) = 0.$

Если новое значение концентрации $c{\kern 1pt} *_{j}^{{\left( {k + 1} \right)}}$ превышает концентрацию насыщенного раствора соли, то часть соли выпадает в твердый осадок таким образом, чтобы выполнилось равенство

(28)

$c_{j}^{*} = {{c}_{{{\text{sat}}}}}.$

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

(29)

${{S}_{s}}_{j}^{{\left( {k + 1} \right)}} = {{S}_{s}}_{j}^{{\left( k \right)}} + {{S}_{b}}_{j}^{{\left( {k + 1} \right)}}\left( {c_{j}^{*} - {{c}_{{{\text{sat}}}}}} \right),$

(30)

${{S}_{i}}_{j}^{{\left( {k + 1} \right)}} = 1 - {{S}_{b}}_{j}^{{\left( {k + 1} \right)}} - {{S}_{s}}_{j}^{{\left( {k + 1} \right)}}.$

Этап 3. Расчет диффузионного переноса растворенной соли:

(31)

$\begin{gathered} \Delta c_{j}^{{\left( k \right)}} = \,\,~\frac{{\Delta t}}{{\Delta r}}\left( {{{D}_{{j + 1/2}}}\frac{{c_{{j + 1}}^{*} - c_{j}^{*}}}{{\Delta r}} - {{D}_{{j - 1/2}}}\frac{{c_{j}^{*} - c_{{j - 1}}^{*}}}{{\Delta r}}} \right) + \\ + \,\,\Delta t\frac{{{{D}_{j}}}}{r}\frac{{c_{{j + 1}}^{*} - c_{{j - 1}}^{*}}}{{2\Delta r}}, \\ \end{gathered} $

(32)

${{D}_{{j \pm 1/2}}} = {{D}_{0}}n\frac{{{{S}_{b}}_{{j \pm 1}}^{{\left( k \right)}} + {{S}_{b}}_{j}^{{\left( k \right)}}}}{2} + {{D}_{{{\text{art}}}}}n.$

Здесь ${{D}_{{{\text{art}}}}}$ – коэффициент искусственной диффузии, м²/с.

Этап 4. Расчет граничных условий. Расчет температуры в граничном узле 0 вблизи замораживающей колонки осуществляется по формуле

(33)

$\lambda \frac{{T_{0}^{{\left( {k + 1} \right)}} - T_{1}^{{\left( {k + 1} \right)}}}}{{\Delta t}} = \alpha \left( {{{T}_{c}} - T_{0}^{{\left( {k + 1} \right)}}} \right).$

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

(34)

$\begin{gathered} T_{N}^{{\left( {k + 1} \right)}} = {{T}_{0}},\,\,\,\,c_{N}^{{\left( {k + 1} \right)}} = {{c}_{0}}, \\ {{S}_{i}}_{N}^{{\left( {k + 1} \right)}} = {{S}_{s}}_{N}^{{\left( {k + 1} \right)}} = 1 - {{S}_{b}}_{N}^{{\left( {k + 1} \right)}} = 0. \\ \end{gathered} $

РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ

Алгоритм (22)–(34) реализован численно на языке C#. Моделирование проводилось для физических параметров задачи, сведенных в табл. 1, 2 . Количество пространственных узлов конечно-разностной сетки принято равным 400, а число Куранта–Фридрихса–Леви равно 0.3.

Таблица 1.

Геометрические и теплофизические характеристики системы в целом

Параметр	Значение
Радиус замораживающей колонки, м	0.073
Радиус расчетной области, м	20
Пористость массива	0.2
Начальная температура породного массива, °С	10
Температура ликвидуса поровой воды, °С	0
Температура солидуса поровой воды, °С	–2
Начальная объемная концентрация рассола в порах, м³/м³	1
Начальная объемная концентрация соли в рассоле, м³/м³	0.05
Объемная концентрация насыщенного раствора соли, м³/м³	0.212
Коэффициент понижения температуры фазового перехода (для раствора соли NaCl), °С	66.7
Коэффициент диффузии соли в рассоле, мм²/с	0.01
Коэффициент искусственной диффузии, мм²/с	0.00025
Температура хладоносителя в замораживающей колонке, °С	–20
Коэффициент теплоотдачи на границе с замораживающей колонкой, Вт/(м² °С)	50
Время моделирования, сут.	200

Таблица 2.

Теплофизические свойства фаз

Параметр	Значения
Параметр	сухой скелет	вода	лед	соль (NaCl)
Плотность, кг/м³	2000	1000	912	2160
Удельная теплоемкость, Дж/(кг °С)	800	4190	2000	855
Теплопроводность, Вт/(м °С)	2	0.5	2.23	0.58

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

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

Рис. 2.

Зависимости температуры (а) и объемных концентраций рассола и льда (б) от пространственной координаты: 1 – лед, 2 – рассол.

Рис. 3.

Зависимости объемной концентрации растворенной в рассоле соли (а), а также масс растворенной и выпавшей в осадок соли (б) от пространственной координаты: 1 – раствор, 2 – осадок.

Рис. 4.

Зависимости масс растворенной (1) и выпавшей в осадок (2) соли (а), а также объемных концентраций льда (1) и рассола (2) (б) от пространственной координаты; коэффициент диффузии соли равен 10^–5 мм²/с.

На рис. 2 можно видеть снижение температур ликвидуса (с 0 до –3.6°С) и солидуса (с –2 до –5.9°С). Это снижение приводит к существенному (около 0.5 м) смещению положения границ полностью замороженного породного массива. Такие значительные изменения граничных температур фазового перехода связаны с высокой минерализацией породного массива (5%-ный рассол NaCl). На практике минерализация пород, как правило, ниже. Так, в условиях Старобинского месторождения подземные воды практически пресные, лишь в верхней части глинисто-мергелистой толщи их минерализация достигает 76 г/л (3.5%). При этом, по данным [1], в отдельных случаях встречаются незначительные скопления рассолов с минерализацией до 400 г/л (18.5%). Интересно отметить, что вследствие неоднородного распределения концентрации соли в рассоле (рис. 3а), ширина температурного интервала фазового перехода увеличилась с 2 до 2.3°С. При более низком коэффициенте диффузии соли в рассоле ширина фазового перехода еще сильнее увеличится (рис. 4б).

В области от 0.7 до 0.9 м (рис. 2б), соответствующей зоне фазового перехода, происходит резкое изменение концентраций льда и рассола. При этом в зоне льда (r < 0.7 м) лед занимает не все поровое пространство – часть его остается занятой солью, выпавшей в осадок (рис. 3б). На рис. 3б также можно отметить наличие “зуба” на зависимостях как массы растворенной соли, так и массы выпавшей в осадок соли – это говорит о том, что в зоне фазового перехода увеличение концентрации соли (рис. 3а) происходит не только за счет частичного замерзания воды в рассоле, но и за счет диффузионного переноса соли из более замороженной области в менее замороженную. В случае близкой к нулю диффузии соли данные зависимости, по сути, упрощаются до вида классической ступенчатой функции Хэвисайда (рис. 4а).

ЗАКЛЮЧЕНИЕ

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

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

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

Работа выполнена при финансовой поддержке Минобрнауки Пермского края в рамках научного проекта № С-26/563.

Список литературы

Ольховиков Ю.П. Крепь капитальных выработок калийных и соляных рудников. М.: Недра, 1984. 238 с.
Yong R.N., Cheung C.H., Sheeran D.E. Prediction of Salt Influence on Unfrozen Water Content in Frozen Soils // Developments in Geotechnical Engineering. 1979. V. 26. P. 137.
Bing H., Ma W. Laboratory Investigation of the Freezing Point of Saline Soil // Cold Regions Sci. Technol. 2011. V. 67. № 1–2. P. 79.
Wan X., Lai Y., Wang C. Experimental Study on the Freezing Temperatures of Saline Silty Soils // Permafrost and Periglacial Processes. 2015. V. 26. № 2. P. 175.
Banin A., Anderson D.M. Effects of Salt Concentration Changes During Freezing on the Unfrozen Water Content of Porous Materials // Water Resources Research. 1974. V. 10. № 1. P. 124.
Frivik P.E. State-of-the-art Report. Ground Freezing: Thermal Properties, Modelling of Processes and Thermal Design // Eng. Geology. 1981. V. 18. № 1–4. P. 115.
Lucas T., Chourot J.-M., Bohuon Ph., Flick D. Freezing of a Porous Medium in Contact with a Concentrated Aqueous Freezant: Numerical Modelling of Coupled Heat and Mass Transport // Int. J. Heat Mass Transfer. 2001. V. 44. № 11. P. 2093.
Plekhov O. et al. The Effect of Cryogenic Suction on the Monitoring Data of Ice Barrier Formation in a Porous Water-saturated Soil // Proc. Structural Integrity. 2019. V. 17. P. 602.
Rouabhi A., Jahangir E., Tounsi H. Modeling Heat and Mass Transfer During Ground Freezing Taking into Account the Salinity of the Saturating Fluid // Int. J. Heat Mass Transfer. 2018. V. 120. P. 523.
Tounsi H., Rouabhi A., Jahangir E. Thermo-hydro-mechanical Modeling of Artificial Ground Freezing Taking into Account the Salinity of the Saturating Fluid // Computers and Geotechnics. 2020. V. 119. 103382.
Zhang X. et al. Numerical Study on the Multifield Mathematical Coupled Model of Hydraulic-thermal-salt-mechanical in Saturated Freezing Saline Soil // Int. J. Geomechanics. 2018. V. 18. № 7. 04018064.
Zhang J. et al. Study on the Mechanism of Crystallization Deformation of Sulfate Saline Soil During the Unidirectional Freezing Process // Permafrost and Periglacial Processes. 2020. V. 31. № 1. P. 102.
Wu D., Lai Y., Zhang M. Thermo-hydro-salt-mecha-nical Coupled Model for Saturated Porous Media Based on Crystallization Kinetics // Cold Regions Sci. Technol. 2017. V. 133. P. 94.
Wu D., Zhou X., Jiang X. Water and Salt Migration with Phase Change in Saline Soil During Freezing and Thawing Processes // Groundwater. 2018. V. 56. № 5. P. 742.
Васильев В.И., Максимов А.М., Петров Е.Е., Цыпкин Г.Г. Математическая модель замерзания-таяния засоленного мерзлого грунта // ПМТФ. 1995. Т. 36. № 5. С. 57.
Galushkin Y.I., Sitar K.A., Frolov S.V. Basin Modelling of Temperature and Heat Flow Distributions and Permafrost Evolution, Urengoy and Kuyumbinskaya Areas, Siberia // Permafrost and Periglacial Processes. 2013. V. 24. № 4. P. 268.
Цытович Н.А. Механика мерзлых грунтов. М.: URSS, 2009.
Де Гроот С.Р. Термодинамика необратимых процессов. М.: ГИТТЛ, 1956. 277 с.
Mortimer R.G., Eyring H. Elementary Transition State Theory of the Soret and Dufour Effects // Proc. National Academy of Sciences. 1980. V. 77. № 4. P. 1728.
Jochem M., Körber C. A Numerical Solution of the Coupled Heat and Mass Transfer Problem of Non-planar Solidification and Melting of Aqueous Solutions // Wärme- und Stoffübertragung. 1993. Bd. 28. № 4. S. 195.
Panteleev I.A. et al. Numerical Simulation of Artificial Ground Freezing in a Fluid-saturated Rock Mass with Account for Filtration and Mechanical Processes // Sciences in Cold and Arid Regions. 2018. V. 9. № 4. P. 363.
Ma G.-Y., Du M.-J., Li D. Numerical Calculation for Temperature Coupled with Moisture and Stress of Soil Around Buried Pipeline in Permafrost Regions // J. China University of Petroleum (Edition of Natural Science). 2011. V. 35. № 3. P. 108.
Ma J., Wang X. Natural Convection and Its Fractal for Liquid Freezing in a Vertical Cavity Filled with Porous Medium // Heat Transfer–Asian Research: Co-sponsored by the Society of Chemical Engineers of Japan and the Heat Transfer Division of ASME. 1999. V. 28. № 3 P. 165.
Semin M.A. et al. Natural Convection in Water-Saturated Rock Mass under Artificial Freezing // J. Mining Sci. 2020. V. 56. № 2. P. 297.
Yong-ji X. Research on Soil Moisture and Thermophoresis in the Course of Seasonal Freeze Thawing [J] // Taiyuan Science & Technology. 2008. V. 3.
Алексютина Д.М., Мотенко Р.Г. Оценка влияния засоления и содержания органического вещества в мерзлых породах на западном побережье Байдарацкой губы, их теплофизические свойства и фазовый состав влаги // Вестн. Моск. ун-та. Сер. 4. Геология. 2016. № 2.
Семин М.А., Богомягков А.В., Левин Л.Ю. Теоретический анализ динамики ледопородного ограждения при переходе на пассивный режим замораживания // Зап. Горн. ин-та. 2020. Т. 243.
Semin M., Levin L. Numerical Simulation of Frozen Wall Formation in Water-saturated Rock Mass by Solving the Darcy-Stefan Problem // Frattura ed Integrità Strutturale. 2019. V. 13. № 49. P. 167.
Tsypkin G.G. Water–Ice Phase Transition in Unsaturated Soil in the Presence of Capillary Pressure // Fluid Dynamics. 2019. V. 54. № 5. P. 681.
Van Genuchten M.T. A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils // Soil Sci. Soc. Amer. J. 1980. V. 44. № 5. P. 892.
Leverett M.C. et al. Capillary Behavior in Porous Solids // Trans. AIME. 1941. V. 142. № 1. P. 152.
Anderson D.M., Tice A.R., McKim H.L. The Unfrozen Water and the Apparent Specific Heat Capacity of Frozen Soils // 7nd Int. Conf. on Permafrost. Yakutsk, USSR. North American Contribution. 1973. P. 289.
Côté J., Konrad J.M. A Generalized Thermal Conductivity Model for Soils and Construction Materials // Canad. Geotech. J. 2005. V. 42. № 2. P. 443.
Kantzas A., Bryan J., Taheri S. Fundamentals of Fluid Flow in Porous Media // Pore Size Distribution. 2012. V. 1. P. 1.
Самарский A.А., Моисеенко Б.Д. Экономичная схема сквозного счета для многомерной задачи Стефана // ЖВМиМФ. 1965. Т. 5. № 5. С. 816.
Трупак Н.Г. Замораживание грунтов в подземном строительстве. М.: Недра, 1974.
Курант Р., Фридрихс К., Леви Г. О разностных уравнениях математической физики // Усп. матем. наук. 1941. № 8. С. 125.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Теплофизика высоких температур

Архивы выпусков Информация о журнале Отправить рукопись в журнал