CC BY
66
УДК 536.24:674.047
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ИНТЕНСИВНЫХ ПРОЦЕССОВ ТЕПЛОМАССОПЕРЕНОСА ПРИ КОНВЕКТИВНОЙ СУШКЕ ДРЕВЕСИНЫ О.Р. Дорняк, В.А. Шамаев
Представлена математическая модель процессов тепломассопереноса в древесине, сформулированная на основе механики гетерогенных систем. Древесина рассматривается как трехфазная ненасыщенная коллоидная капиллярно-пористая среда. Проведено численное исследование модели для случая интенсивной термообработки древесины. Рассмотрены особенности развития полей теплофизических переменных
Ключевые слова: математическое моделирование, тепломассоперенос, древесина
Введение. Большое число физических и структурных факторов, определяющих результат интенсивного теплового воздействия на древесину, делает недоступным исследование теплофизических явлений путем прямого физического эксперимента. Процессы тепло- и массопереноса сопровождаются фазовыми превращениями на поверхности и внутри материала. Древесина имеет сложную капиллярнопористую структуру, а древесинное вещество обладает коллоидными свойствами. Развитие неоднородных полей влагосодержания и температуры ведет к росту напряжений и деформаций, что оказывает влияние на внутренний тепло- массоперенос и качество высушенного материала.
К ограничениям известных исследований по математическому моделированию тепломас-сопереноса в древесине, в частности, [1], следует отнести подход к древесине как к гомогенной среде, одномерный характер изучаемых процессов и недостаточный учет термодинамических особенностей поверхностнокапиллярных явлений. Для изучения тепло- и массопереноса в древесине, имеющей нерегулярную структуру, методологической основой математического моделирования явлений может быть механика многофазных систем [2-3]. В данной работе представлена математическая модель теплофизических процессов в древесине, как в трехфазной системе, с учетом различных механизмов переноса воды, а также анизотропии ее структурных свойств. На основе модели с использованием численных методов проведено исследование характера развития полей макроскопических теплофизических переменных в условиях интенсивного обезвоживания древесины.
x 1
Л
a ■yil
X 2 У
d
Дорняк Ольга Роальдовна - ВГЛТА, д-р техн. наук, ст. науч. сотрудник, e-mail: olga@,dorn.vrn.ru Шамаев Владимир Александрович - ВГЛТА, д-р техн. наук, профессор, e-mail: [email protected]
Рис. 1. Расчетная схема
Математическая модель. Пусть древесный образец имеет форму бруска прямоугольного сечения. Материал считается трансвер-сально-изотропным. Волокна параллельны одной из сторон бруска, например х3 (рис. 1). Увлажненный образец помещен в паровоздушную среду с постоянными значениями температуры Тс, давления рс и влагосодержания фс. Начальная влажность может быть распределена произвольно по сечению образца. Жидкая фаза практически несжимаема, ее плотность постоянна. Влага из образца удаляется преимущественно в виде пара. Перенос пара осуществляется в основном вдоль волокон и
1 3 2 3 3 3
V << у1 ,у1 << у1 . При этом у1 = у1 (х1, х2, х3). Здесь V - скорость. Верхние индексы 1,2,3 обозначают компоненты векторов, нижние - принадлежность к фазам. Индекс 1 относится к газообразной фазе, 2 - к жидкой, 3 - к твердой.
Математическая формулировка задачи включает уравнения, записанные для значений переменных, усредненных по объемам фаз (1 -газообразная, 2 - жидкая, 3 - твердая). Знак усреднения ниже опущен. Т.е. под переменной /2 следует подразумевать </2>2.
Газовая фаза состоит из двух компонент -неконденсирующегося газа и водяного пара. Примем предположение об идеальности паровой и газовой составляющих. Параметры, относящиеся к первой компоненте, обозначим нижним индексом 1§, ко второй - ^. Плотность р, концентрация %, парциальное давление р и удельная внутренняя энергия и компонент парогазовой смеси определяются следующим образом:
p° =p°v+p°g; z=;: -c=—г;
ri
pig
(i)
р
р* = р°т,\; ulg = ^;
Plv = рь^ь; ulv = cvlvTl. (2)
Здесь знак ° означает истинное значение физической величины; В - индивидуальная газовая постоянная, Дж/(кг К); су - теплоемкость при постоянном объеме, Дж/(кг К); T - температура, К.
Давление в парогазовой смеси определяется законом Дальтона для смеси идеальных газов:
Pl = р^Т1В1; Bl =%Б1у + (1 -х^ . (3)
Значения скорости паровой и газовой компонент могут быть различны. Для их описания введены среднемассовая скорость смещений элементарных макрообъемов первой фазы v1 и диффузионные скорости пара и газа ^1у и w1g:
v =cviv+(1 -z)vig; wig = vig - v; wiv = viv- v.
(4)
Относительное движение компонент определяется законом бинарной диффузии Фика:
Pio n dZ .
= ^-D^; w3v = -^-D^,. (5)
Pio n dZ
dx,
Р СЦ рь ^3
Параметр Б - коэффициент бинарной диффузии, зависящий в общем случае от температуры газа.
Уравнения сохранения массы для парогазовой смеси и газовой компоненты при сделанных предположениях имеют вид:
С(р°ад С
Э(рГ«1(1 -х))
- + VpOayx = s12j ;a + a2 + a3 = 1;
(6)
dt
-+v pa (1 - z)(v1 + wg) = о.
(7)
где а - объемное содержание; у -поток массы пара, обусловленный фазовыми переходами, отнесенный к единице времени и единице площади, кг/(м2'с); 512 - удельная поверхность раздела 1 и 2 фазы, м-. Случай у>0, соответствует испарению, у<0 - конденсации.
Уравнение движения и теплопроводности парогазовой фазы при сделанных предположениях записаны следующим образом:
р°а ^1 = а¥п ат^1 • (о)
ра — __а№- __ (8)
Ср,Р,а, ^=у(а1д,ут1)+а1в1т1 р ;
Ш Ш ; (9)
^1512у(Т11^12 - Т1) + 012 + 013
ср1 = хсрК + (1- х)ср1е; Qji = (т |ьу - Т) ;
1,] = 1,2,3. (10)
Здесь К3 - коэффициент проницаемости
третьей фазы первой при полном насыщении (если указаны верхние индексы, то в соответствующем направлении, например, в направлении т: К3тт), м2; ¥(0) - относительная фазовая проницаемость; 01 - насыщенность объема по-рового пространства газообразной фазой; ср -теплоемкость при постоянном давлении, Дж/(кгК); 1 -коэффициент теплопроводности, Вт/(м К); а] - коэффициент теплоотдачи между фазами 1 и у, Вт/м2. Переменные с индексом Еу относятся к границам раздела фаз 1 и у.
Для описания поведения жидкой фазы использованы уравнение сохранения массы, уравнение движения вязкой жидкости, без учета конвективного переноса, массовых сил и динамических эффектов фазовых переходов, а также уравнение теплопроводности [3]:
сра)
dt
' + Vp2a2v2 = s12j ;
о* = m2 (V jv2k + Vkv2); k, j = 1,2;
o dv2 a2m2v2
p°°a2 —- = -a2Vp2------------2^2 2
22 dt 2 2 K23Y2(q2)
dT2
+Va2o2;
c2p2a2 df _ V(a212VT2 ) c2s12 j(T2|s12 T2) +
(11)
(12)
(13)
(14)
+ бі2 + Qз2
Характеристики переноса воды зависят от типа ее связи с твердой фазой. В рассматриваемых условиях термического воздействия на древесный образец дополнительно используются усредненные уравнения баланса масса и количества движения связанной воды в смачивающих пленках воды, аналогичные по форме уравнениям (11-14) [4].
Для твердой фазы уравнение теплопроводности представлено в виде:
c3p3°a3 ^ = Q13 + Q23 + V (a313VT3 ) .
(15)
Функцию распределения давления в жидкой фазе р2 определим, следуя работе [5]. Неравновесный процесс сушки рассматривается как квазиравновесный, когда каждый локальный макрообъем пористого тела проходит через непрерывный ряд мгновенных состояний термодинамического равновесия между фазами. Из равенства химических потенциалов жидкости и пара в состоянии равновесия давление жидкости можно определить по формуле Кельвина [4]:
Р2 = Ps
t (T) + RT- ln -V2 P.
P1v
t (T)
(16)
где v2 - объем моля жидкой фазы, м3/моль; R -универсальная газовая постоянная,
Дж/(моль K), индекс sat относится к состоянию
3
насыщения.
Используя формулу Кельвина (16) и уравнение изотермы сорбции в форме
ръ
(р=-
- = F (W, T),
(17)
р*аг(Т )
имеем зависимость давления воды от влажности Ж и температуры Т в рамках равновесной термодинамики двухфазных многокомпонентных систем.
В среде с двойной пористостью, если радиус макропор много больше радиуса микро-пор, в процессе обезвоживания влагосодержа-ние в элементарных объемах системы различно, и характерное время процесса сушки может оказаться сопоставимым со временем установления полей капиллярного и расклинивающего давления. В этом случае процесс не является равновесным. Для учета неравновесных эффектов приняты различными значения температуры фаз гетерофазной системы. В качестве температуры Т в выражениях (16) и (17) используется температура жидкости в поверхностной фазе Е12.
Математическая модель включает также-уравнения сохранения на межфазных поверхностях. На границе раздела жидкость-пар (в поверхностной фазе Е12) в общем случае следует учитывать неравновесность фазовых переходов, связанную с тем, что количество пара, испаряющегося с поверхности раздела фаз зависит от кинетических возможностей паровой фазы. Кинетика неравновесных фазовых переходов описывается уравнением Герца-
Кнудсена-Ленгмюра:
J =
к
V2P1V
t (T2 Х1?)
Plv (T1 S12)
І
2 ІЛ2
(18)
Здесь к — коэффициент аккомодации.
Взаимосвязь между давлением и температурой вдоль линии насыщения определяется уравнением Клапейрона-Клаузиса:
dPi
dT1
Lpi
Ti(1 -PJ P2)
(19)
Здесь L - удельная теплота фазовых переходов,
2/2 м /С .
Неравновесная схема фазовых переходов предполагает наличие скачка температур в граничном кнудсеновском слое пара [2]
T 1 = T2\si2 - T1
0.45TM,J'
ІХ12
(20)
л/2В1уТ^р1 Уравнение сохранения энергии в поверхностной Е12 фазе без учета тепловой инерции Е12 фазы и влияния искривленности межфазных
поверхностей на ее температуру можно представить в виде
S12jL = а1 S12S12 (T1 — T1 |S12) + a2S12 S12(T2 — T2 | S12) • (21)
Условия теплообмена на остальных поверхностях раздела фаз могут быть записаны в виде балансовых соотношений
s13a3 S13(TS13 — T3) = —s13a1 S13
(Ts13 - T1);
s23a3 S23(TS23 T3) = S23a2S23(TS23 T2) .
(22)
3) °23^2£23(т£23 Т 2) . (23)
Краевые условия на внешних границах бруска
для пара:
дЛ
дп
= bG (р1| r - р >;-f
а I Г =1 -«| Г -«| Г
• 3 dT
;1 иг
для жидкости:
= 0; v.
= 0; -3 dT.
2 г = 0 ; 32 -і
г дп
= g (c|-c);
= aG (T1I Г - Tc); (24)
= аГ (T21 г - Tc);
«21 Г = 0.01W (Tc, j )а3Р3°/р
для тве 3 dT,
АГ
рдой фазы:
= «Г (T3I - Tc).
(25)
(26)
Здесь ау - коэффициент теплоотдачи 7-ой фазы к окружающей среде, Вт/м2; ^Г, ууГ — коэффициенты массообмена 7-ой фазы с окружающей средой, м~1. Нижний индекс Г означает принадлежность к внешним границам древесного образца, индекс с - к окружающей среде.
Начальные условия имеют вид:
V1 = 0 ; а1 = а10 ; X = X 0 ; а 2 = а 20 (Х1, Х2 , Х3 ) ;
v2 = 0, V22 = 0; Т1 = Т2 = Т3 = Тз ( х, , Х2, Х3). (27)
Поставленная начально-краевая задача (127) является нелинейной. Ее исследование может быть проведено только с помощью численных методов. Проведена аппроксимация дифференциальной задачи разностной с использованием метода контрольного объема [5]. Вычислительный алгоритм предусматривает три внутренних итерационных процесса для уравнений переноса в газообразной, конденсированной и твердой фазах, а также один внешний итерационный цикл. Для проверки работы программы была разработана система тестовых задач, имеющих точное аналитическое решение. Блоки уравнений переноса в каждой фазе тестировались отдельно.
Результаты расчетов имеют удовлетворительное совпадение с данными эксперимента [6].
Численные расчеты проведены для образ-
г
г
г
v
2
г
Г
цов из древесины березы размером 25x25x125 мм. Начальная влажность образца задана неод-
нородной по сечению в соответствии с выражением Ж0(х1,х2)=34+328Іп(лх1/а)8Іп(лх2/Ь).
а)
б)
м.%
м. %
в)
г)
Рис. 2. Распределение влагосодержания в центральном поперечном сечении образца (Х3//хар=0.5) в начальный момент времени - а) и при сушке древесины березы и для г/гхар=2 - б)- г). Температура сушильного агента Тс=363К - б), 433К - в), 463К - г)
Данное распределение показано на рисунке 2а. Начальная температура материала -770=77с=293К. Давление паровоздушной среды в сушильной камере равно атмосферному ^с=1.013'105 Па. Расчеты проведены для жестких режимов, когда температура мокрого термометра составляет Тм=357К, а температура сшильного агента варьировалась в диапазоне от 393К до 463К. Условия теплообмена на внешних границах образца заданы для отдельных фаз следующим образом
аГ = 25Вт /(м2 • К),
а3Г = 0.125Вт /(м2 • К),
аГ = 0.0Вт/(м2 • К). Жидкая фаза практически не обменивается теплом с внешней средой, ее температурный режим определяется внутрен-
оо
ним теплообменом с соседними фазами.
В расчетах использовались значения следующих постоянных параметров.
Водяной пар
^=8.3144 Дж/(моль К); ср 1у=2.034 103 Дж/(кг'К); су 1у=1.58103 Дж/(кг К); сУ 1е=7.07‘ 102 Дж/(кг К);
тт' 33-і /л-10 2
л.13 =10 м ; ср 1е=9.05‘ 102 Дж/(кгК);
= 1 • 10 -4 Вт /(м2 • К);
К).
Вода: р2°=103 кг/м3;
у2 =18' 10-6 м3/моль;
и
1 £12 и
1 £13
= 3.3 • 10-4Вт/(м2
и воздух:
Б1у=461.9 Дж/(кг К); 11=0.0248 Вт/(м К);
т1=10-6 Па' с; Б1у=284 Дж/(кг К); р1Г=2104 м“1;
Г і -1
У1 =1 м ;
^2=10 Па с; £=0.072 Н/м.
м. %
м. %
12=0.640 Вт/(м К); С2=4.2' 103 Дж/(кг К);
к=0.04; К2311=210-10 м2;
а= 6-10-5Вт/(м2 • К); Ь=2.25 1 06 м2/с2;
а^ = 0.0017Вт /(м2 • К).
Твердая фаза: р3°=1540 кг/м3;
а^13 = 1.7 10-5Вт/(м2 • К);
13=0.3 Вт/(мград); с3=3.7103 Дж/(кгК);
а1й^13 = 0.0010Вт/(м2 • К).
Характерные значения: 4ар=^, рю=0.7485
кг/м , pxap=pс, Тхар 373 K, ^ар
1хар2/(11/(рюСр1))/^. В расчетах, ^=0.74 (¿хар~15 мин). Для определения изменяющихся значений удельной поверхности раздела фаз в древесном образце 512, 523, 513 сформирована расчетная схема с использованием данных [7]. При этом начальное значение 512=6.7106 м"1.
‘«ха,
Рис. 3. Изменение со временем температуры фаз древесины березы 1 - газообразной, 2 - жидкой, 3 - твердой при сушке в среде с Тс=363К - кривые со штрихом и Тс=463К - без штриха
‘«ха.
Рис. 4. Изменение со временем избыточного давления парогазовой смеси Др=р,-рс и влажности Ж при Тс=463К - 1; 433К - 2, 393К - 3, 363К - 4. Кривые без штриха относятся к центральной точке образца, со штрихом - к точке вблизи поверхности (расположение точек показано на схеме)
Обсуждение результатов. На рисунках
2 - 6 представлены результаты расчетов процессов высокотемпературной сушки образцов, изготовленных из древесины березы. На рисунке 2 б-г показано поле влагосодержания по истечении промежутка времени ¿=24ар для различных значениий температуры воздуха в камере (температура мокрого термометра фиксирована).
Градиенты влажности, формирующиеся в поперечном сечении образца в результате интенсивной термообработки (рисунок 2 в-г) существенно выше, чем в мягких режимах сушки (рисунок 2 б). Рост температуры среды в камере интенсифицирует процесс обезвоживания образцов, так как приводит к снижению относительной влажности паровоздушной смеси, обеспечивая большую разность между парциальным давлением паровой компоненты и давлением насыщенного пара при данной температуре.
Результаты вычислительного эксперимен-
та подтверждают предположение о необходимости использования трехтемпературной модели для прогнозирования процессов тепло- и массопереноса при интенсивной термообработке древесины. Из графиков на рисунке 3 видно, что температура жидкой, газообразной и твердой фаз существенно различна в переходном процессе. Рост температуры жидкой фазы происходит менее интенсивно, чем твердой, и, тем более, газообразной фазы. Это отличие обусловлено существенно различной скоростью переходных тепловых процессов в твердой, жидкой и газообразной фазах, невысокой интенсивностью процессов внутреннего теплообмена жидкой фазы, ограничением температуры жидкой фазы температурой насыщенного пара при соответствующем давлении. Отметим, что распределение поля температуры каждой из фаз по объему материала близко к однородному, что связано, в том числе, и с достаточно малыми размерами образца. Расчетные температурные кривые для для точек в центре и на
периферии (см. схему к рис. 4) практически неразличимы.
При интенсивной термообработке древесины может происходить существенное повышение давления в парогазовой смеси. Из сравнения графиков на рисунке 4 следует, что при более жестком режиме сушки величина избыточного давления первой фазы возрастает, при этом удлиняется период его падения. Рост
давления в центральной точке образца выше, чем в приграничной. Это связано с большей насыщенностью жидкой фазой в данной точке, что может приводить, с одной стороны, к более высокой скорости производства паровой компоненты вследствие большей площади контакта первой и второй фаз, с другой стороны, к более высокому сопротивлению фильтрационного переноса паровоздушной смеси.
А Р/Р:
0.16
0.14
0.12
0.1
0.08
0.06
0.04
0.02
Рис. 5. Поле избыточного давления газообразной фазы Др1=р1(Х1,Х27Хз)-рс при сушке древесины березы в моменты времени г/гхар=0.05 - а), 2 - б) при Тс=463 К
На рисунке 5 представлено трехмерное поле давления в парогазовой фазе при ^хар=0.05 (момент времени, близкий к начальному) и ґ/ґхар=2 (конечный момент). В процессе сушки в каждой точке образца снижается как уровень избыточного давления, так и его градиент. Отвод пара из пористой системы материала обеспечивается градиентом давления парогазовой смеси преимущественно вдоль волокон. Соответствующее поле скоростей этой фазы иллюстрирует рисунок 6а.
Важной характеристикой внутренней паровоздушной среды, которая заполняет порис-
тую систему древесины и оказываетсуществен-ное влияние на кинетику ее обезвоживания, является концентия пара. Рисунок 6 иллюстрирует неравномерное распределение поля концентрации паровой компоненты в бинарной газовой смеси. Следовательно имеет место дополнительный фактор, который может ускорить процесс высушивания. Этот фактор -диффузионный перенос. Таким образом, только численное исследование конкретного процесса позволит оценить роль концентрации пара х как параметра, регулирующего производство пара и интенсивность диффузионного переноса.
Рис. 6. Поле скоростей парогазовой фазы у, (х1гх2,х3) - а) и концентрации пара %(х1гх2,х3) - б) при сушке древесины березы в момент времени г/гхар=2 при Тс=463К
Заключение. Проведенные расчеты показывают, что кинетика высокоинтенсивного процесса обезвоживания древесины в условиях конвективной сушки определяется сложным комплексом теплофизических факторов. Предложенная математическая модель позволяет оценить вклад каждого фактора в развитие пространственных полей влагосодержания, температуры и давления. Вычислительный эксперимент на основе модели позволяет подобрать эффективные режимы сушки древесины в зависимости от внешних условий.
Литература
1. Лыков, А.В. Теория сушки [Текст] / А. В. Лыков. - М. : Энергия, 1960. - 471 с.
1. Нигматулин, Р. И. Основы механики гетерогенных сред [Текст] / Р. И. Нигматулин. - М. :Наука, 1978. - 336 с.
2. Дорняк О.Р. Тепломассоперенос в ненасыщенных коллоидных капиллярно-пористых анизотропных материалах [Текст] / О.Р. Дорняк. - автореф. дис. ... д-ра техн. наук : 01.04.14 / Воронеж, 2007. - 32 с.
3. Дорняк, О. Р. Гидродинамическая задача для процессов модифицирования древесины [Текст] / О. Р. Дорняк // Известия Санкт-Петербургской лесотехнической академии. - 2005. - В. 172. - С. 143-150.
4. Гринчик, Н. Н. Процессы переноса в пористых средах, электролитах и мембранах [Текст] / Н. Н. Гринчик. - Минск: АНК «Институт тепло-и массообмена» АН Беларуси, 1991. - 251 с.
5. Патанкар, С. Численные методы решения задач теплообмена и динамики жидкости [Текст] / С. Па-танкар. - М. : Энергоатомиздат, 1984. - 152 с.
6. Шубин, Г. С. Сушка и тепловая обработка древесины [Текст] / Г. С. Шубин. - М. : Лесн. пром-сть, 1990. - 336 с.
Воронежская государственная лесотехническая академия
MATHEMATICAL MODELLING OF THE INTENSIVE HEAT AND MASS TRANSFER PROCESSES DURING CONVECTIVE DRYING OF WOOD
O.R. Dornyak, V.A. Shamaev
The mathematical model of heat and mass processes in the wood, formulated on the basis of the heterogeneous systems mechanics is presented. Wood is considered as the three-phase nonsaturated colloidal capillary-porous environment. The numerical research of model for a case of the intensive heat treatment of wood is conducted. The features of development of thermophysics variables fields are considered
Key words: heat and mass transfer, mathematical modeling
