Условия плавления тепловыделяющего пористого элемента с активным охлаждением: приближенное аналитическое решение

Полный текст

Введение

Переход к низкоуглеродной энергетике со значительной долей возобновляемых ресурсов требует подготовки больших емкостей для запасания тепловой и электрической энергии. К устройствам для запасания энергии предъявляется ряд требований, среди которых – устойчивость к перегрузкам и низкая скорость деградации. В связи с этими вопросами существенное значение приобретает разработка систем теплового контроля [1, 2]. В частности, необходимо установить предельные значения параметров, соответствующих границам эффективной и безопасной работы аккумуляторов. В настоящей работе рассматривается задача такого рода, а именно – определение критического тепловыделения, соответствующего началу плавления пористого материала и полному его расплавлению. Для этого решаются стационарные уравнения теплообмена в одномерных образцах (с плоской и цилиндрической симметрией) в различных приближениях.

Отметим, что подобные задачи появляются также в химической технологии, обработке материалов, теории горения [3–6]. Определение границы расплавления в пористой среде является важной технической задачей (например, при работе с промерзающими средами) [7]. Полученные результаты для постоянной интенсивности тепловыделения могут быть легко распространены на линейную зависимость тепловыделения от температуры; более сложные зависимости могут приводить к критическим явлениям другого рода [8].

Уравнение теплопроводности в пористой среде с тепловыделением и фазовым переходом

Начнем с рассмотрения уравнения теплопереноса в пористой тепловыделяющей среде при фильтрации несжимаемой жидкости. Примем приближение об изотропности этой среды и слабой зависимости проницаемости от других параметров при температурах ниже температуры плавления (U(r) = const). Тогда распределение скорости в каждом сечении будет однородным. Распределение температуры определяется из уравнения:

$c\rho U\frac{\partial T}{\partial z}=\lambda \left[\frac{1}{r^n}\frac{\partial }{\partial r}\left(r^n\frac{\partial T}{\partial r}\right)+\frac{{\partial }^{2}T}{\partial z^2}\right]+W$. (1)

Подобные задачи исследовались, например, в работе [9] при исследовании режимов охлаждения реакторов. Вопросы существования решения рассматривались в работах [10, 11]. При U = 0 задача (1) сводится к определению профиля температуры в тепловыделяющем образце с охлаждением на внешней границе: эта постановка исследовалась в работах [12–14]. Исследования устойчивости течения в пористых каналах с внешним нагревом проводились в работах [7, 15].

Далее, будем считать, что высота области невелика, а скорость течения достаточна для того, чтобы пренебречь переносом температуры по высоте по сравнению с поперечным переносом. Кроме того, используя характерный температурный перепад ΔT, перейдем к безразмерному уравнению:

