Scientific approaches to solving the problem of joint processes of bubble boiling of freon and its movement in a coil

全文:

Универсальных методов расчета коэффициента теплоотдачи или универсального описания теплофизической картины кипения в потоке жидкости не существует, но создано множество эмпирических методов для различных хладонов в условиях невысоких значений давления (до 2…4 атм) и массовых расходов, а также выведены специальные формулы расчета теплоотдачи на малые каналы (диаметром от 1 мм и выше) и на строго приведенные значения давления [1–5].

В качестве рабочего тела для исследования выбран фреон R407C. Его термодинамические свойства близки к распространенному R22, при этом он не содержит хлора, негорючий, нетоксичный. Особенностью при исследовании являлся способ нагрева фреона – дополнительный нагреватель в виде змеевика, установленного на тыльной стороне фотоэлектрической панели, рис. 1.

 

Рис. 1. Проектируемый экспериментальный стенд.

 

1. Принятые допущения:

– Парожидкостная смесь. Конденсат перед подачей в испаритель подвергается дросселированию. При этом наряду со снижением давления с достаточной точностью можно задаться постоянством энтальпии:

i₁ = i₂. (1)

Процесс дросселирования на lnP-i диаграмме для выбранного хладагента представлен на рис. 2 ниже в процессе 3–4.

 

Рис. 2. Цикл теплового насоса на ln(P)-i диаграмме.

 

Принимая необходимый перепад давления при известных температурах конденсации и испарения, возможно рассчитать степень сухости пара в точке 4, соответствующей состоянию фреона на входе в испаритель. Степень сухости по заданным условиям равна x = 0.4. Таким образом, в качестве обязательного граничного условия задается доля в 40% массового расхода для пара.

•  Кипение. Течение происходит в трубке диаметром 5 мм. Принимается, комбинация пузырькового и пленочного кипения. При кипении и испарении жидкости на стенках образуются пузырьки, которые быстро схлопываются и образуют тонкую пленку жидкости, когда пара становится достаточно много в сечении канала.
•  Режим течения. Такой режим течения соответствует форме смеси, в которой пар образует ядро потока, при этом фаза жидкости протекает в виде пленки на поверхности стенки. Течение парожидкостной смеси в эмульсионном режиме до испарителя, что характерно для начальной стадии парообразования. Ввиду увеличения расхода пара, что видно из диаграммы на рис. 2, до доли 0.4, пузырьки склонны к скоплению и слипанию в потоке в более крупные.

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

 

Рис. 3. (а) Эмульсионный режим течения парожидкостной смеси; (б) Доля жидкостной фазы в сечении трубки; (в) Кольцевой режим течения (пар внутри жидкости); (г) Дисперсный режим течения (капли внутри пара).

 

Интенсивное испарение жидкости по границам сухих пятен (центров парообразования) взято в качестве определяющего механизма пузырькового кипения. Расчет с учетом данного фактора объясняет чрезвычайно высокую интенсивность теплообмена при пузырьковом кипении. В результате рассмотрения физической модели и обобщения экспериментальных данных для расчета теплоотдачи принимается аналитическая формула:

$q_{кип}=3.43\cdot {10}^{-4}\frac{{\lambda }^2\Delta T_{Н}^3}{\nu \sigma T_{\text{Н}}}\left[1+\frac{\left(h_{\sqrt{}}-h_{\Delta }\right){\left(\Delta T_{\text{Н}}\right)}^2}{2R_i{T}_{\text{Н}}}\right]\left(1+\sqrt{1+800B}+400B\right),$ (2)

$B=\frac{\left(h_{\sqrt{}}-h_{\Delta }\right){\left({\rho }_{\sqrt{}}\nu \right)}^{\mathrm{1,5}}}{\sigma {\left(\lambda T_{Н}\right)}^{\mathrm{0,5}}}.$ (3)

где ρ_Г, ν, λ, и σ – параметры жидкого фреона; R_i – газовая постоянная для пара фреона; ΔT_н – перепад между температурой стенки и жидкостью; все параметры берутся при температуре насыщения.

Теплоотдача при кипении в потоке жидкости рассчитывается как комбинация конвективного теплообмена и теплообмена кипением с преобладанием последнего:

$q=q_{\text{конв}}+q_{кип}.$ (4)

Итоговое критериальное уравнение подобия принимает вид:

$Nu=\frac{RePr\left(\xi /8\right)}{1+900/Re+\mathrm{12,7}{\left(\xi /8\right)}^{\mathrm{0,5}}\left(P{r}^{\mathrm{0,67}}-1\right)}.$ (5)

1. Теоретический расчет

Расчет с допущениями состоит из двух этапов:

1) Определение количества теплоты, поступающего в систему с учетом потерь в окружающую среду;

2) Определение температуры фреона на выходе из испарителя.

Первую часть представляется в виде блок-схемы на рис. 4 по методу Гаусса (с вводом исходных данных, прямым и обратным ходом, выводом результатов). Как итог рассчитывается температура на поверхности ФЭП с погрешностью менее 0.1%. Результаты с температурой окружающей среды в 20°С и ветром 1 м/с сводятся в табл. 1.

 

Рис. 4. Блок-схема расчета по методу Гаусса.

 

Таблица 1. Расчетные данные

$t_{\text{пов1}}+273$

Величина			Единицы измерения	Результат
Наименование	Обозначение	Расчетная формула или способ определения
Тепловой поток	Qсол	Дано	Вт/м2	700
Теплопроводность кремния	λкр	Дано	Вт/(м2·К)	148
Поглощательная способность кремния	Акр	Дано	–	0.2
Степень черноты кремния	εкр	Дано	–	0.71
Толщина кремния	δкр	Дано	м	0.001
Коэффициент Стефана-Больцмана	σ	Дано	Вт·м-2·К-4	5.67·10–8
Теплопроводность стекла	λст	Дано	Вт/(м2·К)	0.96
Поглощательная способность стекла	Аст	Дано	–	0.1
Степень черноты стекла	εст	Дано	–	0.93
Толщина стекла	δст	Дано	м	0.002
Длина и ширина ФЭП	l	Дано	м	0.15675
Площадь ФЭП	F	l2	м2	0.025
Отражательная способность кремния	Rкр	Дано	–	0.33
Тепло, поглощенное ФЭП	Q	$Q_{\text{сол}}\cdot F\left(1-R_{\text{кр}}\right)$	Вт	11.524
Тепло, поглощенное кремнием	Qкр.погл	$Q_{сол}\cdot F\cdot A_{\text{кр}}$	Вт	3.44
Тепло, прошедшее через кремний	Qкр.прош	$Q_{\text{сол}}\cdot F\left(1-R_{\text{кр}}-A_{\text{кр}}\right)$	Вт	8.084
Тепло, поглощенное стеклом	Qст.погл	$Q_{\text{сол}}\cdot F\cdot A_{\text{ст}}$	Вт	1.72
Тепло, прошедшее через стекло	Qст.прош	$Q_{\text{сол}}\cdot F\left(1-R_{\text{кр}}-A_{\text{кр}}-A_{\text{ст}}\right)$ $Q_{\text{сол}}\cdot F\left(1-R_{\text{кр}}-A_{\text{кр}}-A_{\text{ст}}\right)$	Вт	6.364
Скорость воздуха	ω	Дано	м2/с	1
Теплопроводность воздуха	λ	Дано	Вт/(м2·К)	2.4·10–2
Вязкость воздуха	νвозд	Дано	м2/с	14.61·10–6
Число Прандтля для воздуха	Pr	Дано	–	0.67
Температура окружающего воздуха	tвозд	Выбирается	℃	20
Температура окружающего воздуха	Tвозд	$t_{\text{возд}}+273$	К	293
Температура наружной поверхности ФЭП	tпов1	Выбирается	℃	29.12
Температура наружной поверхности ФЭП	Tпов1	$t_{\text{пов1}}+273$	К	302.12
Число Рейнольдса	Re1, Re2	$\frac{\omega \cdot 1}{{\nu }_{\text{возд}}}$	–	1.073·104
Число Нуссельта	Nu1, Nu2	$0.67R{e}^{0.5}P{r}^{1/3}$	–	60.727
Коэффициент теплоотдачи	α1, α2	$Nu\frac{\lambda }{1}$	Вт/(м2·К)	9.298
Температура внутренней поверхности ФЭП	tпов2	$t_{\text{пов1}}-\left[Q_{\text{сол}}\cdot F\left(1-R_{\text{кр}}-A_{кр}-A_{ст}\right)\right]\times$$\left(\frac{{\delta }_{кр}}{{\lambda }_{кр}}+\frac{{\delta }_{ст}}{{\lambda }_{ст}}\right)$	°C	29.1
Потери сверху	qсверху	${\alpha }_1\cdot F\left(T_{\text{пов1}}-T_{\text{возд}}\right)$	Вт	2083
Потери излучением кремния	qизл.кр	${\epsilon }_{кр}\cdot \sigma \cdot F\left(T_{\text{пов1}}^4-T_{\text{возд}}^4\right)$	Вт	0.951
Потери снизу	qснизу	${\alpha }_2\cdot F\left(T_{\text{пов2}}-T_{\text{возд}}\right)$	Вт	2.083
Потери излучением стекла	qизл.ст	${\epsilon }_{\text{ñò}}\cdot \sigma \cdot F\left[{\left(T_{\text{пов2}}^{}+273\right)}^4-T_{возд}^4\right]+$$Q_{сол}\cdot F\cdot A_{кр}+Q_{сол}\cdot F\cdot A_{ст}$	Вт	6.399
Тепловые потери в системе	Qпот	$q_{\text{сверху}}+q_{\text{изл.кр}}+q_{\text{снизу}}+q_{\text{изл.ст}}$	Вт	11.516
Погрешность	ΔQ	$\frac{Q-Q_{пот}}{Q}100$	%	0.067

 

