CC BY
28
□
□
УДК 519.63
B.C. Афанасьев
ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧИ ЗАМОРАЖИВАНИЯ РАССОЛА, МОДЕЛЬ С ПРОТЯЖЕННОЙ ОБЛАСТЬЮ ФАЗОВОГО ПЕРЕХОДА
В данной работе рассматривается модель замораживания рассола с протяженной областью фазового перехода, называемой двухфазной зоной. В данной модели существуют две границы разделения фаз: «лед-двухфазная зона» и «двухфазная зона-рассол.». Для численной реализации используется метод ловли фронта «лед-двухфазная зона» в узел пространственно-временной сетки, вторая граница ищется из условия термодинамического равновесия в двухфазной зоне.
Процессы фазовых переходов и тепломассопереноса в значительной степени определяют свойства мерзлых грунтов и пород. Знание закономерности промерзания и протаивания грунтов, умение прогнозировать развитие этих процессов применительно к различным техническим аспектам освоения северных территорий позволяют вырабатывать оптимальные решения по строительству наземных и подземных сооружений в зоне многолетней мерзлоты, применению искусственного замораживания грунтов в строительстве и другим инженерным проблемам.
Можно считать, что грунт представляет собой пористую среду, заполненную водой с растворенными в ней солями. В зависимости от физической постановки задачи применяются два различных способа математического моделирования замораживания растворов: классическая фронтовая задача Стефана с четкой границей раздела жидкой и твердой фаз [1] (и ее модификации) и модель с фазовым переходом в протяженной области.
Как показали экспериментальные и теоретические исследования, замораживание водонасыщенных грунтов сводится не только к распространению фронта фазового перехода, но и сопровождается перераспределением влаги и растворенных в ней примесей, при этом поскольку плотность льда меньше плотности воды, имеет место тепло- и массоперенос в жидкой фазе, деформирование и пучение грунта. Область эффективного применения фронтовой термодиффузионной модели замораживания водного раствора соли оказывается незначительной. То есть фронтовая модель в широком диапазоне изменения параметров оказывается противоречивой [2, 5].
Для адекватного описания процесса замораживания в тех случаях, когда фронтовая модель становится противоречивой, В.М. Ентов, А.М. Максимов и Г.Г. Цыпкин [3, 4, 6] построили корректную математическую модель, в которой фазовый переход вода-лед происходит в протяженной области, названной ими двухфазной зоной.
В двухфазной зоне лед и вода сосуществуют в состоянии термодинамического равновесия. Данная модель основана на фундаментальных законах механики многофазных сред.
В работе представлены вычислительный алгоритм, реализующий модель замораживания с протяженной областью фазового перехода, и результаты расчетов.
Постановка задачи
Рассмотрим модель с двухфазной зоной кристаллизации толщи морской воды.
Закон сохранения энергии в твердой фазе:
Закон сохранения энергии для двухфазной и талой зон:
Ср1г“К*!)-'х^ь}■1 >а (2)
Уравнение диффузии соли в двухфазной зоне:
где Ср1, С11 и Ср2, С12 - объемная теплоемкость и теплопроводность льда и воды соответственно, Ср = Ср1(1 -п) + Ср2п, I = 1(1 - п) + 12п; Б - коэффициент диффузии соли в воде.
Условие термодинамического равновесия для двухфазной зоны:
Т = -ус, х е [<^,^], t > 0. (4)
Значения влажности и концентрации соли в талой зоне:
с(х,t)= С0 , х е[^,Ь\, t > 0 , (5)
у(х, 0 = 1, х е[^,ь\, t > 0 . (6)
Условия на границе раздела твердой и двухфазной зон: —
ііш Т = ііш Т = Т*, ііш с = 0 , ііш с = с*,
х^%-0 х^% + 0 х^%-0 х^| + 0
Ііш V = 0 , Ііш V = V*, ґ > 0.
х^|-0 х^| + 0
(7)
Баланс энергии на границе раздела твердой и двухфаз ной зон:
дТ
А-
дх
Баланс массы соли на границе раздела твердой и двухфазной зон:
Б — + с* = 0 , х = Е (t), I > 0 . (9)
дх Л
Граничные условия для температуры:
Т = Т1, х = 0, t > 0 (10)
Т = Т0, х = Ь . (11)
Условие нахождения правой границы двухфазной зоны:
Т <-ус0, х е(0,^). (12)
Сглаживание распределения влажности в двухфазной зоне целесообразно проводить с помощью введения искусственной вязкости:
д у , х е(£,^), t > 0 (13)
ду
аГ
с граничными условиями: ду
дх
= 0, х = £, ґ> 0,
(14)
у = 1 , х = ^, ґ> 0.
Начальные значения температуры, концентрации соли
и влажности:
Т(х,0) = Т0, х е[0, Ь ], с(х,0) = с0, х є [0, Ь ], у(х,0) = 1, х є [0,Ь].
(15)
(16) (17)
= {*, = х,._! + кг, 7 = 1, п; х0 = 0; кг = ак1_1, 7 = 1, п}
шаги которой образуют возрастающую геометрическую прогрессию со знаменателем а >1. При этом выберем параметры сетки п, И0 и а таким образом, что хп = Ь. Очередной временной шаг определяется в ходе решения поставленной задачи. Вторая граница - между двухфазной и талой зонами находится из условия (12); как только оно перестает выполняться, начинается талая зона - 7 = т, п.
Пусть найдены приближенные решения исходной задачи до момента времени / = Построим вычислитель-
ный алгоритм, осуществляющий переход на следующий временной слой.
Уравнению (1) поставим в соответствие его чисто неявный дискретный аналог
Ъ.Ср, ^ = Тм ~ Т‘ - Т ~ Т‘-1 , 7 = 17^1. (18)
' V ьм к 57
Уравнение (2), с помощью (3) и (4) приведем к виду
Ср--Ьр]— = + —
Т ) дґ дх дх
х є (|, Ь), ґ > 0 ,
который будем использовать для вычисления значений температуры в двухфазной зоне.
Поставим ему в соответствие дискретный аналог
(
\
СРг - уЬР
Т - Т
(
ТТ Т-Т
\
У(І2 -І1 --
ОЬрЛ тм - Т
т ),+1 ъм
2
(19)
БЬрЛ Т - Т_1
т 1 к
Краевая задача (1)-(17) служит двухфазной моделью замораживания толщи морской воды.
Алгоритм решения поставленной задачи
Рассмотрим вычислительный алгоритм, пригодный для численной реализации краевой задачи (1)-( 17). Он основан на идее ловли границы раздела мерзлой и двухфазной зон в узел пространственно-временной сетки. На сегменте [0, Ь] введем квазиравномерную сетку
7 = ] +1, т -1.
Численную реализацию уравнения (2) для нахождения значения влажности в двухфазной зоне проведем в три этапа, используя принцип аддитивности [7].
На первом этапе определяем вспомогательное распределение температуры как решение линейного уравнения теплопроводности:
Ср
и - Т
(ї,гі) ■
На втором этапе для распределения температуры на верхнем временном слое выписываем схему Роте для нахождения промежуточных значений температуры
T-U д ( dT \ vs-V
CP-------- = I-LP------, x Є
г ox ^ ox) т
Используя (3) и (4), получаем
vs.. = -
DLp + D(Cp2 - CA)(- - Ut) + (k,-k1)ycl
(20)
V- 1 — V- v. — v. 1
i+1 і і і—1
л
V hi+1
(21)
Х2т
T
Л
ST
dx
- DLpv, ^T = 0, t > 0. dx
Полученному соотношению поставим в соответствие его дискретный аналог
T
V
t/ - t/ 1
h
-я-
T - T
j______ІІ
\
h
'j+1
-LpDv—J+1 -J = 0.(24)
jh
nj+1
ci = c , і = m, n,
(2З)
Из (9) получаем
7 (vj-jhj+1 - Vj- (hj + j) - ^-1-/-А ) =
T - U vsT + yvc vs -v , ч
CP----= (^2 -^1)----------------------------------------^-LP-, x є УЇЛ).
t Dr t
Из дискретного аналога полученного уравнения определяем распределение вспомогательной влажности в двухфазной зоне:
БЬрї, - бср^т - и,)+(^ - ^)те
= (v, +^) D-jf~~Tl hj+1
(29)
7 = У,т .
После этого решаем систему уравнений (13)-(14) -искусственной вязкости, которое будем использовать для нахождения неколеблющихся значений влажности.
7 = у +1, т -1.
Уравнению (4) поставим в соответствие его дискретный аналог: ___________
Т = ~усi, г = ] +1, т -1. (22)
Значения температуры в талой зоне находим из следующего уравнения, полученного из (2), используя (6)
ГТ1 ГТ7 ГТ1 ГТ1 ГТ1 ГТ1
к1Ср2 Т^Т- = 7+1 ~ 7 —, 7 = т +1, п -1.(23)
Таким образом, для перехода с (У - 1)-го нау'-й временной слой, необходимо решить нелинейную систему (18)-(29). Чтобы определить распределения температуры, влажности и концентрации соли, необходимо знать значение временного шага и границу раздела талой и двухфазной зон. И наоборот, вычисление временного шага и границы возможно лишь по известному распределению температуры. Поэтому для реализации полученной нелинейной системы уравнений приходим к необходимости организации итерационного процесса.
Шаг 1. Положим к = 0 (счетчик итераций) и зададим начальное приближение временного шага /0 (например, т0 = т). Естественно, придется отдельно задавать началь-ное приближение первого временного шага, например, по формуле
г0 = 0.5Ьрй02/[А1 (т0 -Tg)]
Шаг 2. Произведем простейшую линеаризацию системы уравнений (18)-(22), (25)-(28):
rpk+1 ___ 'Г'1
-0 ~ 1 ,
rpk+1 гр
h1CP1 T—T-
Л1Т
rpk+1 у і 1 i+1 _ 1 1
k+1
h
i+1
—k+1 —і
h
k+1
Из уравнения (8), используя (9), исключим слагаемое, содержащее скорость движения границы раздела мерзлой и двухфазной зон:
і = 1, j -1,
(Cpk ~^LP)
— k+1 _____—
f —k+1 —k+1 —k+1 —k+1 ^
1 i+1 ~ 1 і _ -i ~—-1
V hi+1 hi
~k~
DLp)k -++1 - Tk+1
~ 11 “X 2
Естественным образом условия (З), (б), (10) и (11) аппроксимируются уравнениями
-л-
DLp)k -k+1 - T^1
V; = 1, і = m, n .
-0 = T1,
— ___— о
і = j +1, m -1,
rpk+1 nr rpk+1 rpk+1 rpk+1 rpk+1
HCp2 T ~T - T+1 _ T T ~1i-1
X2t
h
h
- vsi
h
і = т +1, п -1,
^ = РЬру, - РС Рісґ 1 - и,+1)+(Я, -^Ж-С, ' РЬр + Р(СР2 - СРі)(Тк+1 - и,+1) + (^ -
ик+1 - у/+1
= *к
, = 7, т,
V:1 - V" ,+1 ,
.к+1
ї V
к
к
7 = у +1, т -1,
Значения вспомогательного распределения температуры и влажности вычисляются на каждой итерации.
Тк+1 = -аск+1, 7 = У +1, т -1,
V
к+1
= 1, і = т, п ,
к+1 0
с, = с
, = т, п ,
Т к+1 т 0
п
Для решения полученной системы линейных трехточечных разностных уравнений относительно Ті (/ = 0, п) целесообразно использовать модифицированный вариант алгоритма встречной прогонки. Его суть заключается в том, что решение рассматриваемой системы ищется в виде
Т,к+1 =(1 -а,)тк+1 +Рі , і =
Т- +1 =(1 -а,)тД+1 + Д , 7 = у + 1, п
Подставляя выписанные соотношения в соответствующие уравнения рассматриваемой системы, получаем рекуррентные формулы для прогоночных коэффициентов
, Д, 7 = 0, у -1, 7 = у +1, п .
Для нахождения значений водонасыщенности у^ = у, т) используется обычная прогонка, то есть решение ищется в виде
у^ =0^ +Д , 7 = т^Ту . ________________
Прогоночные коэффициенты , Д , 7 = у +1, т вычисляются из соответствующего уравнения системы.
Шаг 3. Линеаризуем оставшиеся два уравнения:
(
Тк
V
грк+1 п-гк+1
Т ~ТТ-1
-Я-
грк+1 грк+1 \
ТТ+1 ~ ТТ
(30)
+ ЬрРу
Тк+1 тк+1
к Т+1
= 0,
7к“Т*X. -Тт(кТ + АТ..)-*ТТт)=
=(V;
7+1
)Р
т;+1 грк+1
ТТ+1 _ ТТ
(31)
к
7+1
-гк+1
И используем ИХ ДЛЯ определения Гк+1, Ту
Исключая из уравнения (30) Т^1 и ТД!, получаем
квадратное уравнение относительно температуры фазового перехода Т.. Физический смысл имеет меньший из корней полученного уравнения. Далее из уравнения (31) определяем очередное приближение временного шага.
Шаг 4. Рассмотрим ошибку итерационного процесса
£ = К+1/тк -!|.
Если 5 > е, где е - заданное малое число, то увеличиваем к на единицу и возвращаемся к выполнению шага 2 итерационного процесса. В противном случае считаем, что переход на следующий временной слой осуществлен. В качестве искомых значений временного шага, распределений температуры и концентрации примеси и водонасыщенности берем соответствующие значения, полученные на последней итерации.
Результаты численной реализации
В результате численной реализации получены следующие графики:
Рис. 1. Г1 = 00 С, Т1 = -100 С
7+1
к
Рис. 4. Т° = 3°С, Т1 = -2°С
Рис. 2. Т° = 2°С, Т1 = -10°С
Здесь и далее обозначения кривых: 1 - распределение температуры, 2 - концентрация соли, 3 - распределение влажности. Т0 - начальная температура, Т1 - температура на правой границе.
Рис. 3. Т° = 0°С, Т1 = -2°С
12 24 36
Рис. 5. Т° = 4°С, Т1 = -20С
АО 60 120 160 200 Х.см.
Рис. 6. Г1 = 4°С, Т1 = -З0С
т,°с ю*си
АО 80 120 160 200 X, см
Рис. 7. Т° = 2°С, Т1 = -30С т„°с з*си
_ Рис. 8. Т° = 3°С, Т1 = -13°С
т, с з*си
т.°с з*си
50 ЮО 150 200 250 Х,см.
Рис. 10. Т° = 0°С, Г1 = -15°С
Графики на рис. 1-рис. 1° показывают, что математическая модель адекватно описывает изучаемый процесс. Видно, что толщина полностью промерзшего слоя растет крайне медленно и что он намного тоньше, чем двухфазная зона. Это можно объяснить тем, что диффузия соли в воде происходит довольно медленно, поэтому оставшаяся соль не дает воде полностью замерзнуть.
Литература:
1. Васильев В.И. О фронтовой модели замораживания толщи раствора // Математические заметки ЯГУ. 1994. Т. 1. № 1.
C. 136-144.
2. Васильев В.И. Численная реализация моделей замораживания водонасыщенного грунта // Математическое моделирование. 1995. T. 7. № 8. C.91-1°4.
3. Васильев В.И., Максимов А.М., Петров Е.Е, Цыпкин Г.Г. Тепломассоперенос в промерзающих и протаивающих грунтах. М.: Наука. Физико-математическая литература, 1996. 224 с.
4. Ентов В.М., Максимов AM., Цыпкин Г.Г. Образование двухфазной зоны при промерзании пористой среды. М.: Институт прикладной механики. АН СССР. Препринт №269. 1986. 56 с.
5. Забелина М.П., Фрязинов И.В. Сеточный метод решения задачи Стефанадля бинарной системы // Дифф. уравнения. 1987. Т. 23. № 7. С. 1188-1197.
6. Максимов AM., ЦъткинГ.Г. Образование двухфазной зоны при взаимодействии талых и мерзлых пород с раствором соли. М.: Институт прикладной механики. АН СССР. Препринт № 3°5. 1986. 6° с.
7. Самарский А.А. Теория разностных схем. М.: Наука 1983. 616 с.
VS. Afanasiev
NUMERICAL SOLUTION OF THE SALT-WATER FREEZING TASK, A MODEL WITH EXTENDED REGION OF PHASE TRANSITION
In the article a model of salt water freezing with extended region of phase transition called two-phase region has been studied. The model contains 2 borders of phase division: “ice-two-phase zone” and “two-phase zone-salt water”. For numerical realization a method of battle-front catching “ice-two-phase zone” into the space-time frame juncture is used. The second border is calculated from the thermo-dynamic balance condition in the two-phase zone.
