Валидация численной модели плавления октадекана с помощью лабораторного эксперимента Текст научной статьи по специальности «Физика»

УДК 536.242

DOI: 10.18698/2308-6033-2023-8-2296

Валидация численной модели плавления октадекана с помощью лабораторного эксперимента

© Р.А. Воропаев, И.С. Мацак ПАО «РКК «Энергия», Королёв, 141070, Россия

Представлены результаты исследования процесса плавления теплоаккумулирую-щего вещества — октадекана, используемого в космическом тепловом аккумуляторе. Задача Стефана, формулируемая для границы фазового перехода, имеет аналитическое решение лишь для одномерного случая, а для реальной геометрии конструкции может быть решена только с применением численных методов, в частности, посредством моделирования в универсальной программной системе анализа конечных элементов ANSYS Fluent. Для валидации разработанных численных моделей проводили модельные эксперименты по плавлению октадекана в цилиндрической стеклянной колбе с теплоизолированными боковыми стенками, в нижнем торце которой расположен нагреватель. Были исследованы двумерная и трехмерная численные модели при мощностях нагрева 7,1 Вт и 34,2 Вт. Проведено сравнение динамики доли расплавленного вещества в двумерной и трехмерной моделях с экспериментальными данными. Учтены тепловые потери на нагрев конструкции и конвективный теплообмен с внешней средой. Представлены результаты численных и лабораторных экспериментов, которые можно использовать для проверки модели нагрева теплоаккумулирующего вещества с постоянной мощностью. Полученные результаты показывают, что при абсолютном изменении температуры теплоаккумулирующего вещества на небольшую величину (до 10 °С) приближение Буссинеска дает результаты, хорошо верифицируемые с экспериментальными данными.

Ключевые слова: тепловой аккумулятор, плавление, моделирование, приближение Буссинеска, ANSYS Fluent

Обозначения

g ускорение свободного падения, м/с2 Sh источник энтальпии в уравнении энергии, Вт/м3

ß коэффициент теплового расширения, К-1 Su источник импульса в уравнении движения, Па/м

T температура, °С u вектор скорости, м/с

W ширина, м Ra число Релея

V кинематическая вязкость, м/с2 Ste число Стефана

X температуропроводность, м/с2 Fo число Фурье

c теплоемкость, Дж/(кг . К) Pr число Прандтля

cp изобарная удельная теплоёмкость, Дж/(кг . К) Nu число Нуссельта

L скрытая теплота, Дж/кг Gr число Грасгофа

t время, с Q потери теплоты, Дж

P плотность, кг/м3 F площадь поверхности, м2

k коэффициент теплопроводности, Вт/(м . К) 5 толщина корпуса, м

динамическая вязкость, Па . с m масса, кг

A соотношение сторон V объем, м3

H высота, м e доля расплавленного ТАВ

e единичный вектор в направлении оси OZ У коэффициент отклонения флуктуации

F член плавучести, Па/м h удельная энтальпия, Дж/кг

p давление, Па

Введение. В космической технике существует проблема обеспечения теплового режима аппаратуры, тепловыделение которой превышает несколько сотен ватт, так как удельная производительность холодильников-излучателей находится в диапазоне от 100 до 200 Вт/м2 при характерной температуре 35 °С. Увеличивать размеры холодильника-излучателя для непрерывно-периодического режима работы аппаратуры нецелесообразно [1]. Для решения задач такого класса можно применять тепловые аккумуляторы, осуществляющие отбор и хранение теплоты, которая впоследствии утилизируется излучением с поверхности значительно меньшего размера, например площадью 1 м , при мощности излучения около 100 Вт. Однако при проектировании систем хранения скрытой тепловой энергии возникает проблема, заключающаяся в том, что аналитически можно оценить лишь энергоемкость. В то же время мощность и перепад температур в начале и конце цикла плавления зависят от геометрии, свойств материалов, ориентации и других параметров, поэтому они могут быть достаточно точно определены только экспериментальными методами, а спрогнозированы — численным моделированием. В связи с этим авторами было предложено провести валидацию численной модели плавления теплоаккумулирующего вещества (ТАВ) октадека-на, используемого в космическом тепловом аккумуляторе, с помощью лабораторного эксперимента.

В качестве рабочего тела тепловых аккумуляторов можно применять большое количество ТАВ, температура плавления которых соответствует выбранной задаче [2]. Теплоаккумулирующие вещества используются для систем тепловых насосов, солнечной техники, систем терморегулирования космических аппаратов, а также в строительстве [3], что делает актуальной задачу исследования процессов теплопередачи во время плавления. Исследования проводят с применением численных [4, 5] и экспериментальных методов [6]. Особое внимание уделяется валидации численных моделей. Но следует отметить, что экспериментально-расчетные исследования, соответствующие трехмерным моделям, пока еще немногочисленны, хотя и представляют интерес для проектирования и оптимизации работы тепловых аккумуляторов.

Ранее проведенные исследования фазовых переходов у конструктивно различающихся вариантов тепловых аккумуляторов показывают, что необходимо учитывать особенности происходящих тепловых источников. Например, в работе [7] представлено экспериментальное исследование влияния угла наклона прямоугольной камеры на естественную конвекцию при плавлении материала с фазовым переходом. Было показано, что теплопередача для случая горизонтального расположения нагревателя усиливается более чем в 2 раза по сравнению

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