На основе расчета температуры стенки ФЭП в зависимости от двух входных аргументов: скорости ветра и температуры окружающей среды необходимо построить 3D-диаграмму. Результат представлен на рис. 5.

 

Рис. 5. Зависимость температуры стенки ФЭП от температуры окружающей среды и скорости ветра.

 

Полученный результат используется для расчета температуры на поверхности трубки испарителя:

$t_{\text{ст.исп}}=t_{\text{ст.}}-q\cdot \sum _{i=1}^n{R}_i,$ (6)

где t_ст.исп – температура стенки испарителя; t_ст.ФЭП – температура на поверхности ФЭП из расчета выше; q – тепловой поток через ФЭП; R – термическое сопротивление слоев (защитного стекла, кремниевых пластин и т.д.).

Температура на стенке испарителя является граничным условием для моделирования процесса кипения фреона. Для дальнейшего теоретического расчета она является проверочной величиной при расчете коэффициента теплоотдачи.

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

$\alpha =0.075\left[1+10{\left(\frac{{\rho }^{\text{'}\text{'}}}{\rho -{\rho }^{\text{'}\text{'}}}\right)}^{2/3}\right]{\left(\frac{{\lambda }^2}{\nu \sigma T_H}\right)}^{1/3}{q}^{2/3}.$ (7)

По справочникам определяется температура насыщения T_н при давлении p в точке 4 на диаграмме по рис. 2. Выполняется поиск физических свойств насыщенного жидкого фреона на полученную температуру по справочным данным [6]:

•  теплопроводность λ;
•  число Прандтля Pr;
•  динамическую вязкость μ;
•  коэффициент поверхностного натяжения σ;
•  плотность насыщенной жидкости ρ и пара ρ››

Рассчитывается число Рейнольдса при заданной скорости течения ω:

$Re=\frac{\omega \cdot d\cdot \rho }{\mu }.$ (8)

Число Нуссельта рассчитывается по формуле Петухова (подразумевается развитый турбулентный режим, который возникнет из-за наличия поворотов испарителя и высокой скорости течения хладагента):

$Nu=0.021R{e}^{0.8}P{r}^{0.43}\epsilon ,$ (9)

где ε – поправочный коэффициент, учитывающий нагрев жидкости от стенки.

Далее сравнивается конвективная теплоотдача и теплоотдача кипением путем сравнения соответствующих коэффициентов. Расчет коэффициента теплоотдачи для процесса вынужденного течения в трубке определяется следующим образом:

${\alpha }_{\text{конв}}=\frac{Nu\cdot \lambda }{d}.$ (10)

Коэффициент теплоотдачи для кипения фреона α_кип можно определить по уравнению, соответствующему геометрическому и теплофизическому подобию по формуле (7) при условии давления кипения выше 1 атм. Если α_кип / α_конв > 3, то коэффициент теплоотдачи при кипении в трубке принимается равной α_кип, в противном случае – α_конв.

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

$t_{\text{ст. исп}}=\left(T_{\text{í}}-273.15\right)+\frac{q}{{\alpha }_{\text{кип}}}+\frac{q}{2{\lambda }_{ст}}d\cdot \ln \left(\frac{d}{d+2\delta }\right),$ (11)

где λ_ст – теплопроводность стенки; δ – толщина стенки (d равняется внутреннему диаметру трубки).