$P{e}^{\text{'}}\frac{D}{H}\frac{\partial \theta }{\partial \eta }=\frac{1}{{\xi }^n}\frac{\partial }{\partial \xi }\left({\xi }^n\frac{\partial \theta }{\partial \xi }\right)+Q$. (2)

Здесь D и H – поперечный и продольный размеры области. Комбинируя их и заменяя в (2) производную по высоте конечной разностью, приходим к обыкновенному дифференциальному уравнению:

$\frac{1}{{\xi }^n}\frac{d}{d\xi }\left({\xi }^n\frac{d\theta }{d\xi }\right)-Pe\theta +Q=0$. (3)

Здесь число Pe включает отношение размеров D и H. Для образца с плоской симметрией n = 0, для цилиндрической симметрии n = 1.

Уравнение (3) определяет распределение температуры в тепловыделяющем пористом материале при его фильтрационном охлаждении с интенсивностью Pe. При заданной температуре оболочки граничными условиями для уравнения (3) будут следующие равенства:

$\frac{d\theta }{d\xi }\left(0\right)=0$; $\theta \left(1\right)=0$. (4а)

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

$\frac{d\theta }{d\xi }\left(0\right)=0$; $-\frac{d\theta }{d\xi }\left(1\right)=Bi\theta \left(1\right)$. (4б)

Если в пористом материале может происходить фазовый переход при температуре T_m, то ΔT = T_m – T₀, и в области с θ > 1 течение становится невозможным из-за коллапса пор за счет агломерации. Из работ по тепловой устойчивости тепловыделяющих элементов известно, что в стационарном состоянии формируется поверхность раздела фаз. Обозначим расстояние от оси симметрии до поверхности раздела ξ_m так что θ(ξ_m) = 1. Тогда параметр Pe зависит от переменной ξ, и в первом приближении является кусочно-постоянной функцией, так что:

$Pe\left(\xi \right)=\left\{\begin{array}{l}\mathrm{0,}\xi \le {\xi }_m\\ P{e}_0,\xi >{\xi }_m\end{array}\right.$. (5)

Уравнение энергии разбивается при этом на отдельные уравнения для фазы расплава и фазы оболочки:

$\frac{1}{{\xi }^n}\frac{d}{d\xi }\left({\xi }^n\frac{d{\theta }_s}{d\xi }\right)-P{e}_0{\theta }_s+Q=0$ , $\xi >{\xi }_m$,(6а)

$\frac{\kappa }{{\xi }^n}\frac{d}{d\xi }\left({\xi }^n\frac{d{\theta }_l}{d\xi }\right)+Q=0$, $\xi <{\xi }_m$. (6б)

Здесь κ – отношение теплопроводностей фаз. Тогда распределение температур в области ξ < ξ_m известно из работ [12, 13]:

${\theta }_l\left(\xi \right)=1+\frac{Q}{2\left(n+1\right)\kappa }\left({\xi }_m^2-{\xi }^2\right)$. (7а)

Далее необходимо достроить распределение, решив уравнение (3) с новыми граничными условиями. Эти условия определяются равенством температур и тепловых потоков на границе раздела фаз:

${\theta }_l\left({\xi }_m\right)={\theta }_s\left({\xi }_m\right)=1$, $\kappa \frac{d\theta }{d\xi }\left({\xi }_m^{-}\right)=\frac{d\theta }{d\xi }\left({\xi }_m^{+}\right)$. (7б)

Первая задача состоит в определении критического значения интенсивности тепловыделения при заданном потоке жидкости (либо критического значения потока жидкости при заданной интенсивности тепловыделения), т.е. такого Q, при котором в рассматриваемой области начинается расплавление материала (θ_max = 1), либо происходит полное расплавление материала (θ_min = 1). Для определения условий начала плавления достаточно решить однофазную задачу и найти условия достижения температуры плавления. Полное расплавление становится возможным только при конечной скорости теплообмена на внешней границе: как будет показано ниже, проточное охлаждение накладывает дополнительное ограничение на условие полного расплавления.

Условие начала расплавления

Решение для плоского канала

Для плоской симметрии уравнение (3) для однофазной области с заданной температурой границы (4а) может быть решено простым образом [16]:

$\theta =\frac{Q}{P{e}_0}\left[1-\frac{\cosh \left(\xi \sqrt{P{e}_0}\right)}{\cosh \left(\sqrt{P{e}_0}\right)}\right]$. (8)

Достижение температуры θ(0) = 1 соответствует началу плавления. При этом выполняется следующее соотношение между интенсивностью тепловыделения и числом Пекле:

$Q_{min}=\frac{P{e}_0}{1-\frac{1}{\cosh \left(\sqrt{P{e}_0}\right)}}$. (9)

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

$W=\frac{c\rho U\left(T_m-T_0\right)D^2}{L}Q$. (10)

Критическое значение интенсивности тепловыделения для конечных чисел Био (4б) может быть найдено в виде:

$Q_{min}=\frac{P{e}_0}{1-\frac{1}{\cosh \left(\sqrt{P{e}_0}\right)+\frac{1}{Bi}\sinh \left(\sqrt{P{e}_0}\right)}}$. (11)

При близких к нулю числах Пекле критическое значение интенсивности тепловыделения близко к 2 (это точное решение для задачи теплопроводности с объемным тепловыделением в [12]).

$Q_{min}\simeq 2+P{e}_0$. (12)

Решение для цилиндрического канала

Для цилиндрической симметрии уравнение (3) с условиями (4а) может быть представлено с помощью модифицированных функций Бесселя [16]:

$\theta =\frac{Q}{P{e}_0}\left[1-\frac{I_0\left(\xi \sqrt{P{e}_0}\right)}{I_0\left(\sqrt{P{e}_0}\right)}\right]$. (13)

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

$Q_{min}=\frac{P{e}_0}{1-\frac{1}{I_0\left(\sqrt{P{e}_0}\right)}}$. (14)

С граничными условиями (4б) получаем поправку на число Био:

$Q_{min}=\frac{P{e}_0}{1-\frac{1}{I_0\left(\sqrt{P{e}_0}\right)+\frac{1}{Bi}{I}_{-1}\left(\sqrt{P{e}_0}\right)}}$. (15)

Функция I₀ является растущей функцией модуля аргумента, поэтому в пределе больших числе Пекле зависимость в формуле превращается в прямую пропорциональность, как и для плоской симметрии. Это совпадение ожидаемо: при больших расходах конвективное охлаждение играет более важную роль по сравнению с поперечным переносом, поэтому геометрия области становится незначимой. Для малых чисел Пекле разложение до первых двух слагаемых дает приближение:

$Q_{min}\simeq 4+P{e}_0$. (16)

Условие полного расплавления

При полном расплавлении образца на внешней границе достигается температура фазового перехода. Поскольку при этом, согласно (5), число Пекле равно нулю во всей области, то условие полного расплавления для задачи с проточным охлаждением совпадает с условием расплавления образца с внешним охлаждением [13]:

$Q_{max}=\left(n+1\right)Bi$. (17)

Численные расчеты проводились с помощью оборудования ЦКП “Высокотемпературный контур” (ИСЭМ СО РАН). Для вычислений использовалась численная схема из работы [17]. Результаты численного моделирования показывают, что условие полного расплавления (17) действительно выполняется при малых числах Пекле. Однако при больших числах Пекле, а именно при

$P{e}_0>\left(n+1\right)Bi$, (18)

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

Постоянный расход охлаждающей жидкости

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

$Pe\left(\xi \right)=\left\{\begin{array}{l}\mathrm{0,}\xi \le {\xi }_m\\ \frac{P{e}_0}{1-{\xi }_m^n},\xi >{\xi }_m\end{array}\right.$. (19)

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

Расчетная зависимость координаты поверхности раздела фаз от интенсивности тепловыделения и числа Пекле для цилиндрической области приведены на рис. 1. Видны характерные скачки в области критических значений интенсивности тепловыделения, связанные с изменением теплового баланса: такие резкие переходы между режимами теплообмена напоминают критические явления в химических реакторах [8]. Учет кинетики фазовых переходов может привести к большему разнообразию особенностей поведения [18].

 

Рис. 1. Зависимость положения границы фазового перехода в тепловыделяющем цилиндре от интенсивности тепловыделения (Q) и интенсивности охлаждения (числа рядом с кривыми, Pe₀) для Bi = 20 и κ = 1: (а) постоянная скорость течения; (б) постоянный расход.

 

Рассмотренные одномерные задачи и их решения могут послужить основой для исследования двумерной задачи Стефана для тепловыделяющих пористых сред. Можно ожидать, что наблюдаемые в одномерном случае критические явления будут проявляться и при переходе к более реалистичным постановкам. При этом, однако, граница расплава будет уже зависеть от высоты слоя. При малых расходах и больших температурных перепадах становятся значимыми процессы естественной конвекции [19, 20]. Эти направления могут стать объектом следующих работ.

Заключение

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

Работа выполнена при поддержке Российского научного фонда, проект № 23-2900406 (https://rscf.ru/project/23-29-00406/).

Обозначения

Bi – число Био

D – поперечный размер области (полуширина, радиус), м

H – высота области, м

I_k – модифицированная функция Бесселя k-го порядка

Pe – число Пекле

Pe₀ – число Пекле в проницаемой части сечения

Q – безразмерная интенсивность тепловыделения

Q_min – безразмерная интенсивность тепловыделения, соответствующая началу плавления

Q_max – безразмерная интенсивность тепловыделения, соответствующая полному расплавлению области

T – температура

T_m – температура плавления

T₀ – температура охлаждения

U – скорость фильтрации, м/с

W – интенсивность тепловыделения, Вт/м³

c – теплоемкость, Дж/кг/К

n – показатель симметрии области

r – поперечная координата, м

z – продольная координата, м

κ – отношение теплопроводностей фаз

λ – коэффициент теплопроводности, Вт/м/К

ρ – плотность, кг/м³

θ – безразмерная температура

θ_l – температура в области расплава

θ_m – безразмерная температура плавления

θ_s – температура в пристенной области

η – безразмерная продольная координата

ξ – безразмерная поперечная координата

ξ_m – безразмерное положение фронта плавления

Об авторах

И. Г. Донской

Федеральное государственное бюджетное учреждение науки Институт систем энергетики им. Л.А. Мелентьева Сибирского отделения Российской академии наук

Автор, ответственный за переписку.
Email: donskoy.chem@mail.ru
Россия, Иркутск

Список литературы

Дополнительные файлы