Для определения характеристик теплоотвода при фазовом переходе предполагается решение задачи Стефана для движущейся границы ТАВ. Сложность решения таких задач обусловлена тем, что закон движения раздела фаз определяется конвекцией в жидком ТАВ, граничными условиями и скоростью поглощения теплоты при фазовом переходе. При этом положение границы фазового перехода заранее неизвестно и является частью решения. Аналитическое решение данной задачи существует только для одномерного случая [9]. Для остальных вариантов геометрии используют численные методы, поэтому тепловой расчет проводится обычно в среде ANSYS Fluent или в других системах моделирования.

В работах [6, 10] приведены результаты сравнения двух подходов к моделированию процесса плавления для случая нагрева ТАВ от боковой стенки при постоянной температуре. Первый подход, именуемый методом объема жидкости (Volume of fluid method — VOF), основан на модели, которая учитывает зависимость теплофизических свойств от температуры, в частности, заполнение воздушной прослойки расплавленным веществом при его расширении. Во втором подходе используется приближение Буссинеска для плотности, в то время как для остальных теплофизических свойств принимаются постоянные значения. Результаты валидации метода VOF показывают незначительные отклонения при расчете (до 4 %) доли жидкой фазы и суммарного теплового потока, однако для скоростей и температур различия более заметны (до 20 %). По мнению авторов [10], метод VOF более предпочтителен для расчета систем с высоким градиентом температур, однако использование приближения Буссинеска дает достаточно точные результаты на этапе проектирования. Поэтому в качестве основного метода выбрано приближение Буссинеска как более быстрый способ расчетов тепловых аккумуляторов в цилиндрических корпусах.

Цель работы — анализ адекватности результатов при использовании различных подходов к моделированию процесса плавления путем сравнения с натурным экспериментом, а также оценка точности численных расчетов с учетом естественной конвекции.

Экспериментальная установка. Для валидации численных моделей был изготовлен лабораторно-исследовательский стенд (рис. 1) в виде цилиндрической стеклянной колбы с теплоизолированными стенками, заполненной ТАВ. В нее помещен нагреватель постоянной мощности, состоящий из нихромовой нити и алюминиевой пластины, обеспечивающей равномерность теплового потока. С помощью лабораторного автотрансформатора (ЛАТР) задается мощность нагрева. Для измерения температуры в объеме ТАВ размещены цифровые датчики температуры (ДТ) диаметром 4 мм модели DS18B20+ с пластмассовым покрытием, расположенные на одной вертикали (рис. 2).

Датчики опрашивались с помощью программного продукта Arduino Uno. Для обработки данных использовалось программное обеспечение MATLAB. Доли расплавленного ТАВ и скорости в жидкой фазе определяли с помощью видеокамеры. Благодаря равномерному тепловому потоку от пластины ТАВ плавилось равномерно по вертикали во всем объеме, благодаря чему коэффициент расплавления или границу фазового перехода можно было определить по нанесенной на колбе шкале. Для вычисления коэффициента теплоотдачи с поверхности теплоизоляции и ТАВ использовались значения температуры, регистрируемые тепловизором.

Рис. 1. Схема лабораторно-исследовательского стенда:

1 — тепловизор; 2 — стеклянная колба; 3 — Arduino Uno; 4 — ПК; 5 — камера;

6 — алюминиевая пластина; 7 — нагреватель; 8 — ЛАТР; 9 — теплоизоляция

Габаритные размеры установки указаны на рис. 2. ТАВ размещено в вертикально ориентированной колбе, имеющей внутренний диаметр W, равный 63,5 мм, и высоту 403,8 мм. Корпус колбы, сделанный из стекла, имеет толщину стенки 2,5 мм. Круглый нагреватель из ни-хромовой нити толщиной 2 мм и диаметром 63,5 мм, закрепленный на алюминиевой пластине толщиной 7,8 мм, обеспечивал равномерный

нагрев ТАВ. Питание нагревателя осуществлялось от лабораторного автотрансформатора с возможностью регулировки мощности нагрева. Для снижения потерь теплоты поверхность стеклянной колбы была теплоизолирована пятью слоями вспененного полиэтилена, покрытого с одной стороны слоем фольги, который имеет общую толщину 13,75 мм и высоту ТАВ в твердом состоянии H, равную 43 мм, а в жидком, вследствие объемного расширения, составляющую 49 мм.

Z, мм

403,80

202,80 200,80

g

158,80 156,80

149,00 147,00

48,0 31,75

о irT

ОО

верт

Рис. 2. Габаритные размеры лабораторной модели:

— ось симметрии;

— стеклянный цилиндр;

— теплоизоляция;

— ТАВ;

— круглая алюминиевая пластина;

— нагреватель;

— монтажная пена;

— воздух;

— термопара

2,50 0

с-,

34,25

г, мм

В начальный момент времени ТАВ находится в твердом состоянии, близком к изотермическому, без наблюдаемых воздушных пузырей внутри объема. Для него после полного расплавления в ходе эксперимента теплоизоляция постепенно удалялась с нижней части колбы, чтобы процесс остывания проходил снизу вверх и растворенные газы успевали выйти наружу из объема ТАВ. Перед каждым экспериментом установка остывала до комнатной температуры (23.. .25 °С) в течение 24 ч, чтобы исключить дополнительный теплообмен в ходе эксперимента.

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