Для определения температуры фреона на выходе из испарителя определяется количество теплоты на доведение пароводяной смеси до состояния насыщенного пара, а также на его перегрев:

$Q_{кип}=\frac{\omega \cdot \rho \cdot \pi \cdot d^2}{4}r\left(1-x\right).$ (12)

Длина участка испарения:

$l_{исп}=\frac{Q_{кип}}{q\cdot \pi \cdot d}.$ (13)

Оставшийся участок является перегревателем хладагента:

$l_{перегр}=l_{\text{общ}}-l_{исп}.$ (14)

Температура фреона на выходе составляет:

$t_{\text{фр.вых}}=\frac{4q\cdot l_{перегр}}{\omega \cdot \rho \cdot d\cdot c_p^{\text{'}\text{'}}}+T_{Н}-\mathrm{273.15.}$ (15)

где c_p'' – теплоемкость пара фреона.

2. Дифференциальные уравнения при решении численным методом

Уравнения энергии:

$\frac{\partial \left(\rho u_i\right)}{\partial t}+\frac{\partial }{\partial x_j}\left(\rho u_i\rho u_j\right)+\frac{\partial P}{\partial x_i}=\frac{\partial }{\partial x}\left({\tau }_{ij}+{\tau }_{ij}^R\right)+S_i,$ (16)

$\frac{\partial \left(\rho E\right)}{\partial t}+\nabla \left[V\left(\rho E+\rho \right)\right]=\nabla \left[k_{\text{эф}}\Delta T-\sum _j{h}_j{J}_j+{\tau }_{эф}V\right]+S_h,$ (17)

где $\frac{\partial \left(\rho E\right)}{\partial t}$ – частный дифференциал нестабильного теплопереноса; $\nabla \left[V(\rho E+\rho )\right]$ – слагаемое, учитывающее теплопроводимость; $k_{эф}\Delta T$ – слагаемое, учитывающее диффузию смесей; $\sum _j{h}_j{J}_j$ – сумма, учитывающая вязкостную диссипацию; $S_h$ – слагаемое, учитывающее энтальпию притоков/ стоков.

В общем виде используемые в моделях с одним или двумя дифференциальными уравнениями, уравнение переноса можно записать в следующем виде:

