1. На главную
  2. Электронные версии
  3. Лёд и Снег
  4. T. 63, № 1, 2023

Лёд и Снег, 2023, T. 63, № 1, стр. 130-140

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

С. В. Попов 1, 2, 3, *

1 Полярная морская геологоразведочная экспедиция
Санкт-Петербург, Россия

2 Санкт-Петербургский государственный университет
Санкт-Петербург, Россия

3 Институт мерзлотоведения им. П.И. Мельникова СО РАН
Якутск, Россия

* E-mail: spopov67@yandex.ru

Поступила в редакцию 21.07.2022
После доработки 19.09.2022
Принята к публикации 06.03.2023

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

Аннотация

Представлено численное решение одномерной задачи Стефана с двумя фазовыми границами в виде конечно-разностных схем, реализованных на неравномерной сетке. Уравнения записаны в наиболее общей форме, то есть включают в себя не только кондуктивный, но также конвективный и диссипативный члены. В качестве примера выполнено оценочное моделирование процесса замерзания трещины в леднике, заполненной водой. Получено, что для ледников с температурой ниже –5°C время замерзания 30-сантиметровой трещины составляет менее трёх месяцев.

Ключевые слова: математическое моделирование, численное решение, задача Стефана, конечно-разностные схемы, ледниковые трещины, Антарктида

ВВЕДЕНИЕ

Основой познания окружающего мира считаются не только результаты полевых исследований, но и математическое моделирование. Именно оно позволяет наиболее эффективно применять полученные данные. Современные модели описания гляциальных и субгляциальных процессов, а также процессов, протекающих в подледниковых водоемах от начала их формирования до завершения всего цикла существования, очень сложны. Как правило, они основаны на фундаментальных законах природы и описываются уравнениями математической физики, например, модель тепломассопереноса в ледниках (Huybrechts, 1992; Greve, 1997; Greve, Blatter, 2009; Рыбак, Рыбак, 2010), или движения воды в подледниковом озере Восток (Thoma et al., 2007; Казко и др., 2012). При этом формируется целый класс процессов, который может быть описан решением задачи Стефана о движении фазовой границы. Они весьма актуальны, применительно к гляциальным и субгляциальным гидрологическим процессам, поскольку именно задача Стефана (если точнее, то определенный тип граничных условий) описывает изменение конфигурации поверхности льда на границе с водой.

Теоретически, описывать эти процессы можно по-разному. Например, Дж. Най, моделируя прорывы подледниковых водоёмов, применял закон сохранения энергии и эмпирические соотношения (Nye, 1976). Сходным образом поступил и Ф. Паттин (Pattyn, 2003). Это было обосновано, прежде всего тем, что в указанных работах описывались процессы, протекающие под ледником значительной мощности. Поскольку лёд – это пластичное вещество, то конфигурация границы раздела льда и воды, в соответствии с законом Глена (Патерсон, 1984), определяется также вышележащим весом. Расчёты автора показывают, что при мощности ледника 1 км канал сечением 1 м2 полностью затекает примерно за 2.5 часа, поэтому моделировать конфигурацию границы между льдом и водой особого смысла не имеет. С другой стороны, имеется достаточно много задач, для которых это не только целесообразно, но и важно. В частности, одной из них считается описание процесса прорыва приледниковых водоемов. Поскольку мощность ледника или снежника относительно невелика, то и его деформацией за счёт растекания можно пренебречь. В этом случае конфигурация образующегося тоннеля описывается решением задачи Стефана. В целом она хорошо разработана, и существуют множество монографий по этой тематике. Однако есть известная разница между теоретическими построениями и их практической реализацией, которая не всегда очевидна, особенно при решении конкретных и сложных задач.

Исходя из общих соображений, применение моделирования поведения двух фазовых границ вполне достаточно для рассмотрения абсолютного большинства даже самых сложных случаев, применительно к гляциальным и субгляциальным процессам от формирования и развития подледниковых водоёмов, их прорывов, развития таликов в многолетней мерзлоте, или донного льда, до моделирования подледниковых водоёмов полярных шапок Марса. Основой расчетов для них служит решение задачи Стефана. Аналитически она может быть реализована лишь в весьма ограниченном количестве случаев. Поэтому целесообразно рассмотреть численное решение для произвольных начальных и граничных условий. Опыт автора показывает, что имеет также смысл применять неравномерную сетку, которая позволит существенно повысить точность построений на границах раздела сред. Особенно это полезно, когда моделируемый процесс предполагает формирование или исчезновение слоев. Именно эта задача и решается в рамках настоящей работы. В ней предлагается вариант конечно-разностной схемы, который может быть впоследствии применен во всевозможных задачах моделирования. Распространение на двумерный или трёхмерный случай достигается применением известных способов, в частности, метода переменных направлений (Самарский, 1977; Кузнецов, Шеремет, 2007; Кольцова и др., 2020). Решение иллюстрируется на примере весьма упрощенного одномерного процесса замерзания воды между двумя стенками трещины в леднике.

ПОСТАНОВКА ЗАДАЧИ

Для определённости смоделируем упрощенный процесс замерзания трещины исключительно как иллюстрацию к решению задачи Стефана с двумя фазовыми границами. Предположим, что существует безграничный однородный и изотропный ледник. В некоторый момент времени в нем образуется трещина бесконечной глубины, которая заполняется талой ледниковой водой и начинает постепенно замерзать. Реальным примером может служить процесс залечивания ледниковых трещин, сформированных на так называемом голубом льду, то есть на леднике, где в силу различных причин не накапливается снежный покров. В течение относительно тёплого антарктического лета на его поверхности происходит интенсивное таяние. Талые воды стекают в трещины и впоследствии замерзают, залечивая их. Подобная ситуация наблюдалась автором в районе станции Мирный на леднике недалеко от сопки Ветров (рис. 1, а). В ходе недавних инженерных изысканий, выполненных на посадочной площадке станции Новолазаревская (Попов и др., 2022), выявлены линейные дайкообразные структуры голубого цвета, которые контрастно выделяются на фоне более светлого массива льда (см. рис. 1, б). Они считаются трещинами, которые были залечены аналогичным образом.