могут преодолеть силы вязкости, обусловленные пограничными слоями в жидкости. Оценить характер движения жидкого расплавленного ТАВ помогает использование при расчетах значений чисел Релея (Яа), Стефана (Б1е), Фурье (Бо) и Прандтля (Рг), а также такого параметра подобия, как соотношение сторон А, которые к концу процесса плавления представлены в табл. 1. Данные по температурам взяты из эксперимента.

Таблица 1

Значение безразмерных чисел подобия в конце процесса плавления для двух разных мощностей

Мощность, Вт Ra Ste Fo (4 j Pr A

7,14 2,7 • 106 0,139 0,94 18,75 1,21

34,23 — время 8,0 • 106 i плавления. 0,423 0,29 18,75 1,21

Числа подобия рассчитываются по следующим формулам:

Ra = giA^, Ste = ^, Fo(,) = ХжW, Pr = V, A = H, (1) V-жХж L ' V ' *ж W Хж W W

где g — ускорение свободного падения; Р — коэффициент теплового расширения; AT — вариация температуры; W — ширина (или диаметр); — кинематическая вязкость жидкой фракции ТАВ; хж — температуропроводность жидкой фракции ТАВ; cp ж — удельная

теплоемкость жидкой фракции ТАВ; L — скрытая теплота; t — время; H — высота колбы.

Согласно критерию преобладанию естественной конвекции и условию ламинарного течения [11], в процессе плавления конвекция считалась естественной после того, как будет достигнута определенная доля жидкости, при которой движение слоев жидкого ТАВ становилось ламинарным.

Условие преобладание естественной конвекции:

Ra

A

4 > 500; (2)

условие ламинарного течения:

Ra < 1010. (3)

Экспериментальный стенд был спроектирован таким образом, чтобы теплообмен с окружающей средой сводился к минимуму. Несмотря на это, для правильной настройки граничных условия моделей в ANSYS Fluent требуется осуществить оценку тепловых потерь. При оценке учитывались потери теплоты вследствие естественной конвекции с боковой поверхности теплоизоляции колбы и верхней торцевой поверхности ТАВ. Температура поверхности теплоизоляции измерялась телевизором (рис. 3). По ее полученным значениям в разные моменты времени выполнялся расчет потерь теплоты через эмпирические числа Нуссельта (Nu).

Рис. 3. Общий вид боковой поверхности колбы при измерении ее температуры телевизором

В соответствии с [11], для вертикально расположенного цилиндра большого диаметра среднее значение Киц оценивается по формуле

Nu ц =0,0246 (GrPr ))

1 Pr6

1 + 0,494Pr3

(4)

Для визуализации фронта расплава ТАВ использовалось неизолированное окно размерами 6x2 см, которое периодически закрывалось теплоизоляцией. Поскольку время просмотра профиля проплав-ления мало по сравнению с общим временем эксперимента, потери энергии при этом не учитывались.

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

2

температур, поэтому для такого случая, в соответствии с [12], число Нуссельта Nuд было рассчитано по формуле

Nu д =0,54 (GrPr )1/4. (5)

В формулах (4) и (5) число Грасгофа Gr определяется так:

Gr = Л-П, (6)

V 2

где l — характерный линейный размер поверхности теплообмена; Тп — температура поверхности.

В случае нагреваемого цилиндра l = 1верт (см. рис. 2), а для торцевой поверхности ТАВ размер l определяется как отношение площади поверхности к ее периметру.

Общие потери теплоты QIl в результате конвективного теплообмена можно рассчитать по формуле

1=п N11 • k i=n Nu • k

Qп (t) = EF; AT; At, + £ —U^F AT, At,, (7)

i=1 верт i=1 гор

где AT; — разница между температурой поверхности изоляции и температурой окружающей среды в течение интервала времени Ati ; F i — площадь нагретой поверхности в зависимости от момента времени t;; FЦ — площадь внешней поверхности теплоизоляции; F^ — верхняя площадь поверхности ТАВ; k — коэффициент теплопроводности.

На начальной стадии процесса (от 2 мин до 5 мин для различных мощностей нагревателя) большая часть тепловой энергии расходуется на нагрев установленной на нагревателе алюминиевой пластины Q& (t) и на нагрев монтажной пены Qп (t). В процессе проплав-

ления ТАВ энергия также тратится на нагрев корпуса из стекла Qс (t) толщиной 5, равной 2,5 мм (см. рис. 2). Тепловые потери энергии на нагрев колбы из стекла, алюминиевой пластины толщиной 7,8 мм и монтажной пены рассчитываются по формулам соответственно:

Qс (t) = Ссm (t)(ТЖ (t)-To), (8)

Qа (t) = cama (T2 (t)- To ), (9)

Qп (t) = cm (T2 (t)-To), (10)

где сс, тс (t) — теплоемкость и масса стеклянной колбы, контактирующей с расплавленным ТАВ; са, та — теплоемкость и масса алюминиевой пластины; сп, тп — теплоемкость и масса монтажной пены; Тж — приведенная средняя температура жидкой фракции ТАВ; Т — начальная температура; Т2 (t) — температура вблизи поверхности нагревателя.

