CC BY
32
Вестник МГТУ, том 16, № 4, 2013 г.
стр.641-643
УДК 622.413.4 : 551.34
Проверка кода для численного моделирования тепловых процессов в пористой среде с учётом фазового перехода "лёд - вода”
П.В. Амосов1,2
1 Горный институт КНЦ РАН
Л
Физико-энергетический факультет КФ ПетрГУ, кафедра теплофизики
Аннотация. В статье представлены предварительные результаты сравнения численных расчётов по моделированию тепловых процессов в пористой среде с учётом фазового перехода "лёд - вода". Расчёты выполнены с помощью двух верифицированных иностранных кодов и доработки автора к программе FFM, используемой в Горном институте КНЦ РАН. Выполненная проверка подтверждает возможность использования в будущих расчётах усовершенствованного кода.
Abstract. The paper presents some preliminary results of comparison of numerical calculations relating to the thermal processes modelling in a porous medium with phase transition "ice - water" taken into account. The computations have been carried out using two verified international codes and author's improvement to FFM software applied in the Mining Institute KSC RAS. The verification performed confirms possibility of the updated code using in future calculations.
Ключевые слова: численное моделирование, фазовый переход, "лёд - вода"
Key words: numerical modelling, phase transition, "ice - water"
1. Введение
Цель исследования - проверить доработку кода, разрабатываемого для численного моделирования тепловых процессов в пористой среде при наличии фазового перехода "лёд - вода". Код создаётся для использования при решении следующих проблем:
1) исследование теплового воздействия подземных атомных станций малой мощности (АСММ), размещаемых в труднодоступных районах Восточной Сибири на многолетнемёрзлые горные породы (ММГП);
2) оценка тепловой безопасности ММГП на потенциальном объекте окончательной изоляции ОЯТ Билибинской АЭС.
Первая проблема решается в рамках НИР "Научно-техническое обоснование и разработка концепции создания заглублённых и подземных АСММ модульного типа для энергоснабжения горнопромышленных предприятий и населённых пунктов в труднодоступных регионах России" (научные руководители: академик РАН Н.Н. Мельников, профессор В.П. Конухин).
Вторая проблема исследуется в рамках дипломной работы студентки физико-энергетического факультета КФ ПетрГУ Е.В. Резец (научный руководитель: к.т.н., доцент П.В. Амосов).
2. Алгоритм усовершенствования и результаты сравнения
В принципе, исследования теплового состояния ММГП можно было бы проводить на базе известных программных продуктов, например, PORFLOW (метод конечных разностей) или COMSOL (метод конечных элементов), которые позволяют решать подобные или достаточно близкие к ним задачи. Вместе с тем, иметь собственный "инструмент", который можно модернизировать под новые близкие проблемы (например, замёрзшие отвалы или хвостохранилища в условиях Крайнего Севера), несомненно лучше. Тем более, что собственный проверенный программный продукт можно использовать не только в научных исследованиях.
Модернизации подвергается одна из версий программы FFM, которая более 20 лет используется в Горном институте КНЦ РАН (разработчик: с.н.с. А.В. Наумов) при решении уравнения теплопроводности (метод конечных разностей, алгоритм Патанкара) (Мельников и др., 2001). Учёт фазового перехода выполнен по алгоритму, изложенному в русскоязычной литературе, например, в коллективной монографии сотрудников ИГД Севера им. Н.В. Черского СО РАН (Курилко и др., 2011).
Проверка программы осуществлялась на следующем модельном примере. Имеем замёрзшую (-20 °С) пористую (10 %) породу размером 1x1 м. Граничные условия следующие: на трёх границах (нижней и боковых) условия изоляции, на верхней - фиксированное значение температуры (100 °С).
641
Амосов П.В. Проверка кода для численного моделирования...
Сравнительному анализу подвергается центральное вертикальное распределение температуры на 24 часа процесса оттаивания.
Естественно, что первоначально описанная задача была решена с помощью указанных выше верифицированных кодов, и было достигнуто, на взгляд автора, хорошее совпадение результатов (см. таблицу). При этом параметры моделей, описывающих фазовый переход, приняты по умолчанию в соответствии с рекомендациями разработчиков.
Коротко опишем подход, используемый д.т.н. Курилко А.С. с коллегами (Курилко и др., 2011). В двухмерной постановке процесс распространения тепла в массиве горных пород с учётом фазовых переходов описывается следующей системой уравнений и условий:
[C(T) + LWpS(T- T )](dT/dt) = д[Л(Т)(дТ/дх)]/дх + д[Л(Т)(дТ/ду)]/ду,
\ClPl, Т < Т,
С(Т) = \
Т > Т,
\ Л l, Т < Т*,
3 II
\ Л2, Т > Т,
0 < x < Hx, 0 < У < Ну,
t > 0,
где Т - температура породы, К; Т* - температура фазовых переходов влаги в породе, К; t и x,y -временная (с) и пространственные координаты (м); L - скрытая теплота плавления (замерзания) льда (воды), Дж/кг; W - влажность породы, доли единицы; р - плотность воды, кг/м3; сь р\. Л\ (с2, р2. /-г) -удельная теплоёмкость (Дж/(кт-К)), плотность (кг/м3) и коэффициент теплопроводности (Вт/(м К)), соответственно для мёрзлых (талых) пород; 8(Т-Т) - S-функция Дирака; Hx и Hy - границы области моделирования по осям координат, м.
Из множества существующих вариантов аппроксимации эффективной теплоёмкости широкое применение получили методы, в которых влияние S-функции Дирака распространяется только на две смежные точки пространственной сетки (в нашем случае попарно вдоль продольной и поперечной осям) и вариант, когда эффективная теплоёмкость непрерывна в точках Т -АТ и Т + АТ.
В соответствии с методом сглаживания S-функция Дирака приближённо заменяется S-образной функцией 8(Т-Т ,АТ), которая отлична от нуля на интервале (Т -АТ, Т +АТ) и удовлетворяет условию нормировки
Т*+АТ
\8(Т - Т*,АТ) = 1.
Т* -АТ
Записывается эффективная теплоёмкость [С(Т)] = С(Т) + LWpS(Т- Т ,АТ), которая далее интегрируется в пределах (Т-АТ, Т +АТ). В результате получаем следующее условие постоянства энтальпии на интервале сглаживания
Т*+АТ
\[С(Т)]ёТ = LWp + (ср + СР2)АТ.
Т*-АТ
Предлагается эффективную теплоёмкость записывать в следующем виде
\С\Р\,
[С{Т)] = \(с1р1 + с2р2)/2 + (с1р1- с2р2)(Т- Т*)/2/АТ + LWp(l -\Т- Т*\/АТ)/ АТ,
\ c2p2,
Т < Т*- АТ; \Т- Т*\ < АТ; Т > Т* + АТ.
Разрывность коэффициента теплопроводности устраняется посредством соединения точек
(Л1,Т*- АТ) и (Л2,Т* + АТ) прямой линией
\ЛХ, Т < Т*- АТ;
[Л{Т)] = М + Л2)/2 + (Л2-Л1)(Т- Т*)/2/АТ, \Т- Т*\ < АТ;
\Л2, Т > Т* + АТ.
Таким образом, начальное уравнение теплопроводности при численной реализации заменяется уравнением
[С(Т)](дТ/д() = д[[Л(Т)](дТ/дх)]/дх + д[[Л(Т)](дТ/ду)]/ду, 0 < x < Hx, 0 < y < Hy, t > 0.
В начальный момент времени задаётся распределение температуры:
642
Вестник МГТУ, том 16, № 4, 2013 г.
стр.641-643
T(x,y,0) = T0, 0 < х < Hx, О < y < Hy.
В задаче, выбранной для проверки, использованы следующие граничные условия:
на верхней границе - T(x,Hy,0) = Tup, 0 < х < Hx (условие Дирихле);
на нижней и боковых границах - X(dTldy) \y=0 = 0, 0 < x < Hx и
X(dT/dx)\x=0x=Hx = 0, 0 <y <Hy (условие Неймана, в частности, нулевой поток тепла).
Результаты численных экспериментов с использованием всех программных продуктов сведены в
таблицу.
По мнению автора, полученные результаты позволяют с определённым оптимизмом смотреть на возможное использование в будущих расчётах собственного кода. Действительно, согласие неплохое, хотя некоторое расхождение в цифрах наблюдается. Например, собственный код выдаёт заниженные положительные температуры (относительная ошибка по модулю 7,8 %), но в области фазового перехода собственные расчётные значения попадают в середину интервала верифицированных кодов.
В качестве примера потенциального применения программы может быть задача по определению скорости протаивания замёрзшей горной породы при наступлении оттепели. В плоской геометрии на модели размером 30x20 м прослежено продвижение "фронта" нулевой температуры вглубь массива при задании на поверхности условия 3-го рода. Выполненная оценка позволила получить значение скорости протаивания на уровне 3-4 см/сутки за время порядка 60 суток. Результат достаточно физичен и отвечает набору исходных параметров.
Таблица. Пространственное распределение температуры на 24 часа процесса оттаивания, °С
Координата, м Код COMSOL Код PORFLOW Собственный код
0,0 -17,2 -17,1 -17,6
0,1 -16,7 -16,6 -17,1
0,2 -14,7 -15,0 -15,6
0,3 -12,3 -12,1 -12,7
0,4 -7,9 -7,7 -8,1
0,5 -1,8 -0,9 -1,3
0,6 10,7 10,9 9,7
0,7 28,2 28,2 26,0
0,8 49,6 49,6 47,2
0,9 74,1 74,0 72,5
1,0 100,0 100,0 100,0
3. Заключение
Достигнутая точность расчётных значений температуры в сравнении с результатами, полученными с использованием верифицированных программ POFRLOW и COMSOL, позволяет автору рекомендовать созданную программу для численного моделирования процессов переноса тепла при наличии фазового перехода "лёд - вода".
Литература
Курилко А.С., Ермаков С.А., Хохолов Ю.А., Каймонов М.В., Бураков А.М. Моделирование тепловых процессов в горном массиве при открытой разработке россыпей криолитозоны. Новосибирск, Академическое изд-во "Гео", 139 с., 2011.
Мельников Н.Н., Наумов В.А., Конухин В.П., Амосов П.В., Гусак С.А., Наумов А.В.
Радиогеоэкологические аспекты безопасности подземного захоронения радиоактивных отходов и отработавшего ядерного топлива на Европейском Севере России. Апатиты, КНЦ РАН, 194 с., 2001.
643
