Задачи теплопроводности, включающие в себя фазовые превращения, встречаются в различных областях прикладной науки. Одним из важных и актуальных направлений, связанных с решением таких задач, является моделирование гляциальных процессов. Так, изучение и мониторинг состояния ледникового покрова в районах полярных станций Антарктиды является одним из ключевых факторов для поддержания безопасности инфраструктуры и авиационного сообщения, что имеет важнейшее значение для обеспечения деятельности Российской антарктической экспедиции (РАЭ). В представленной работе рассматривается процесс замерзания талой воды в ледниковых трещинах (т.н. «залечивание» трещин). Залеченные трещины наблюдались в районе взлетно-посадочной полосы антарктических станций Новолазаревская [1], Прогресс, полевой базы Оазис Бангера и других. Исследование эволюции трещин востребовано и для района сопки Ветров на станции Мирный, где возможно создание новой авиационной площадки.
Математическое описание процесса сводится к краевой задаче Стефана с подвижными фазовыми границами. Для численного решения такой задачи обычно рассматривают два типа методов [2]: с выделением границы фазового перехода (ловля фронта в узел пространственной сетки, метод выпрямления фронта) и методы сквозного счета (схема со сглаженными коэффициентами, энтальпийная формулировка). Для модели, представленной в данной работе, был выбран метод выпрямления фронта, так как он позволяет сформулировать задачу в статической области. Этот метод использовался для моделирования замерзания трещины в работе С.В. Попова [3] для одномерного случая, когда стенки трещины параллельны и пространство между ними заполнено водой. В настоящей работе рассматривается более общий двумерный случай, что позволяет смоделировать процесс замерзания в основании трещины, учесть распределение температуры льда по глубине и температурные колебания у поверхности льда, связанные с суточными или годичными изменениями температуры воздуха. Двухфазная модель строится для заполненной водой трещины под тонким поверхностным слоем талой воды на леднике. Такое явление наблюдается на плоских низменных участках в весенне-летний период.
Постановка задачи. Рассматривается двухфазная двумерная задача Стефана. Нижняя область Ω1 ∶{0≤x≤W, 0≤y< f(x,t)} занята льдом, верхняя Ω2 ∶{0≤x≤W,f(x,t)<y≤H} – водой. Граница фазового перехода описывается функцией f(x,t): 0≤x≤W,t>0, которая при каждом фиксированном t является однозначной, без самопересечений. В начальный момент времени значение f задано для всех x∈[0,W]. Предполагается, что среда является однородной и теплоемкость c, плотность ρ и коэффициент теплопроводности k не зависят от температуры. Также в рассматриваемой модели не учитывается конвективный перенос тепла. В таком случае, изменение функции температуры u(x, y, t) описывается уравнением теплопроводности следующего вида:
(1)
где m = 1 соответствует области со льдом Ω1, m = 2 – области с водой Ω2. Кроме того, заданы граничные и начальные условия:
(2)
где um,0 – начальные распределения температуры для льда и воды, g(t) – функция, задающая зависимость температуры на поверхности воды от времени, uice – некоторая фиксированная температура, up.t. – температура фазового перехода вода-лёд.
Постановку краевой задачи замыкает дополнительное условие на неразрывность теплового потока на границе раздела фаз (условие Стефана) [4], [5]:
(3)
Численная модель. Выпрямление фронта было проведено по координате y, с помощью преобразования:
(4)
Теперь переменная y′ изменяется в пределах от 0 до 2, причем твердой фазе соответствует область 0 ≤ y' < 1, жидкой фазе: 1 < y' ≤ 2, а y'=1 – координата границы фазового перехода. Также проведем масштабирование координаты x, чтобы ее значения изменялись от 0 до 1: x = x / W. Пусть u(x,y,t) = v(x',y',t), тогда после преобразования координат и обезразмеривания (t' = t/t0, v' = v/v0, где t0, v0 – характерные масштабы времени и температуры соответственно), краевая задача (1)-(3) может быть переписана следующим образом (штрихи у новых переменных опущены):
(5)
(6)
(7)
Разностная схема. Параболическое уравнение в частных производных вида (5), содержит смешанную частную производную, а также частную производную первого порядка, причем, коэффициенты при частных производных могут зависеть как от координат, так и от времени. Важно подобрать такую разностную схему, которая учитывала бы все эти факторы и при этом была устойчива, не ресурсоемка и имела достаточную точность. В работе реализована неявная схема переменных направлений [6], удовлетворяющая этим критериям. Она состоит из двух шагов, на каждом из них решается система с трехдиагональной матрицей. В случае равномерной сетки порядок точности схемы O(∆x2,∆y2, ∆t).
Неоднородная сетка. Для данной задачи целесообразно задать неоднородную сетку по координате y так, чтобы она сгущалась у границы фазового перехода. Это сделано отображением равномерной сетки на отрезке [0, 1] с помощью функции:
(8)
где s – фиксированный параметр, отвечающий за степень сгущения сетки. Полученные координаты симметрично относительно y = 1 отображаются на отрезок [1, 2].
На неравномерной сетке (рис. 1) реализованная схема имеет порядок точности O(∆x2,∆y, ∆t).
Программная реализация. Описанная выше модель была реализована на языке программирования Python с библиотекой NumPy. Следует отметить, что использование JIT-компилятора Numba [7] позволило на порядок повысить производительность расчетов.
Проверка корректности численной модели. Первый тест численной модели заключается в проверке согласованности результатов с полуэмпирической формулой для толщины намерзшего льда, приведенной в работе [8] в простейшем случае плоской границы фазового перехода, когда задача сводится к одномерной и однофазной. Параметры моделирования: начальная температура льда -10 °С, толщина льда в начальный момент времени – 10м, время моделирования – 35 суток. Численное решение со временем начинает отставать от вычисленного по формуле, но относительная погрешность остается практически неизменной и не превышает 2 %.
Второй тест представляет собой качественную проверку модели в двумерном двухфазном случае. Параметры моделирования: H = 1 м, W = 1 м, s = 8, число узлов – 200х200, шаг по времени – 0,5 часа, время моделирования – 200 дней. Лед в начальный момент времени полностью находился при температуре -5 °С, вода – при +5 °С. Температура воды на верхней границе и льда на нижней поддерживается постоянной. Как видно из временных срезов, для разной начальной формы фазовой границы, фронт асимптотически приходит к некоторому постоянному равновесному значению, рис. 2. Таким образом, было продемонстрировано, что при заданных стационарных условиях система приходит к равновесному состоянию, а конечное положение и форма границы не зависят от начальной конфигурации.
Моделирование замерзания воды в трещине. На рис. 3 представлен результат численного эксперимента, который проводился со следующими параметрами. Глубина щели – 5 м, ширина у поверхности ~ 10 см. Толщина поверхностного слоя воды – 5 см. Начальная температура воды равна 0 °C. Начальная температура льда линейно уменьшается от 0 °C на поверхности до -20 °C на глубине 10 м. На поверхности воды температура изменяется по синусоиде от 0 °C до 4 °C с периодом в 1 сутки. Начальное положение границы фазового перехода задано функцией Число узлов: 1000 по координате x и 200 по координате y. Параметр сгущения сетки: s = 8. Шаг по времени – 1 сек. Время моделирования – 82 часа.
За время моделирования глубина трещины уменьшилась на 4 метра. На рис. 3 также можно наблюдать области, где скорость замерзания изменялась из-за температурных колебаний воздуха.
Выводы. Представленная в настоящей работе простая двумерная модель позволяет получить адекватное описание процесса замерзания воды в ледниковой трещине. Дальнейшее развитие модели предполагает учет ряда дополнительных факторов, имеющих место в реальном физическом процессе, в частности, конвективный теплоперенос, солнечная активность и текучесть льда.