Во время плавления ТАВ энергия тратится на нагрев (явная) Q^ (t)

и на фазовый переход (скрытая) QCK (t). Таким образом, полная поглощенная энергия ТАВ Q^^ (t), представляющая собой сумму скрытой и явной энергии, рассчитывается, согласно [13], по формуле

Q™™ (t) = Qяв (t)+Qск (t), (11)

где

Qяв (t)= Jp(T)c (T)( - To )) + (12)

Уж (t)

+ j p(T)c(Т)(Тж (t)-ТфП2) + j p(T)c(Т)(Ттв (t)-To))в (12)

МО ^тв (t)

и

Qcc (t) = РфП¿ФПУполн0 (t) . (13)

Здесь Уж (t) и Утв (t) — объемы жидкой и твердой фракции ТАВ соответственно; ТфП1 и Тш2 — начальная и конечная температуры плавления ТАВ соответственно; Тж (t) и Ттв (t) — приведенные средние температуры жидкой и твердой фазы соответственно; 7) — начальная температура; рФП — плотность ТАВ при температуре фазовых превращений (ФП); ЬФП — теплота ФП; Уполн — заполняемый объем с учетом расширения ТАВ в процессе нагрева; 0 (t) — доля расплавленного ТАВ; температурные зависимости плотности и теплопроводности р(Т) и с (Т) взяты из [14].

Экспериментальные зависимости (4) и (5) позволяют оценить коэффициенты теплоотдачи с поверхности теплоизоляции и наружной плоскости ТАВ, чтобы использовать их в качестве граничных условий для модели, построенной с помощью ANSYS Fluent. Были получены значения коэффициентов теплоотдачи с поверхности

теплоизоляции, равные 1,29 Вт / м2, и с верхней поверхности ТАВ, составляющие 9,01 Вт / м2. С помощью формулы (11) можно оценить поглощенные потери энергии, а также тепловые потери, которые идут на нагрев конструкции, и сравнить их значения с полученными путем численного расчета.

Теплофизические свойства материалов конструкции лабораторной модели приведены в табл. 2. В качестве ТАВ был использован октадекан (C18H38) [10, 14, 15], так как его температура плавления (28,18 ° С) входит в температурный диапазон работы большинства лазерных устройств и другой электронной техники.

Таблица 2

Теплофизические свойства материалов корпуса макета и окружающего воздуха

Параметр Материал

Стекло Теплоизоляция Алюминий Силикон Пена монтажная Воздух

Плотность р, кг/м3 2500 30 2719 1300 45 1,225

Изобарная удельная теплоемкость ср, Дж/(кг • К) 841 1,95 871 1420 1470 1006

Коэффициент теплопроводности к, Вт/(м • К) 1 0,0037 222,4 0,08 0,02 0,0242

Динамическая вязкость ц, Па • с - - - - - 1,7-10-5

Численная модель. Авторами были разработаны три численные модели с различными граничными условиями и размерностью. В модели I использовался двумерный подход к моделированию и не учитывались внешние потери теплоты при нагреве и плавлении ТАВ. Модель II от модели I отличается лишь учетом теплопотерь на нагрев конструкции и конвективным теплообменом с окружающей средой. В модели III использован трехмерный подход к моделированию и учитывались потери теплоты в окружающее пространство. Для всех случаев применялись приближение Буссинеска и постоянные свойства ТАВ. Постоянная тепловая мощность подводилась от нагревателя (см. рис. 2), ее соответствующие дискретные значения равны 7,14 Вт и 34,23 Вт.

Для всех моделей были приняты следующие упрощения: во-первых, течение в жидкой фазе ТАВ было несжимаемым и ламинарным; во-вторых, вязкая диссипация в жидкой фазе ТАВ была незначительной.

Математическое описание задачи. В случае, когда на процесс теплопередачи накладывается естественная конвекция для жидкой фазы ТАВ, процесс теплоотдачи определяется не только тепловыми, но и гидродинамическими явлениями. Их совокупность описывается системой дифференциальных уравнений, в которую входят уравнение теплопроводности, уравнение движения и уравнение непрерывности [16].

Упрощенная модель Буссинеска предполагает включать в уравнение импульса линеаризованную плавучесть, зависящую от температуры. Согласно [17], уравнение непрерывности с постоянной плотностью имеет вид

V-u = 0, (14)

где u — вектор скорости.

Уравнение движения с членом плавучести Fb и членом источника импульса Su можно представить следующим образом:

р fu + р (u -V)u = |V2u -Vp + Fb + Su. (15)

Gt

Здесь Vp — изменение давления; Su — источник импульса в уравнении движения; член плавучести в приближении Буссинеска определяется по формуле

Fb = еЩТ - Тпл) ez, (16)

где р — средняя приведенная плотность; ß — коэффициент теплового расширения; Тпл — температура плавления.

Плавучесть зависит от ß и разности температур. Член плавучести, входящий в уравнение движения (15), позволяет учитывать естественную конвекцию в жидкой фазе ТАВ.

Закон сохранения энергии для удельной энтальпии h с источником энтальпии Sh имеет вид

Gh

р — + р V - (uh) - V - (kVT) = Sh. (17)