$\rho \frac{\partial \Phi }{\partial t}+\rho \overline{u_j}\frac{\partial \Phi }{\partial x_j}=P-D+\frac{\partial }{\partial x_j}\left[(\mu +{\Gamma }_{\text{'}})\frac{\partial \Phi }{\partial x_j}\right]+A.$ (18)

Расшифровка параметров для каждого вида уравнений представлена в табл. 2.

 

Таблица 2. Расшифровка параметров

Наименование параметра	Ф	P	D	Гф
Кинетическая энергия	k	${\tau }_y\frac{\partial {\overline{u}}_i}{\partial x_j}$	pε или β*pkw	$\frac{{\mu }_t}{{\sigma }_k}$
Скорость диссипации кинетической энергии	ε	$c_{\epsilon 1}\frac{\epsilon }{k}{\tau }_{ij}\frac{\partial {\overline{u}}_i}{\partial x_j}$	$c_{\epsilon 2}p\frac{{\epsilon }^2}{k}$	$\frac{u_t}{{\sigma }_e}$
Удельная скорость диссипации	ω	$\alpha \frac{\omega }{k}{\tau }_{ij}\frac{\partial {\overline{u}}_i}{\partial x_j}$	βpω2	σωut

 

В расчетах использована модель турбулентности Transition k-kl-ω.

Моделирование в программной среде ANSYS Fluent

Геометрическая модель, созданная с помощью компьютерного моделирования процесса течения фреона R407C в канале латунного испарителя с внешним диаметром 10 мм и толщиной стенки 1 мм при дисперсно-кольцевом течении, показана на рис. 6а.

 

Рис. 6. (а) Геометрия модели испарителя; (б) Сеточная структура трубки испарителя.

 

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

Моделирование включает в себя условие течения жидкости в поле тяжести с образованием паровых пузырьков. Группа задач воздействия греющей стенки на испаряемую жидкость с удовлетворительной точностью решается нестационарной постановкой при задействовании модели турбулентности Transition k-kl-ω. Это трехпараметрическая модель группы RANS, подходящая для узкого круга задач, позволяющая моделировать переход от ламинарного течения к турбулентному, в которой используются пристеночные функции. Основные решаемые уравнения, закладываемые в математическую модель: уравнение энергии и уравнения модели турбулентности и диссипации турбулентности.

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

$\frac{\partial \tilde{u}}{\partial t}+\left(\tilde{u}\cdot \nabla \right)\tilde{u}=-\frac{1}{{\rho }_0}\nabla p+v\nabla \tilde{u}-\beta \theta \overrightarrow{g},$ (19)

$\frac{\partial \theta }{\partial t}+\tilde{u}\cdot \nabla \theta =\chi \Delta \theta ,$ (20)

$\text{div}\tilde{u}=\mathrm{0,}$ (21)

$\rho \left(T\right)={\rho }_0(1-\beta \theta ),$ (22)

где $\tilde{u}$ – скорость течения; ρ₀ – плотность жидкости при некоторой равновесной температуре T₀; Θ – отклонение температуры от равновесного состояния (Θ = = T – T₀); p – давление; ν – кинематическая вязкость; χ − коэффициент температуропроводности; $\overrightarrow{g}$ − ускорение свободного падения; β − коэффициент объемного расширения.

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

На рис. 7 показаны основные результаты численного моделирования в пост-процессинге. На рис.7а видно, что из-за нестабильности процесса сходимость протекает медленно. Однако на рис.7б, в, г, д показаны итоги по температурам и долям жидкой и паровой фаз фреона.

 

Рис. 7. (а) Колебания невязок расчета дисперсно-кольцевого режима; (б) Схождение степени сухости на выходе из испарителя; (в) Перемешивание парожидкостной смеси; (г) Поле температур; (д) Доля пара в объеме.

 

Сравнение температуры на тыльной поверхности ФЭП полученные экспериментально, теоретически и численным моделированием показало расхождение данных оказалось приемлемым для инженерных методик. Наиболее близкие результаты получились при теоретическом расчете и численном моделировании. Разность составила 0.1°С. При сравнении с экспериментальным данными по термодатчикам – расхождение уже более значительное, примерно, 0.3°С. Если учитывать погрешность измерительных приборов – 0.5°С, но даже при таких значениях отклонений температуры необходимо понимать порядок измеряемых величин – от –10 до + 90°С. То есть при диапазоне измерений в 100°С погрешность оказывается незначительной.

Заключение

При оптимизации математической модели теплопередачи к трубке испарителя и ряда численных расчетов было установлено, что режим течения в испарителе дисперсно-кольцевой. Пар в эмульсионной смеси перераспределяется так, что основная фаза (жидкая с долей в 60% расхода) движется по поверхности трубки в виде мелких капель в паре – в кольцевом режиме (рис. 2в).

По мере движения двухфазного потока происходит смена на противоположный предельный случай, когда вдоль стенки протекает пар, а жидкость движется в виде мелких капель – в дисперсном режиме (рис. 3г). Влажный пар доходит до состояния насыщения, а затем сухой пар перегревается (процесс 4–1 на рис. 2).

Также выявлено, что для получения корректного коэффициента теплоотдачи фреона для испарителя типа ”змеевик“ необходим вывод критериального уравнения подобия на основании экспериментальных данных.

Результаты расчета испарителя с моделью турбулентности Transition k-kl-ω показали достаточно приемлемую для инженерных расчетов сходимость с результатами теоретического расчета и экспериментальными данными. Расхождение между значениями полученных температур составили от 0.1 до 0.5°С в измеряемом диапазоне температур 100°С.

Исследование выполнено за счет гранта Российского научного фонда № 22-19-20011.

作者简介

A. Kovalev

Federal State Autonomous Educational Institution of Higher Education “South Ural State University (National Research University)”

Email: korniakovaoi@susu.ru
俄罗斯联邦, Chelyabinsk

Ya. Bolkov

Federal State Autonomous Educational Institution of Higher Education “South Ural State University (National Research University)”

Email: korniakovaoi@susu.ru
俄罗斯联邦, Chelyabinsk

K. Osintsev

Federal State Autonomous Educational Institution of Higher Education “South Ural State University (National Research University)”

Email: korniakovaoi@susu.ru
俄罗斯联邦, Chelyabinsk

O. Kornyakova

Federal State Autonomous Educational Institution of Higher Education “South Ural State University (National Research University)”

编辑信件的主要联系方式.
Email: korniakovaoi@susu.ru
俄罗斯联邦, Chelyabinsk

参考

补充文件