Рис. 1.

Фотографии залеченных трещин: а – частично замерзшие трещины с талой водой в приповерхностной части ледника в районе сопки Ветров, станция Мирный (фото автора, февраль 2017 г.); б – полностью замерзшие ледниковые трещины на взлётно-посадочной полосе станции Новолазаревская (фото А.С. Борониной, ноябрь 2021 г.).

Fig. 1. Photos of crevasses: а – partially frozen crevasse with melt water in the near-surface part on the glacier in the Sopka Vetrov area, Mirny Station (photo by the author, February 2017); б – completely frozen crevasse on the Novo Runway (photo by A. Boronina, November 2021).

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

Таким образом, модель включает в себя три области: Ω1 (лёд), Ω2 (вода) и Ω3 (лёд), между которыми формируются подвижные границы Γ12 = = Ω1 ∩ Ω2 и Γ23 = Ω2 ∩ Ω3. Ось абсцисс направлена вправо и совместима с началом координат с левой границей Γ1 области Ω1 (рис. 2). Все три среды характеризуются разными плотностями ρ, а также коэффициентами теплопроводности λ и теплоемкости c. Нижними индексами 1, 2 и 3 обозначены номера областей. Если a – это коэффициент температуропроводности, a = λ/(с · ρ), то для одномерного случая предлагается система уравнений, описывающая изменение температуры θ в средах в зависимости от времени t и расстояния x вдоль оси абсцисс (Смирнов, 1974; Тихонов, Самарский, 1977):

(1)
$\begin{gathered} \frac{{\partial {{{{\theta }}}_{k}}}}{{\partial t}} = {{a}_{k}}\frac{{{{\partial }^{2}}{{{{\theta }}}_{k}}}}{{\partial {{x}^{2}}}} + {{v}_{k}}\frac{{\partial {{{{\theta }}}_{k}}}}{{\partial x}} + {{\Phi }_{k}}, \\ {{x}_{k}} \in {{\Omega }_{k}},\quad k = 1,2,3, \\ \end{gathered} $
где ${{{v}}_{k}}$ – вертикальные скорости движения сред, а Φk – диссипативный член, связанный с выделением тепла при трении слоев сред друг о друга при перемещении. Отмечено, что ρ, λ и c не считаются константами (Патерсон, 1984), но они неизменны в пределах своих сред. Рассмотрен упрощенный пример моделирования, который не предполагает физического перемещения вещества сред. В этом случае следовало бы убрать второй и третий члены в правой части уравнения (1), однако они оставлены для того, чтобы получить наиболее общее решение задачи Стефана.

Рис. 2.

Иллюстрация к математической модели.

Fig. 2. Illustration to the mathematical model.

На левой и правой границах Γ1 и Γ2 (см. рис. 2) заданы условия I рода (Дирихле), соответствующие средней температуре в средах ${{{{\bar {\theta }}}}_{k}}$ и считающиеся начальными условиями:

(2)
$\begin{gathered} {{\left. {{{{{\theta }}}_{1}}(t)} \right|}_{{{{\Gamma }_{1}}}}} = {{{{{\bar {\theta }}}}}_{1}},\quad {{\left. {{{{{\theta }}}_{{\text{3}}}}(t)} \right|}_{{{{\Gamma }_{3}}}}} = {{{{{\bar {\theta }}}}}_{3}},\quad {{\left. {{{{{\theta }}}_{1}}(x)} \right|}_{{t = 0}}} = {{{{{\bar {\theta }}}}}_{1}}, \\ {{\left. {{{{{\theta }}}_{2}}(x)} \right|}_{{t = 0}}} = {{{{{\bar {\theta }}}}}_{2}},\quad {{\left. {{{{{\theta }}}_{3}}(x)} \right|}_{{t = 0}}} = {{{{{\bar {\theta }}}}}_{3}}. \\ \end{gathered} $

Это вполне логично, исходя из постановки задачи: сначала существовал безграничный ледник, затем в некоторый момент времени t = t0 произошел его раскол и сформировалась трещина. Согласно условиям решаемой задачи, на внутренних подвижных границах Γ12 и Γ23 задаются условия Стефана:

(3)
$\begin{gathered} {{\left. {{{{{\lambda }}}_{1}}\frac{{\partial {{{{\theta }}}_{1}}}}{{\partial x}}} \right|}_{{{{\Gamma }_{{12}}}}}} - {{\left. {{{{{\lambda }}}_{2}}\frac{{\partial {{{{\theta }}}_{2}}}}{{\partial x}}} \right|}_{{{{\Gamma }_{{12}}}}}} = {{q}_{{F12}}}{{{{\rho }}}_{2}}\frac{{d{{s}_{1}}}}{{dt}}, \\ {{\left. {{{{{\lambda }}}_{2}}\frac{{\partial {{{{\theta }}}_{2}}}}{{\partial x}}} \right|}_{{{{\Gamma }_{{23}}}}}} - {{\left. {{{{{\lambda }}}_{3}}\frac{{\partial {{{{\theta }}}_{3}}}}{{\partial x}}} \right|}_{{{{\Gamma }_{{23}}}}}} = {{q}_{{F23}}}{{{{\rho }}}_{2}}\frac{{d{{s}_{2}}}}{{dt}}, \\ \end{gathered} $
где qF12 и qF23 – удельная теплота плавления, требуемая для перехода вещества из одного агрегатного состояния в другое на границах раздела сред Γ12 и Γ23 (в общем случае они могут отличаться), а s1 и s2 – их положение на оси абсцисс.

Поскольку плотности агрегатных состояний различны, смещения границ также приводят и к изменению размеров сред. Для примера рассмотрена граница Γ12. Если она сдвинулась из области Ω1 в Ω2 на ds1, то фактически в силу закона сохранения массы это означает, что область Ω1 увеличится на эту величину, но Ω2 уменьшится на ds1 + Δs1:

$\Delta {{s}_{1}} = \frac{{d{{s}_{1}}}}{{{{{{\rho }}}_{2}}}}\left( {{{{{\rho }}}_{1}} - {{{{\rho }}}_{2}}} \right).$

В противоположном случае сдвиг на s1 уменьшит область Ω1 на эту величину, но Ω2 увеличится на ds1 + Δs1. Аналогична ситуация и со смещением границы Γ23. Если переходить в практическую плоскость, то эффективный учёт вышеизложенного достигается путем изменения размеров соответствующих сред на величины Δs1 и Δs2.

Сформулированная выше краевая задачу в виде уравнений (1), (2), (3) решена численно методом спрямления фронтов. Его достоинство заключается в том, что он позволяет рассматривать не только одномерные, но и многомерные задачи Стефана с произвольным количеством границ (Самарский, 1977; Кузнецов, Шеремет, 2007; Кольцова и др., 2020). Произведён переход из системы координат (t, x) в новую (τ, η) согласно правилу

(4)
${{\eta }}(x) = \left\{ \begin{gathered} \frac{x}{{{{s}_{1}}({{\tau }})}},\quad {{\Omega }_{1}},\;0 \leqslant x \leqslant {{s}_{1}}({{\tau }}) \hfill \\ \frac{{x - {{s}_{2}}({{\tau }})}}{{{{s}_{1}}({{\tau }}) - {{s}_{2}}({{\tau }})}},\quad {{\Omega }_{2}},\;{{s}_{1}}({{\tau }}) \leqslant x \leqslant {{s}_{2}}({{\tau }}) \hfill \\ \frac{{x - {{s}_{2}}({{\tau }})}}{{L - {{s}_{2}}({{\tau }})}},\quad {{\Omega }_{3}},\;{{s}_{2}}({{\tau }}) \leqslant x \leqslant L \hfill \\ \end{gathered} \right.,\quad {{\tau }} = t.$

Если требуется рассмотреть большее количество границ раздела сред, то η(x) записывается аналогичным образом. В нашем случае для области Ω1: η(0) = 0 и η(s1) = 1, для Ω2: η(s1) = 1 и η(s2) = = 0 и для Ω3: η(s2) = 0 и η(L) = 1, где L – это положение границы Γ3 (правого края Ω3) на оси абсцисс (см. рис. 2). Соотношения (1), (2) и (3) должны быть преобразованы с учётом (4). Для этого записаны частные производные, входящие в указанные соотношения в новых координатах для каждой области. Производные по времени для них выглядят следующим образом:

$\frac{{\partial {{{{\theta }}}_{1}}}}{{\partial t}} = \frac{{\partial {{{{\theta }}}_{1}}}}{{\partial {{\tau }}}}\frac{{\partial {{\tau }}}}{{\partial t}} + \frac{{\partial {{{{\theta }}}_{1}}}}{{\partial {{\eta }}}}\frac{{\partial {{\eta }}}}{{\partial {{s}_{1}}}}\frac{{\partial {{s}_{1}}}}{{\partial {{\tau }}}}\frac{{\partial {{\tau }}}}{{\partial t}},\quad {\text{при}}\quad {{\eta }} \in {{\Omega }_{1}},$
$\begin{gathered} \frac{{\partial {{{{\theta }}}_{2}}}}{{\partial t}} = \frac{{\partial {{{{\theta }}}_{2}}}}{{\partial {{\tau }}}}\frac{{\partial {{\tau }}}}{{\partial t}} + \frac{{\partial {{{{\theta }}}_{2}}}}{{\partial {{\eta }}}}\frac{{\partial {{\eta }}}}{{\partial {{s}_{1}}}}\frac{{\partial {{s}_{1}}}}{{\partial {{\tau }}}}\frac{{\partial {{\tau }}}}{{\partial t}} + \frac{{\partial {{{{\theta }}}_{2}}}}{{\partial {{\eta }}}}\frac{{\partial {{\eta }}}}{{\partial {{s}_{2}}}}\frac{{\partial {{s}_{2}}}}{{\partial {{\tau }}}}\frac{{\partial {{\tau }}}}{{\partial t}}, \\ {\text{при}}\quad {{\eta }} \in {{\Omega }_{2}}, \\ \end{gathered} $
$\frac{{\partial {{{{\theta }}}_{3}}}}{{\partial t}} = \frac{{\partial {{{{\theta }}}_{3}}}}{{\partial {{\tau }}}}\frac{{\partial {{\tau }}}}{{\partial t}} + \frac{{\partial {{{{\theta }}}_{3}}}}{{\partial {{\eta }}}}\frac{{\partial {{\eta }}}}{{\partial {{s}_{2}}}}\frac{{\partial {{s}_{2}}}}{{\partial {{\tau }}}}\frac{{\partial {{\tau }}}}{{\partial t}},\quad {\text{при}}\quad {{\eta }} \in {{\Omega }_{3}}.$

Поскольку

$\begin{gathered} \frac{{\partial {{{{\theta }}}_{k}}}}{{\partial x}} = \frac{{\partial {{{{\theta }}}_{k}}}}{{\partial {{\eta }}}}\frac{{\partial {{\eta }}}}{{\partial x}}, \\ {\text{и}}\quad \frac{{{{\partial }^{2}}{{{{\theta }}}_{k}}}}{{\partial {{x}^{2}}}} = \frac{\partial }{{\partial x}}\left( {\frac{{\partial {{{{\theta }}}_{k}}}}{{\partial x}}} \right) = \frac{\partial }{{\partial {{\eta }}}}\left( {\frac{{\partial {{{{\theta }}}_{k}}}}{{\partial {{\eta }}}}\frac{{\partial {{\eta }}}}{{\partial x}}} \right)\frac{{\partial {{\eta }}}}{{\partial x}}, \\ k = 1,2,3, \\ \end{gathered} $
то окончательно производные по плановым координатам примут вид:
$\frac{{\partial {{\eta }}}}{{\partial {{s}_{1}}}} = \frac{{ - {{\eta }}}}{{{{s}_{1}}(\tau )}}\quad {\text{и}}\quad \frac{{\partial {{\eta }}}}{{\partial x}} = \frac{1}{{{{s}_{1}}({{\tau }})}},\quad {\text{при}}\quad {{\eta }} \in {{\Omega }_{1}},$
$\begin{gathered} \frac{{\partial {{\eta }}}}{{\partial {{s}_{1}}}} = \frac{{ - {{\eta }}}}{{{{s}_{1}}({{\tau }}) - {{s}_{2}}({{\tau }})}},\quad \frac{{\partial {{\eta }}}}{{\partial {{s}_{2}}}} = \frac{{{{\eta }} - 1}}{{{{s}_{1}}({{\tau }}) - {{s}_{2}}({{\tau }})}} \\ {\text{и}}\quad \frac{{\partial {{\eta }}}}{{\partial x}} = \frac{1}{{{{s}_{1}}({{\tau }}) - {{s}_{2}}({{\tau }})}},\quad {\text{при}}\quad {{\eta }} \in {{\Omega }_{2}}, \\ \end{gathered} $
$\begin{gathered} \frac{{\partial {{\eta }}}}{{\partial {{s}_{2}}}} = \frac{{{{\eta }} - 1}}{{L - {{s}_{2}}({{\tau }})}}\quad {\text{и}}\quad \frac{{\partial {{\eta }}}}{{\partial x}} = \frac{1}{{L - {{s}_{2}}({{\tau }})}}, \\ {\text{при}}\quad {{\eta }} \in {{\Omega }_{3}}, \\ \end{gathered} $
с учётом того, что $\frac{{\partial {{\tau }}}}{{\partial t}} = 1$ уравнения теплопроводности (1) могут быть записаны следующим образом:

(5)
$\frac{{\partial {{{{\theta }}}_{1}}}}{{\partial {{\tau }}}} = \frac{{{{a}_{1}}}}{{s_{1}^{2}}}\frac{{{{\partial }^{2}}{{{{\theta }}}_{1}}}}{{\partial {{{{\eta }}}^{2}}}} + \frac{{{\eta }}}{{{{s}_{1}}}}\frac{{\partial {{{{\theta }}}_{1}}}}{{\partial {{\eta }}}}\frac{{\partial {{s}_{1}}}}{{\partial {{\tau }}}} + \frac{{{{v}_{1}}}}{{{{s}_{1}}}}\frac{{\partial {{{{\theta }}}_{1}}}}{{\partial {{\eta }}}} + {{\Phi }_{1}},\quad {{\eta }} \in {{\Omega }_{1}},$
(6)
$\begin{gathered} \frac{{\partial {{{{\theta }}}_{2}}}}{{\partial {{\tau }}}} = \frac{{{{a}_{2}}}}{{{{{\left( {{{s}_{1}} - {{s}_{2}}} \right)}}^{2}}}}\frac{{{{\partial }^{2}}{{{{\theta }}}_{2}}}}{{\partial {{{{\eta }}}^{2}}}} + \frac{1}{{{{s}_{1}} - {{s}_{2}}}} \times \\ \times \;\frac{{\partial {{{{\theta }}}_{1}}}}{{\partial {{\eta }}}}\left[ {{{\eta }}\frac{{\partial {{s}_{1}}}}{{\partial {{\tau }}}} + \left( {1 - {{\eta }}} \right)\frac{{\partial {{s}_{2}}}}{{\partial {{\tau }}}} + {{v}_{2}}} \right] + {{\Phi }_{2}}, \\ {{\eta }} \in {{\Omega }_{2}}, \\ \end{gathered} $
(7)
$\begin{gathered} \frac{{\partial {{{{\theta }}}_{3}}}}{{\partial {{\tau }}}} = \frac{{{{a}_{3}}}}{{{{{\left( {L - {{s}_{2}}} \right)}}^{2}}}}\frac{{{{\partial }^{2}}{{{{\theta }}}_{2}}}}{{\partial {{{{\eta }}}^{2}}}} + \frac{1}{{L - {{s}_{2}}}} \times \\ \times \;\frac{{\partial {{{{\theta }}}_{1}}}}{{\partial {{\eta }}}}\left[ {\left( {1 - {{\eta }}} \right)\frac{{\partial {{s}_{2}}}}{{\partial {{\tau }}}} + {{v}_{3}}} \right] + {{\Phi }_{3}}, \\ {{\eta }} \in {{\Omega }_{3}}. \\ \end{gathered} $

Начальные условия и условия на границах Γ1 и Γ3 сохранятся, а условия Стефана примут новый вид

(8)
${{\left. {\frac{{{{{{\lambda }}}_{1}}}}{{{{s}_{1}}}}\frac{{\partial {{{{\theta }}}_{1}}}}{{\partial {{\eta }}}}} \right|}_{{{{\Gamma }_{{12}}}}}} - {{\left. {\frac{{{{{{\lambda }}}_{2}}}}{{{{s}_{1}} - {{s}_{2}}}}\frac{{\partial {{{{\theta }}}_{2}}}}{{\partial {{\eta }}}}} \right|}_{{{{\Gamma }_{{12}}}}}} = {{q}_{{F12}}}{{{{\rho }}}_{2}}\frac{{d{{s}_{1}}}}{{d{{\tau }}}},$
(9)
${{\left. {\frac{{{{{{\lambda }}}_{2}}}}{{{{s}_{1}} - {{s}_{2}}}}\frac{{\partial {{{{\theta }}}_{2}}}}{{\partial {{\eta }}}}} \right|}_{{{{\Gamma }_{{23}}}}}} - {{\left. {\frac{{{{{{\lambda }}}_{3}}}}{{L - {{s}_{2}}}}\frac{{\partial {{{{\theta }}}_{3}}}}{{\partial {{\eta }}}}} \right|}_{{{{\Gamma }_{{23}}}}}} = {{q}_{{F23}}}{{{{\rho }}}_{2}}\frac{{d{{s}_{2}}}}{{d{{\tau }}}}.$

Таким образом, уравнения (5), (6) и (7), а также начальные и граничные условия (2), (8) и (9) представляют собой краевую задачу, которая может быть решена численно.

МЕТОД ЧИСЛЕННОГО РЕШЕНИЯ

Известные работы, в частности (Самарский, 1977; Кузнецов, Шеремет, 2007; Кольцова и др., 2020), а также недавнее исследование автора показывают, что результаты моделирования, выполненного с применением основных конечно-разностных схем (явной, неявной и Кранка–Николсона), практически не отличаются друг от друга. Естественно, речь идёт о корректном решении исходя из условия устойчивости Куранта. Отмечено, что параметры конечно-разностной схемы, при которых число Куранта ζ ≈ 0.1 (ζ ≡ 2Δt/Δx2), где Δt и Δx – шаг по времени и по расстоянию) считаются универсальными, и дают наименьшую погрешность при различных соотношениях Δt и Δx.