Доля расплавленного ТАВ, обозначенного 9, вводится следующим образом:

0 =

0, если T < Тфщ;

T ТФПХ ТФП2 - ТФПХ

если Тфп < Т < Тфп2; (18)

если Т > Тфп2.

Как отмечалось выше, Тфщ и Тш2 означают нижний и верхний

пределы температуры плавления ТАВ соответственно. В уравнении (18) значения 0 = 0 и 0 = 1 представляют собой твердую и жидкую фазы ТАВ соответственно, а значение между 0 и 1 указывает на смесь двух фаз.

Приближение Буссинеска. Добавим еще несколько замечаний о приближении Буссинеска. Во-первых, поскольку уравнения (14) и (15) — нелинейные, они часто могут оказаться неустойчивыми и склонными к расхождению, если попытаться решить их численно. Во-вторых, требования к памяти довольно высоки, так как решаются три (или два) уравнения импульса для двух (или трех) компонент поля скоростей. Необходимо также решить уравнения для плотности и давления. Поэтому, для того чтобы улучшить стабильность этих уравнений, следует придать некоторую форму для аппроксимации плотности в уравнениях Навье — Стокса.

Если предположить, что изменение плотности достаточно мало, то можно заменить полиномиальную зависимость плотности от температуры простой суммой эталонной плотности Р0 и некоторой небольшой флуктуации Ар, так что Ар ^ Р0:

Р = Р0 +Ар. (19)

Поэтому после постановки (19) в уравнения Навье — Стокса флуктуация плотности Ар будет влиять только на член плавучести (16). Для того чтобы заменить зависимость р (Т) в уравнении (16)

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

в = -± if 1—1 f* (20)

Р0 VdTJp Р0 Т-Т0

где Т0 — опорное значение температуры.

В этом случае зависимость р (Т) с точки зрения эталонной плотности Р0 будет выглядеть следующим образом:

Р(Т) = Р0 (1 -Р(Т-Т0)). (21)

Подстановка (21) в выражение силы плавучести (16) позволяет избавиться от явной зависимости плотности р(Т). Следовательно,

приближение Буссинеска позволяет не записывать и не хранить информацию о плотности для каждой ячейки вычисленной сетки.

Принципы моделирования процесса плавления ТАВ. В качестве опорных значений постоянной средней плотности р и средней удельной теплоемкости Ср выбраны значения вблизи точки плавления. Для коэффициента теплопроводности к было принято значение в жидком состоянии кж, так как температурные градиенты и теплопроводность имели место в жидкой фазе. Данные с датчиков температур показывают, что твердая фаза была практически при постоянной температуре на протяжении всего процесса плавления. Динамическая вязкость д выбрана постоянной, она оценивалась по средней температуре жидкого ТАВ. Значения теплофизических свойств выбирались исходя из температурных градиентов для каждой мощности нагрева, полученных при моделировании модели I и из экспериментальных данных.

Теплофизические свойства октадекана, используемые в процессе моделирования, приведены ниже:

Мощность нагревателя, Вт...................................... 7,14 34,23

Плотность р, кг/м3............................................... 820

Удельная теплоемкость Ср, Дж/(кг • К)...................... 2100

Теплопроводность в жидком состоянии кж, Вт/(м • К) .... 0,15

Динамическая вязкость д, мПа • с.............................. 4,82 3,40

Коэффициент теплового расширения р, 1/К................ 8,9 • 10-4

Начальные и граничные условия для численных моделей.

Для сравнения скоростей расплавления ТАВ лабораторного эксперимента с расчетными данными были разработаны три численные модели. Габаритные размеры для моделей были выбраны согласно приведенным на рис. 2. С учетом объемного расширения ТАВ при нагреве было принято среднее значение 46 мм по высоте, соответствующее половинному расплавлению ТАВ.

Двумерные модели I и II представлены на рис. 4, а и б соответственно. Их различие заключается в конструктивных особенностях, а именно расчет для модели II выполнялся с учетом теплопотерь, которые идут на нагрев материалов конструкции, помимо нагрева самого ТАВ.

Модель I не учитывала также потери теплоты за счет излучения и конвекции с внешних поверхностей в окружающее пространство: стенки выбирались адиабатическими. Для модели II тепловые граничные условия учитывали теплообмен с окружающей средой. Значения коэффициентов теплоотдачи взяты из табл. 2.

к= 1,3 Вт/(м -К)

¿=1,3 Вт/(м -К)

Адиабатические

к =9,0 Вт/(м -К)

2 мм

2 мм

Q

Тепловой поток

Тепловой поток б

Рис. 4. Начальные и граничные условия для моделей I и II при То = 24,5 °С; Up = 0; P0 = 0

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

Размеры области плавления выбраны одинаковыми. Инициализация начальных условий ТАВ выбиралась следующим образом: все градиенты скоростей и давлений равнялись нулю, начальная температура всех частей для каждой модели принималась равной 24,5 °С. Начальные и граничные условия для каждой расчетной модели представлены в табл. 3.

Таблица 3

Граничные и начальные условия для трех расчетных моделей

Модель Граничные условия Начальные условия

I (2D) Стенки адиабатические: Q (t, x, у) = 0, при (x, у) е Г и Vt T0 (t = 0, x, у) = 24,5 °С, grad U0 (t = 0, x, у) = 0, grad p0 (t = 0, x, у) = 0, / ч ** где ( x, у) е Q

II (2D) Конвективный теплообмен с поверхности: - теплоизоляции, к = 1,3 Вт/(м2 • K); - ТАВ при контакте с воздухом, к = 9,0 Вт/(м2 • K)

III (3D) *Г(х, у) **Q — Конвективный теплообмен с поверхности: - теплоизоляции, к = 1,3 Вт/(м2 • K); - ТАВ при контакте с воздухом, к = 9,0 Вт/(м2 • K) - граница ТАВ (см. рис. 4). область внутри ТАВ (см. рис. 4).

Настройка ANSYS Fluent. Связь между давлением и скоростью осуществляется с помощью алгоритма PISO, а PRESTO применяется для интерполяции давления. Расчеты выполняются с использованием ПО ANSYS Fluent [18]. Сходимость достигается тогда, когда невязка уравнений энергии, импульса и неразрывности уменьшается до зна-

—15 _ 8 —8

чений менее 10 , 10 и 10 соответственно. Шаг по времени был установлен равным 0,1 с. Число итераций на шаг по времени установлено равным 10. Моделирование проводилось на ПК с процессором Intel Core i7 с частотой 3,4 ГГц и с 32 Гб ОЗУ.

Результаты сравнения. Сопоставим результаты численного расчета с экспериментальными данными, а именно сравним зависимости объемной доли расплавленного ТАВ для процесса плавления. Анализ проводился при мощностях нагрева 7,1 Вт и 34,2 Вт для трех численных моделей.

10 20 30 40 50 60 70 80 Время, мин а

10 15 Время, мин б

Рис. 5. Доли жидкой фазы для трех численных моделей и эксперимента при мощности 7,1 Вт (а) и 34,2 Вт (б): ----модель I;----модель II;-— модель III; — эксперимент

Зависимости долей расплавленного ТАВ для различных мощностей представлены на рис. 5. Характер движения границы фазового перехода имеет линейный характер для каждого случая, средние скорости плавления для моделей I и II обеих мощностей схожи. Различие заключается во временной задержке, равной 96 с для случая 34,2 Вт, и 228 с для случая 7,1 Вт, обусловленной нагревом конструкции в начальный момент времени и конвективным теплообменом с воздухом в конце процесса плавления.

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

Таблица 4

Средняя скорость движения границы фазового перехода и время расплавления

Вариант Средняя скорость, мм/мин Отличие модели от эксперимента, %

34,2 Вт

Модель I 2,52 9,1

Модель II 2,43 5,2

Модель III 3,03 31,2

Эксперимент 2,31 -

7,1 Вт

Модель I 0,62 5,4

Модель II 0,61 6,4

Модель III 0,65 0,3

Эксперимент 0,65 -

Обсуждение результатов. Большая расходимость при сравнении средних скоростей движения границы фазового перехода для мощности нагрева 34,2 Вт показывает, что выбранная для описания процесса плавления с большими градиентами температур модель с приближением Буссинеска дает значительное отклонение от эксперимента. Возможно, модель УОБ, учитывающая объемные расширения ТАВ, даст лучшую верификацию с экспериментальными данными, однако для нее потребуется больше вычислительной мощности.

Для меньшей мощности нагрева, составляющей 7,1 Вт, результат расчета лучше верифицируется с экспериментом, его отличие можно отнести к погрешности измерений.

При использовании приближения Буссинеска было сделано предположение, что оно справедливо в том случае, когда изменение плотности достаточно мало, т. е.

Ар<<р0, или Р-Р0 = в(Т-Т0)<< 1, (22)

Р0

где Т, Т0 — начальная и конечная температуры в процессе плавления.

При мощности нагрева 7,1 Вт разность этих температур составляет (Т - Т0 ) = 9 ° С, а при 34,2 Вт — (Т - Т0 ) = 45 °С. Следовательно,

Ар будет 0,8 % и 3,9 % от Р0 для случая 7,1 Вт и 34,2 Вт соответственно.

Видимо, коэффициент отклонения флуктуации у = (Ар / Р0) х х 100 % <1 % является границей применимости упрощенной модели с приближением Буссинеска. Можно ожидать, что для у > 1 % подробная модель с методом VOF будет более точной.

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

В работе [9] высказывалось предположение о том, что приближение Буссинеска будет давать некачественные результаты при больших градиентах температуры из-за учета явной зависимости плотности р(Т) в уравнениях движения жидкости, что и было получено

в представленной работе.

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

Заключение. Сравнение экспериментальных и расчетных данных показало различия в средних скоростях движения границы фазового перехода только при высокой мощности нагрева — 34,2 Вт. Однако для меньшей мощности, равной 7,1 Вт, расхождение с численным моделированием можно отнести к погрешности измерений. Несоответствие для большей мощности нагрева вызвано значительным изменением плотности ТАВ в начале и конце процесса плавления. Это заставляет рассматривать для процесса плавления с высокими градиентами температур модель VOF вместо упрощенного приближения Буссинеска. Сделан вывод, что при соответствующих настройках моделирование процесса плавления ТАВ для накопления скрытой теплоты может быть выполнено с помощью программного обеспечения ANSYS Fluent.