Решение сформулированной выше задачи реализовано в виде четырехточечной неявной конечно-разностной схемы (Самарский, 1977; Кузнецов, Шеремет, 2007; Кольцова и др., 2020). Поскольку границы Γ12 и Γ23 двигаются, и в предельном случае центральная область Ω2 полностью вырождается, то для повышения точности моделирования целесообразно применять неравномерную сетку (Самарский, 1977; Кузнецов, Шеремет, 2007; Кольцова и др., 2020). Её формирование различно, поэтому в общем случае при построении конечно-разностной схемы не нужно привязываться к какому-то конкретному её виду. Однако не следует при этом забывать про критерий устойчивости Куранта. Пусть неравномерная сетка ξ(η) формируется по некоторому закону таким образом, что ξ(0) = 0 и ξ(1) = 1.

Для удобства последующей записи считается, что все три области разделены на одинаковое количество отрезков, равное N. Их границами считаются N + 1 точек, нумеруемые от 0 до N. Традиционно верхним индексом n обозначается номер временнóго слоя, а нижним j – номер точки по оси абсцисс. В итоге задача решается методом прогонки (Самарский, 1977; Кузнецов, Шеремет, 2007; Кольцова и др., 2020), а конечно-разностная схема приведена к виду:

${{A}_{j}}{{\theta }}_{{j + 1}}^{{n + 1}} - {{B}_{j}}{{\theta }}_{j}^{{n + 1}} + {{C}_{j}}{{\theta }}_{{j - 1}}^{{n + 1}} = {{F}_{j}},$
где ${{A}_{j}}$, ${{B}_{j}}$, ${{C}_{j}}$ и ${{F}_{j}}$ – некоторые коэффициенты. Пусть $\Delta {{{{\xi }}}_{j}} \equiv {{{{\xi }}}_{{j + 1}}} - {{{{\xi }}}_{j}}$, тогда выражения (7), (8) и (9) примут общий вид
(10)
$\begin{gathered} {{\theta }}_{{k,j + 1}}^{{n + 1}}\left[ {{{{{\varphi }}}_{{k,j}}}\Delta {{{{\xi }}}_{{j - 1}}} + {{{{\psi }}}_{{k,j}}}} \right] - {{\theta }}_{{k,j}}^{{n + 1}}\left[ {{{{{\varphi }}}_{{k,j}}}{{{{\xi }}}_{j}} + 1} \right] + \\ + \;{{\theta }}_{{j - 1}}^{{n + 1}}\left[ {{{{{\varphi }}}_{{k,j}}}\Delta {{{{\xi }}}_{j}} - {{{{\psi }}}_{{k,j}}}} \right] = F_{{k,j}}^{n}, \\ k = 1,2,3, \\ \end{gathered} $
(11)
$\begin{gathered} {{{{\varphi }}}_{{k,j}}} \equiv \frac{{2{{a}_{k}}\Delta {{\tau }}}}{{{{{(T_{k}^{n})}}^{2}}(\Delta {{{{\xi }}}_{{j - 1}}} + \Delta {{{{\xi }}}_{j}})\Delta {{{{\xi }}}_{{j - 1}}}\Delta {{{{\xi }}}_{j}}}}, \\ {{{{\psi }}}_{{k,j}}} \equiv \frac{{G_{{k,j}}^{n} + {{v}_{k}}\Delta {{\tau }}}}{{T_{k}^{n}(\Delta {{{{\xi }}}_{{j - 1}}} + \Delta {{{{\xi }}}_{j}})}},\quad F_{{k,j}}^{n} \equiv - ({{\theta }}_{j}^{n} + {{\Phi }_{k}}\Delta {{\tau }}), \\ \end{gathered} $
где $T_{k}^{n}$ – мощность k-го слоя, а $G_{{k,j}}^{n}$ – некоторый параметр:

(12)
$\begin{gathered} T_{k}^{n} = \left\{ \begin{gathered} s_{1}^{n},\quad k = 1 \hfill \\ s_{1}^{n} - s_{2}^{n},\quad k = 2 \hfill \\ L - s_{2}^{n},\quad k = 3 \hfill \\ \end{gathered} \right.\quad {\text{и}} \\ G_{{k,j}}^{n} \equiv \left\{ \begin{gathered} {{{{\xi }}}_{j}}(s_{1}^{{n + 1}} - s_{1}^{n}),\quad k = 1 \hfill \\ {{{{\xi }}}_{j}}(s_{1}^{{n + 1}} - s_{1}^{n}) + (1 - {{{{\xi }}}_{j}})(s_{2}^{{n + 1}} - s_{2}^{n}),\quad k = 2 \hfill \\ (1 - {{{{\xi }}}_{j}})(s_{2}^{{n + 1}} - s_{2}^{n}),\quad k = 3 \hfill \\ \end{gathered} \right.. \\ \end{gathered} $

Величины смещений границ $s_{1}^{{n + 1}} - s_{1}^{n}$ и $s_{2}^{{n + 1}} - s_{2}^{n}$ за время Δτ могут быть получены из условий Стефана для Γ12 и Γ23 (8), (9). Применительно к конечно-разностной схеме они выглядят следующим образом:

(13)
$\begin{gathered} s_{1}^{{n + 1}} - s_{1}^{n} = \\ = \frac{{\Delta {{\tau }}}}{{{{q}_{{F12}}}{{{{\rho }}}_{{12}}}}}\left[ {\frac{{{{{{\lambda }}}_{1}}}}{{s_{1}^{n}}}\frac{{{{\theta }}_{{1,N}}^{n} - {{\theta }}_{{1,N - 1}}^{n}}}{{\Delta {{{{\xi }}}_{{N - 1}}}}}} \right. - \left. {\frac{{{{{{\lambda }}}_{2}}}}{{s_{1}^{n} - s_{2}^{n}}}\frac{{{{\theta }}_{{2,1}}^{n} - {{\theta }}_{{2,0}}^{n}}}{{\Delta {{{{\xi }}}_{0}}}}} \right], \\ {{\theta }}_{{1,N}}^{n} = {{\theta }}_{{2,0}}^{n} = {{{{\theta }}}_{{F12}}}, \\ \end{gathered} $
(14)
$\begin{gathered} s_{2}^{{n + 1}} - s_{2}^{n} = \\ = \frac{{\Delta {{\tau }}}}{{{{q}_{{F23}}}{{{{\rho }}}_{{23}}}}}\left[ {\frac{{{{{{\lambda }}}_{2}}}}{{s_{1}^{n} - s_{2}^{n}}}\frac{{{{\theta }}_{{2,N}}^{n} - {{\theta }}_{{2,N - 1}}^{n}}}{{\Delta {{{{\xi }}}_{{N - 1}}}}} - \frac{{{{{{\lambda }}}_{3}}}}{{L - s_{2}^{n}}}\frac{{{{\theta }}_{{3,1}}^{n} - {{\theta }}_{{3,0}}^{n}}}{{\Delta {{{{\xi }}}_{0}}}}} \right], \\ {{\theta }}_{{2,N}}^{n} = {{\theta }}_{{3,0}}^{n} = {{{{\theta }}}_{{F23}}}, \\ \end{gathered} $
где θF12 и θF23 – температуры фазового перехода на границах Γ12 и Γ23, которые могут отличаться в силу тех или иных причин.

Таким образом, соотношения (10), (11), (12), (13) и (14) с учётом граничных и начальных условий (2) полностью описывают решение поставленной задачи. Переход от неравномерной сетки ξ к равномерной η выполняется исходя из задаваемой аналитической или табличной функции пересчёта координат, а преобразование от η к реальным координатам x вдоль оси абсцисс выполняется по соотношению (4).

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

(15)
$\begin{gathered} {{\xi }}\left( {{\eta }} \right) = \frac{{f({{\eta }}) - f(0)}}{{f(1) - f(0)}}, \\ f\left( {{\eta }} \right) \equiv {{\left[ {1 + \exp \left\{ { - {{\kappa }}\left( {{{\eta }} - 1{\text{/}}2} \right)} \right\}} \right]}^{{ - 1}}}, \\ \end{gathered} $
где κ – некоторый коэффициент, который определяет степень сгущения точек по краям. График зависимости ξ от η при различных κ, построенный по выражению (15), показан на рис. 3.

Рис. 3.

Расстояния между узлами неравномерной сетки при различных значениях κ и N = 5000: 1 – равномерная сетка; 2 – κ = 5; 3 – κ = 7; 4 – κ = 10; 5 – κ = 12; 6 – κ = 15; 7 – κ = 20.

Fig. 3. Distances between the nodes of an uneven grid for various values κ and N = 5000: 1 – uniform grid; 2 – κ = 5; 3 – κ = 7; 4 – κ = 10; 5 – κ = 12; 6 – κ = 15; 7 – κ = 20.

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

МОДЕЛИРОВАНИЕ ПРОЦЕССА ЗАМЕРЗАНИЯ ТРЕЩИНЫ

Сделаны приблизительные оценки времени замерзания трещины при различных параметрах. Выполнено моделирование, основанное на решении задачи Стефана с двумя фазовыми границами. Рассмотрен лишь процесс теплопереноса без учета движения воды, которое связано с её вытеснением из-за изменения объёма при замерзании. Считается, что температура фазового перехода не меняется. В реальности это не совсем так, поскольку она известным образом зависит от давления (Патерсон, 1984) в ходе замерзания воды будет постепенно нарастать на стенки трещины. Однако в рассматриваемом случае этого не происходит, поскольку вода вытесняется на поверхность избыточным давлением, приводя его к начальным значениям.

Расчёты выполнялись из предположения, что среды Ω1 и Ω3 одинаковы и представляют собой лед плотностью ρ1 = ρ3 = 910 кг/м3, теплопроводностью λ1 = λ3 = 2.22 Вт/(м К) и теплоёмкостью с1 = с3 = 2060 Дж/(кг К). Среда Ω2 считается водой ρ2 = 1000 кг/м3, λ2 = 0.569 Вт/(м К), с3 = 4212 Дж/(кг К). Скрытая теплота плавления льда составляет qF12 = qF23 = 3.32 × 105 Дж/кг (Патерсон, 1984). Поскольку изначально предполагается, что трещина заполнена талой водой, образованной на поверхности ледника, то её начальная температура соответствует температуре фазового перехода, то есть ${{{{\bar {\theta }}}}_{2}}$ = 0°C. Температуры сред Ω1 и Ω3 одинаковы, следовательно ${{{{\bar {\theta }}}}_{1}}$ = ${{{{\bar {\theta }}}}_{3}}$. Моделирование выполнялось при различных температурах ледника и начальной ширине трещины (то есть области Ω2).

На рис. 4 показан характер изменения распределения температуры в трёхслойной среде с течением времени. Для определённости выбрана температура ледника, которая может соответствовать району посадочных площадок станций Мирный или Новолазаревская в начале летнего периода, то есть ${{{{\bar {\theta }}}}_{1}}$ = ${{{{\bar {\theta }}}}_{3}}$ = –8°C (Попов и др., 2017). Начальная ширина трещины выбрана 10 см (см. рис. 1, б). Размеры областей Ω1 и Ω3 достаточно велики и составляют по 250 м каждая. Это позволяет во временных масштабах моделирования считать, что ледник безграничен в обе стороны от трещины, то есть существуют области в леднике, где он не успевает нагреться от тепла трещины.

Рис. 4.

Температура в трёхслойной среде с течением времени: белыми линиями показано положение стенок трещины.

Fig. 4. Temperature in a three-layer media with the time: position of the crevasse walls is depicted by white lines.

На рис. 4 продемонстрированы два одновременно протекающих процесса: растепление ледника от краев трещин вглубь и постепенное замерзание воды, начиная от стенок. На начальном этапе из-за большого перепада температур (8°C) скорость замерзания достаточно большая, но по мере нагревания стенок и уменьшения температурного градиента она естественным образом падает. Общее время замерзания трещины при заданных параметрах составляет 7.57 суток, что соответствует средней скорости 0.55 мм/час.

На рис. 5 показаны графики зависимости времени полного замерзания трещины от её начальной ширины и температуры окружающего льда в рамках рассмотренной модели. В тёплых горных ледниках, температура которых составляет первые градусы ниже нуля (Глазовский, Мачерет, 2014), трещины в течение летнего периода не замерзнут, но это произойдет зимой, когда температура приповерхностной части ледника понизится. Применительно к холодному антарктическому леднику видно, что в течение тёплого периода вода в трещине разумных размеров полностью замерзнет, что и наблюдается в районе станций Мирный и, по-видимому, на посадочной площадке станции Новолазаревская (см. рис. 1, а–б). Моделирование показало, что для прибрежной части холодного антарктического ледника со средней температурой –10°C и ниже (где в силу особенности его динамики и наблюдаются трещины), 5–10 см трещины замерзают менее чем за неделю. Более широкие трещины замерзают чуть дольше: замерзание 30 см требует от двух до трёх недель, в зависимости от температуры ледника.

Рис. 5.

Время замерзания воды в трещине при различной ширине и температуре ледяной толщи: температура льда: ‒0.5°C (1); –1.0°C (2); –5.0°C (3); –10.0°C (4); –15.0°C (5). Утолщенной линией показаны значения, полученные в рамках представленной модели, а пунктирной – результаты расчетов по соотношению (16).

Fig. 5. Freezing time of water in a crevasse with a different wide and temperatures of the ice: ice temperatures: –0.5°C (1); –1.0°C (2); –5.0°C (3); –10.0°C (4); –15.0°C (5). The values obtained by the presented model is shown by solid lines; the results of calculations according (16) is depicted by dotted lines.

Опубликовано множество статей, в которых приводится как аналитическое, так и численное решение задачи об определении скорости замерзания трещин, например, (Краслоу, Едгер, 1964; Тихонов, Самарский, 1977; Alley et al., 2005; Poinar et al., 2017). В работе (van der Veen, 2007) приводится соотношение для толщины намерзшего льда d за время t:

$d\left( t \right) = \frac{{2{{c}_{1}}\left( {{{{{\theta }}}_{{F12}}} - {{{{\theta }}}_{1}}} \right)}}{{{{q}_{{F12}}}\sqrt {{\pi }} }}\sqrt {{{a}_{1}}t} .$

Поскольку рассчитывается скорость намерзания от одной стенки, то для получения значения общей толщины эту величину нужно удвоить. Соотношения t, которое станет временем замерзания трещины заданной ширины D, D = 2d:

(16)
$t = \frac{{{\pi }}}{{16{{a}_{1}}}}{{\left[ {\frac{{{{q}_{{F12}}}d}}{{{{c}_{1}}\left( {{{{{\theta }}}_{{F12}}} - {{{{\theta }}}_{1}}} \right)}}} \right]}^{2}}.$

Расчёты, выполненные по соотношению (16), для сравнения с приведенной выше моделью, нанесены на рис. 5. Они показывают, что процесс замерзания трещины, описываемый моделью, происходит большей частью чуть быстрее, чем по уравнению (16). При этом с уменьшением температуры ледника уменьшается и степень расхождения результатов расчёта. В частности, 10-метровая трещина при разности температур льда и воды в 15°C замерзает за 2.153 суток по оценкам (16) и за 2.268 суток (то есть на 5.3% медленнее). Для разности температур в 1°C эти значения составляют 498.4 и 459.2 суток соответственно. То есть во втором случае процесс, описываемый уравнением (16), протекал на 7.9% медленнее. Изменение Δτ для моделирования (в пределах разумных значений числа Куранта) существенных изменений не дало. Указанные расчёты выполнены на неравномерной сетке при степени сгущения точек κ = 10 в выражении (15). Вычисления на равномерной сетке (при κ = 1) лишь незначительно изменяли расчётные значения.

Предположено, что несовпадение результатов не связано с ошибками моделирования, поскольку модель достаточно простая. Все дело в соотношении (16), про которое автор публикации (van der Veen, 2007) пишет, что это приближенная оценка.

ОБСУЖДЕНИЕ

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

Применительно к этой задаче, сотрудники Российской антарктической экспедиции (РЭА), обеспечивающие работу посадочной площадки станции Новолазаревской, рассказывали о том, что в тёплый период антарктического лета на взлётно-посадочной полосе (ВПП) протекают потоки талой воды и образуются небольшие трещины. Сам автор этого не видел, поскольку его работы в составе 67-й РАЭ (2021/22 г.) закончились до начала сезона таяния. В конце антарктической весны и начале лета открытых трещин на посадочной площадке обнаружено не было. Тем не менее, дайкообразные тела в пределах ВПП наблюдались (см. рис. 1, б). Это подтверждает рассказы сотрудников РАЭ. Таким образом, можно предположить, что трещины образуются главным образом в наиболее тёплую часть лета. Известно, что характер течения ледника, который зависит от вязкости льда, резко изменяется при температуре выше –10°C (Budd, 1969). Однако мощность ледника в районе ВПП составляет сотни метров, и сложно предположить, что в летний период температура всей толщи резко меняется, что приводит к существенному увеличению скоростей деформаций и возникновению трещин. Такая ситуация скорее возможна в районе сопки Ветров, где мощность ледника гораздо меньше.

Возможны два варианта: первый, когда трещины начинают возникать зимой. Это приводит к тому, что они, с одной стороны, постепенно расширяются, а с другой – забиваются снегом в результате выпадения твёрдых осадков и метелевого переноса. Затем, с приближением тёплого времени, начинается приповерхностное таяние и быстрое наполнение трещин водой с последующим полным или частичным замерзанием; второй вариант – это формирование трещины в тёплое время в период интенсивного таяния. В этом случае они, образуясь, заполняются талой водой и, судя по представленным графикам, практически сразу замерзают. Оценки времени замерзания трещины для обоих вариантов выполнены в рамках изложенной выше модели. Однако это всего лишь оценки, и для более корректного решения требуется гораздо более сложная двухмерная или трёхмерная модель, которая учитывает конфигурацию границ трещины.

ЗАКЛЮЧЕНИЕ

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

Благодарности. Автор выражает благодарность А.С. Борониной и двум анонимным рецензентам за конструктивную критику настоящей работы. Исследование выполнено при финансовой поддержке Российского научного фонда № 22-27-00266.

Acknowledgments. The author is grateful to A.S. Boronina and two anonymous reviewers for fruitful criticism. The study was financially supported by the Russian Science Foundation No. 22-27-00266.

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

  1. Глазовский А.Ф., Мачерет Ю.Я. Вода в ледниках. Методы и результаты геофизических и дистанционных исследований. М.: ГЕОС, 2014. 528 с.

  2. Казко Г.В., Саватюгин Л.М., Сократова И.Н. Моделирование циркуляции воды в антарктическом подледниковом озере Восток // Лёд и Снег. 2012. № 4. С. 86–91.https://doi.org/10.15356/2076-6734-2012-4-86-91

  3. Краслоу Г., Едгер Д. Теплопроводность твёрдых тел. М.: Наука, 1964. 488 с.

  4. Кольцова Э., Скичко А., Женса А. Численные методы решения уравнений математической физики и химии. М.: Юрайт, 2020. 220 с.

  5. Кузнецов Г.В., Шеремет М.А. Разностные методы решения задач теплопроводности: учебное пособие. Томск: Изд-во ТПУ, 2007. 172 с.

  6. Патерсон У.С.Б. Физика ледников. М.: Мир, 1984. 472 с.

  7. Попов С.В., Кашкевич М.П., Боронина А.С. Состояние взлетно-посадочной полосы станции Новолазаревская (Восточная Антарктида) и оценка безопасности её эксплуатации по данным исследований 2021 г. // Лёд и Снег. 2022. Т. 62. № 4. С. 621–636.

  8. Попов С.В., Поляков С.П., Пряхин С.С., Мартьянов В.Л., Лукин В.В. Строение верхней части ледника в районе планируемой взлетно-посадочной полосы станции Мирный, Восточная Антарктида (по материалам работ 2014/15 года) // Криосфера Земли. 2017. Т. XXI. № 1. С. 73–84.

  9. Рыбак О.О., Рыбак Е.А. Алгоритм решения системы уравнений течения льда в трёхмерной математической модели // Известия высших учебных заведений. Северо-Кавказский регион. Естественные науки. 2010. № 6. С. 117–121.

  10. Самарский А.А. Теория разносных схем. М.: Наука, 1977. 656 с.

  11. Смирнов В.И. Курс высшей математики. Т. 2. М.: Наука, 1974. 656 с.

  12. Тихонов А.Н., Самарский А.А. Уравнения математической физики. М.: Наука, 1977. 736 с.

  13. Alley R.B., Dupont T.K., Parizek B.R., Anandakrishnan S. Access of surface meltwater to beds of sub-freezing glaciers: preliminary insights // Annals of Glaciology. 2005. V. 40. P. 8–14.

  14. Budd W.F. The dynamics of ice masses. ANARE Sci. Rep. Publ. 1969. V. 108, 212 p.

  15. Greve R. A continuum-mechanical formulation for shallow polythermal ice sheets // Philos. Transaction Royal Society. London, 1997. V. 355. № 1726. P. 921–974.

  16. Greve R., Blatter H. Dynamics of ice sheets and glaciers. Springer Science & Business Media, 2009. 300 p.

  17. Huybrechts P. The Antarctic ice sheet and environmental change: a three-dimensional modelling study // Ber. Polarforsch. 1992. V. 99. 244 p.

  18. Nye J.F. Water flow in glaciers: jökulhlaups, tunnels, and veins // Journ. of Glaciology. 1976. V. 17. № 76. P. 181–207.

  19. Pattyn F. A new three-dimensional higher-order thermomechanical ice sheet model: Basic sensitivity, ice stream development, and ice flow across subglacial lakes // Journ. of Geophys. Research. 2003. V. 108. № B8. 2382 p.

  20. Poinar K., Joughin I., Lilien D., Brucker L., Kehrl L., Nowicki S. Drainage of Southeast Greenland Firn Aquifer Water through Crevasses to the Bed // Journ. of Front. Earth Sci. 2017. V. 5. 5 p. https://doi.org/10.3389/feart.2017.00005

  21. Thoma M., Grosfeld K., Mayer C. Modelling mixing and circulation in subglacial Lake Vostok, Antarctica // Ocean Dynamics. 2007. V. 57. № 6. P. 531–540.

  22. van der Veen C.J. Fracture propagation as means of rapidly transferring surface meltwater to the base of glaciers // Geophys. Research Letters. 2007. № 34. L01501. https://doi.org/10.1029/2006GL028385

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

Инструменты

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

Лёд и Снег

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