ЛИТЕРАТУРА

[1] Евдокимов Р.А., Тугаенко В.Ю., Щербенко Н.В. Перспективы применения и отработка технологии беспроводной передачи электрической энергии

между космическими аппаратами. Инженерный журнал: наука и инновации, 2022, вып. 7 (127). DOI: 10.18698/2308-6033-2022-7-2196

[2] Sharma A., Tyagi V.V., Chen C.R., Buddhi D. Review on thermal energy storage with phase change materials and applications. Renewable and Sustainable Energy Reviews, 2009, vol. 13 (2), pp. 318-345.

[3] Zalba B., Marin J., Cabeza L., Mehling H. Review on thermal energy storage with phase change: materials, heat transfer analysis and applications. Applied Thermal Engineering, 2003, vol. 23 (3), pp. 251-283.

[4] Asyraf W.M., Vasu A., Hagos F.Y., Noor M.M., Mamat R. Transient modelling of heat loading of phase change material for energy storage. MATEC Web of Conferences, 2017, vol. 90. pp. 1-12.

[5] Bondareva N.S., Sheremet M.A. Natural convection melting influence on the thermal resistance of a brick partially filled with phase change material. Fluids, 2021, vol. 6 (7), p. 258.

[6] Vogel J., Bauer D. Phase state and velocity measurements with high temporal and spatial resolution during melting of n-octadecane in a rectangular enclosure with two heated vertical sides. International Journal of Heat and Mass Transfer,

2018, vol. 127, pp. 1264-1276.

[7] Sharifi N., Robak C.W., Bergman T.L., Faghri A. Three-dimensional PCM melting in a vertical cylindrical enclosure including the effects of tilting. International Journal of Heat and Mass Transfer, 2013, vol. 65, pp. 798-806.

[8] Kamkari B., Shokouhmand H., Bruno F. Experimental investigation of the effect of inclination angle on convection-driven melting of phase change material in a rectangular enclosure. International Journal of Heat and Mass Transfer, 2014, vol. 72, pp. 186-200.

[9] Faghri A., Zhang Y. Fundamentals of Multiphase Heat Transfer and Flow. Springer, 2020.

[10] Vogel J., Thess A. Validation of a numerical model with a benchmark experiment for melting governed by natural convection in latent thermal energy storage. Applied Thermal Engineering, 2019, vol. 148, pp. 147-159.

[11] Batchelor G. K. Heat transfer by free convection across a closed cavity between vertical boundaries at different temperatures. Quarterly of Applied Mathematics, 1954, vol. 12, pp. 209-233.

[12] Wong H.Y. Heat transfer for engineers. Longman London and New York, 1979, 212 p.

[13] Shokouhmand H., Kamkari B. Experimental investigation on melting heat transfer characteristics of lauric acid in a rectangular thermal storage unit. Experimental Thermal and Fluid Science, 2013, vol. 50, pp. 201-212.

[14] Faden M., Hohlein S., Wanner J., Konig-Haagen A., Bruggemann D., Review of thermophysical property data of octadecane for phase-change studies. Materials,

2019, vol. 12 (2974), pp. 1-23.

[15] Galione P.A., Lehmkuhl B.O., Rigola S.J., Oliva L.A. Fixed-grid numerical modeling of melting and solidification using variable thermo-physical properties — Application to the melting of n-Octadecane inside a spherical capsule. International Journal of Heat and Mass Transfer, 2015, vol. 86, pp. 721-743.

[16] Михеев М.А., Михеева И.М. Основы теплопередачи. Москва, Энергия, 1977.

[17] ANSYS 15 Fluent Theory Guide. ANSYS, Inc, 2015, 814 p.

[18] Vikas A.Y., Vadav A., Soni S.K. Simulation of melting process of a phase change material (PCM) using ANSYS (Fluent). International Research Journal of Engineering and Technology (IRJET), 2017, vol. 4 (5), 7 p.

Статья поступила в редакцию 31.03.2023

Ссылку на эту статью просим оформлять следующим образом: Воропаев Р. А., Мацак И. С. Валидация численной модели плавления октадекана с помощью лабораторного эксперимента. Инженерный журнал: наука и инновации, 2023, вып. 8. http://dx.doi.org/10.18698/2308-6033-2023-8-2296

Воропаев Роман Андреевич — инженер II категории ПАО «РКК «Энергия». e-mail: roman.voropaev@rsce.ru

Мацак Иван Сергеевич — канд. техн. наук, ведущий научный сотрудник ПАО «РКК «Энергия». e-mail: Ivan.Macak@rsce.ru

Validation of the octadecane melting numerical model using a laboratory experiment

© R.A. Voropaev, I.S. Matsak

S.P. Korolev Rocket and Space Public Corporation Energia, Korolyov, 141070, Russia

The paper presents results of studying the melting process of octadecane, the heat-accumulating substance used in a space heat accumulator. The Stefan problem formulated for the phase transition boundary has an analytical solution only for the one-dimensional case. And in the case of the structure real geometry, it could be solved only using the numerical methods, in particular, by simulation in the ANSYS Fluent universal software system for the finite element analysis. To validate the developed numerical models, model experiments were carried out on octadecane melting in the cylindrical glass flask with heat-insulated side walls, the heater was positioned at the flask lower end. Two-dimensional and three-dimensional numerical models were investigated at the heating powers of 7.1 W and 34.2 W. The molten substance fraction dynamics in the two-dimensional and three-dimensional models was compared to the experimental data. Heat losses in heating the structure and convective heat exchange with the external media were taken into account. Results of numerical and laboratory experiments are presented, which could be used in testing the constant power heating model of a heat-accumulating substance. The results obtained show that with absolute temperature alteration in the heat-accumulating substance by an insignificant value (up to 10°C), the Boussinesq approximation provides results well verified by the experimental data.

Keywords: heat accumulator, melting, simulation, Boussinesq approximation, ANSYS Fluent

REFERENCES

[1] Evdokimov R.A., Tugaenko V.Yu., Scherbenko N.V. Perspektivy primeneniya i otrabotka tekhnologii besprovodnoy peredachi elektricheskoy energii mezhdu kosmicheskimi apparatami [Application and development prospects for spacecraft-to-spacecraft wireless transmission of electric energy]. Inzhenerny zhurnal: nauka i innovatsii — Engineering Journal: Science and Innovation, 2022, iss. 7 (127). https://doi.org/10.18698/2308-6033-2022-7-2196

[2] Sharma A., Tyagi V.V., Chen C.R., Buddhi D. Review on thermal energy storage with phase change materials and applications. Renewable and Sustainable Energy Reviews, 2009, vol. 13 (2), pp. 318-345.

[3] Zalba B., Marin J., Cabeza L., Mehling H. Review on thermal energy storage with phase change: materials, heat transfer analysis and applications. Applied Thermal Engineering, 2003, vol. 23 (3), pp. 251-283.

[4] Asyraf W.M., Vasu A., Hagos F.Y., Noor M.M., Mamat R. Transient modelling of heat loading of phase change material for energy storage. MATEC Web of Conferences, 2017, vol. 90, pp. 1-12.

[5] Bondareva N.S., Sheremet M.A. Natural convection melting influence on the thermal resistance of a brick partially filled with phase change material. Fluids, 2021, vol. 6 (7), p. 258.

[6] Vogel J., Bauer D. Phase state and velocity measurements with high temporal and spatial resolution during melting of n-octadecane in a rectangular enclosure with two heated vertical sides. International Journal of Heat and Mass Transfer, 2018, vol. 127, pp. 1264-1276.

Validation of the octadecane melting numerical model using a laboratory experiment

[7] Sharifi N., Robak C.W., Bergman T.L., Faghri A. Three-dimensional PCM melting in a vertical cylindrical enclosure including the effects of tilting. International Journal of Heat and Mass Transfer, 2013, vol. 65, pp. 798-806.

[8] Kamkari B., Shokouhmand H., Bruno F. Experimental investigation of the effect of inclination angle on convection-driven melting of phase change material in a rectangular enclosure. International Journal of Heat and Mass Transfer, 2014, vol. 72, pp. 186-200.

[9] Faghri A., Zhang Y. Fundamentals of Multiphase Heat Transfer and Flow. Springer, 2020.

[10] Vogel J., Thess A. Validation of a numerical model with a benchmark experiment for melting governed by natural convection in latent thermal energy storage. Applied Thermal Engineering, 2019, vol. 148, pp. 147-159.

[11] Batchelor G. K. Heat transfer by free convection across a closed cavity between vertical boundaries at different temperatures. Quarterly of Applied Mathematics, 1954, vol. 12, pp. 209-233.

[12] Wong H.Y. Heat transfer for engineers. Longman, London and New York, 1979, 212 p.

[13] Shokouhmand H., Kamkari B. Experimental investigation on melting heat transfer characteristics of lauric acid in a rectangular thermal storage unit. Experimental Thermal and Fluid Science, 2013, vol. 50, pp. 201-212.

[14] Faden M., Höhlein S., Wanner J., König-Haagen A., Brüggemann D. Review of thermophysical property data of octadecane for phase-change studies. Materials, 2019, vol. 12 (2974), pp. 1-23.

[15] Galione P.A., Lehmkuhl B.O., Rigola S.J., Oliva L.A. Fixed-grid numerical modeling of melting and solidification using variable thermo-physical properties — Application to the melting of n-Octadecane inside a spherical capsule. International Journal of Heat and Mass Transfer, 2015, vol. 86, pp. 721-743.

[16] Mikheev M.A., Mikheeva I.M. Osnovy teploperedachi [Basics of heat transfer]. Moscow, Energiya Publ., 1977.

[17] ANSYS 15 Fluent Theory Guide. ANSYS, Inc., 2015, 814 p.

[18] Vikas A.Y., Vadav A., Soni S. K. Simulation of melting process of a phase change material (PCM) using ANSYS (Fluent). International Research Journal of Engineering and Technology (IRJET), 2017, vol. 4 (5), 7.

Voropaev R.A., Engineer of the II category, S.P. Korolev Rocket and Space Public Corporation Energia. e-mail: roman.voropaev@rsce.ru

Matsak I.S., Cand. Sc. (Eng.), Leading Researcher, S.P. Korolev Rocket and Space Public Corporation Energia. e-mail: Ivan.Macak@rsce.ru