CC BY
22
ПРОБЛЕМЫ АГРОПРОМЫШЛЕННОГО КОМПЛЕКСА
УДК 536.24:674.047
ДИФФУЗИОННЫИ ПЕРЕНОС В ПАРОГАЗОВОЙ ФАЗЕ ПРИ СУШКЕ ДРЕВЕСИНЫ. Ч. I. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ
© 2006 г. О.Р. Дорняк
Введение
Прогресс в теории и практике сушки требует дальнейшего изучения фундаментальных закономерностей тепломассообмена, гидродинамики, технологии процесса на базе современных методов анализа, использования методов физического и математического моделирования [1]. При изучении процессов сушки материалов, допускающих усадку, подверженных короблению и растрескиванию, важен расчет локальных значений влагосодержания и температуры, что возможно осуществить лишь путем создания и исследования математических моделей, описывающих внутренний тепло- и массоперенос. В [2-3] разработана математическая модель процессов тепло- и мас-сопереноса в древесине как анизотропной ненасыщенной трехфазной среде с учетом фазовых переходов и поверхностных явлений. Континуальные уравнения получены на основе методологии механики гетерогенных систем [4], позволяющей сформулировать единый подход к теоретическому исследованию и созданию математических моделей процессов гидротермической и механической обработки древесины.
В данной работе на основе базовой модели [2] изучаются особенности диффузионного переноса пара в парогазовой смеси, влияние внутренних и внешних условий массообмена на интенсивность процесса сушки.
Математическая модель
Уравнения сохранения для макроскопических параметров системы получены методом объемного усреднения по отдельным фазам микроуравнений для макроскопических параметров каждой фазы. Обозначим параметры переноса первой (газообразной) фазы нижним индексом 1, второй фазы (воды) - 2, третьей фазы (твердого скелета, древесинного вещества -комплекса природных полимеров) - 3. Знак усреднения по соответствующей фазе будет опущен.
Газовая фаза является двухкомпонентной гомогенной системой, содержащей неконденсирующийся газ и водяной пар. Параметры, относящиеся к первой компоненте, отмечены нижним индексом ко второй компоненте - Плотность парогазовой смеси р1
и концентрации составляющих определяются следующим образом:
„О п"
p l _р lv +p lg . X_ o . 1 X_ o P1 P1
(1)
где знак ° означает истинное значение физической величины; х - массовая концентрация пара.
Паровая и газовая компоненты рассматриваются как идеальные газы. Давление в смеси p1 определяется законом Дальтона
Р! =р М ; Б1 =х£1у + (1 - X)Б^ . (3)
Для парциального давления и удельной внутренней энергии и имеем
plg = рlgTlБlg ; = су;
Р\у =Р°УТ1Б1У ; и IV = су1ут1 , (4)
где В - индивидуальная газовая постоянная, Дж/(кгК); cv - теплоемкость при постоянном объеме, Дж/(кгК); Т - температура, К.
Значения скорости паровой и газовой компонент могут быть различны. Для их описания введены сред-немассовая скорость смещений элементарных макрообъемов первой фазы v1 и диффузионные скорости пара и газа w1v и w1g:
V: = XV¡V +(1 -Х)Vlg ; w1g = v1g - V];
^ = ^ - v1. (5)
Относительное движение компонент определяется законом бинарной диффузии Фика:
w3 _flDix . w,3v _-PlDix, (6) 1g P olg dx / 1v p o°v dx 3 W
где параметр В - коэффициент бинарной диффузии, зависящий в общем случае от температуры газа. Декартова координата х3 направлена вдоль волокон образца.
Уравнения сохранения массы для парогазовой смеси и газовой компоненты при сделанных предпо-
ложениях имеют вид
д(рОа,) д(рОау3) ,
-—— + —= s12j ; а 1 + а2 + а3 = 1;(7)
dt
- + -
dx 3
д(р Га 1(1 -X)) + д (Р 1а 1(1 -X)(v1 + )) = 0 (8)
dt
дх 3
где / -поток массы пара, обусловленный фазовыми переходами, отнесенный к единице времени и единице площади, кг/(м2с); 512 - удельная поверхность раздела 1 и 2 фазы, м-1. Случай/ > 0 соответствует испарению, / < 0 - конденсации.
Уравнение движения и теплопроводности парогазовой фазы записаны следующим образом:
Р1 а 1
ср1Р 1а 1
ду1 3 ду^
3
дt
- + у
дх 3
= -а
дР1 1 дх 3
а у
f дГ 3 дГ 1 — + у{—1
дt дх
= а 1В1Г1
K1333Т 1(91)
дрГ |у 3 др О
(9)
дt
-+ у
дх1
д f дГ1 ^ д f дГ1 ^ д
+- а 1—1
дх1 дх1
дх
f
+а ВТ
др О + у3 др О
1
дt дх1
а 1
дх 2
дх 3
. дГ1
а 1Л1—1 дх 3
+ cy1s1^.j'(r1| v12 - Г) + 021 + Ö31;
(10)
cp1 = Xcp1y +(1 -X) c
p1 g
= //(Ц^ -Тг); / = 1,2,3.
Здесь К13тт - коэффициент проницаемости 1-ой фазы при полном насыщении в направлении т, м2; Т(9) - относительная фазовая проницаемость; 91 -насыщенность объема порового пространства газообразной фазой; ср - теплоемкость при постоянном давлении, Дж/(кгК); X -коэффициент теплопроводности, Вт/(мК); а/ - коэффициент теплоотдачи между фазами г и /, Вт/м2; переменные с индексом Е/ относятся к границам раздела фаз г и/
Для описания поведения жидкой фазы использованы уравнение сохранения массы, уравнение движения жидкости без учета конвективного переноса, массовых сил и динамических эффектов фазовых переходов, а также уравнение теплопроводности. Характеристики переноса воды зависят от типа ее связи с твердой фазой. В рассматриваемых условиях термического воздействия на древесный образец дополнительно применены усредненные уравнения баланса массы и количества движения связанной воды в смачивающих пленках.
Для твердой фазы использовано уравнение теплопроводности.
Функция распределения давления в жидкой фазе определена, следуя работе [5]. Неравновесный про-
цесс сушки рассматривается как квазиравновесный, когда каждый локальный макрообъем пористого тела проходит через непрерывный ряд мгновенных состояний термодинамического равновесия между фазами. Из равенства химических потенциалов жидкости и пара в состоянии равновесия давление жидкости определяется по формуле Кельвина [5]. Используя формулу Кельвина и уравнение изотермы сорбции, можно получить зависимость давления воды от влажности и температуры в рамках равновесной термодинамики двухфазных многокомпонентных систем.
Уравнения переноса массы, количества движения и внутренней энергии отдельных фаз дополняют уравнения сохранения на межфазных поверхностях, записанные в виде балансовых соотношений. На границе раздела жидкость - пар (в поверхностной фазе) в общем случае следует учитывать неравновесность фазовых переходов, связанную с тем, что количество пара, испаряющегося с поверхности раздела фаз, зависит от кинетических возможностей паровой фазы. Кинетика неравновесных фазовых переходов описывается уравнением Герца-Кнудсена-Ленгмюра. Неравновесная схема фазовых переходов предполагает наличие скачка температур в граничном кнудсенов-ском слое пара, что может быть учтено в рамках предлагаемой математической модели. Взаимосвязь между давлением и температурой вдоль линии насыщения определяется уравнением Клапейрона-Клаузи-са.
В процессе сушки происходит изменение площадей поверхностей контактного взаимодействия трех фаз изучаемой гетерогенной системы Расчетная схема определения удельной площади поверхности раздела фаз жидкость - парогазовая смесь, жидкость -твердая фаза, парогазовая смесь - твердая фаза построена на базе устоявшихся в древесиноведении представлений о капиллярно-пористой структуре древесины и формах связи влаги с древесиной [6], а также опытных данных.
Система дифференциальных уравнений в частных производных дополнена начальными и граничными условиями. Объемная концентрация жидкой фазы в поверхностной зоне соответствует равновесной. Краевые условия на внешних границах бруска для температуры отдельных фаз, концентрации и давления пара записаны в форме условий 3 рода. В частности, для переменных, относящихся к первой, фазе граничные условия записаны в виде Эр!
дп
= ß Г (Р11Г - p c);
- D
дх
дп
= Y
ir (х| Г -X c):
(11)
а1| Г =1 -а 2 | Г -а3 |Г;
-А
ът1
дп
= а }(Г1| г - Г с),
(12)
г
г
г
где а1 — коэффициент теплоотдачи /-ой фазы к окружающей среде, Вт/м2; Р1Г, у1Г — коэффициенты массообмена 1-ой фазы с окружающей средой, м-1; п - внешняя нормаль к поверхности образца. Нижний индекс Г означает принадлежность к внешним границам древесного образца, индекс с - к окружающей среде.
Начальные условия для величины х предполагают равновесное распределение пара по объему образца х = Х 0(х 1,х 2,х 3), соответствующее имеющимся полям температуры Т1 = Т0( х1, х 2, х 3) и влажности материала. Паровоздушная смесь неподвижна в начальный момент времени V3 = 0. Ее объемное содержание определяется по известным формулам в зависимости от породы и влажности древесины а 1 = а 10(х1,х2,х3).
Задача эволюции концентрации паровой компоненты в процессе термического воздействия на образец, определяемая соотношениями (1)-(12), может быть приведена к безразмерному виду. В качестве единиц величин с независимой размерностью используем характерный размер заготовки /хар, характерную скорость течения vxap, давление рхар, температуру Тхар. Безразмерные переменные определим следующим образом:
Э(р O*a i (1 -X)) d(p O*a 2 (1 -X)(v 3* + w g*))
. /=1,2,3; t* =-L
. B* = _ Bi
xap
B
хар
* pn
T
'12 _ д 12'xap
р m . * ; Cv1 Cv1
р xap B xap
j ; t xap l xaр
j xap V xap
'1
dt *
rdv3* 3* dv3* ^
—l— + vf —L
dt dx 3
dx *
= 0;
3*
dp!* -__l__
dxз* Re 1 Da33 ^2(62)
P1 a 1 1
+-
Pe: 1
+Pe1
( dT-f 3* dT-f ^
—* + v3 —*
dt dx *
** = a! Bj Tj
f o* o* Л
Эр J + v 3* Эр J +
dt *
dx *
d ( dT*) d ( dT*) d ( dT*^ 1 1 a 1 1 1 ™ 1
dx1
a
dx*
dx
dx
dx *
a
dx *
•^Nu?m(T;|^ -T1*) + s^Nuf£13(71*1 ^ -T1*)
лЛ * . * x/T-T * I m * ч * * ч *
+ -üL s 12 j (T1 I 12 - T1 ); с p1 = Xе p1v +(1 -X)c p1 g
с p1 Im ' ^
dpL
dn *
оГ*, H *4 dx
= ßi (p1 I -p*);
T dn
= Nur
Лх| г -x с);
aJГ=1-a2 Г-a
dT1*
Г
dn *
T^ Л А
= Nur(T1 |Г -T* );
v13*(x1 *,х2 *,х3* ,0) = 0 ; р1(х1",х2 *,х3* ,0) = р*;
ТГ^*, х 2*, х 3*,0) = Т0*. (13)
Безразмерные комплексы, характеризующие теп-ломассоперенос в двухкомпонентной газообразной фазе, введены следующим образом:
p p V2
р = У xap . j = У хар . в = . хаР М xap 2 ' J xap v ' xap
Vx2
xaр
T
xap
Система безразмерных уравнений переноса парогазовой фазы, а также относящихся к ним граничных и начальных условий, имеет вид
* * *
О* o* , о* К 1v 1 ' 1g * * т> *
р 1 =р iv +р ig. x = -0*; 1 -x = 7^; u 1g = *v 1gT1 ; р1
p1g. рО
u 1v = cv *vt1*; p*g = POgT** b1*g; p*v = p 1vt* B*v;
p* =р OT*B*; B* = xB*v + (1 - X)B*g ;
v* = xv*v + (1 - X)v*g ; w*g = vjg - v*; w* = v* - v* ;
w3* 1 р0 dX . w3» 1 р0 dX .
W л —---Wi —----
1g Pe1D р О* dx **' 1v Pe1D р О* dx 3' ^O*a 1) + d(pO*a*v3*) = s1*2j*; a 1 + a2 + a3 = 1;
dt *
dx 3
N г Y^ NuJD =-
1D D
V l р
xap
Pe 1D =
l2
xap . ßr* =ßrl
d ' M1 _ M1 ' тар ■ t xap D
Re1 =
Nu
xap xa^^ xap
ц 1
a h 1
i Yij1 xap
Pe1 =
l
2
xap
А1
t xap a 1
Cp*P xap
i £ij
St
a
i £ij
i £ij
Nu
i £ij .
C iр xap Vxap
Pe i
arl
Nu[ = a 11 xap . -
А1
i-1,2; Dai7 =
nn
m = 1, 2, 3.
l2
xaр
Сформулированная модель диффузионного переноса пара (13) позволяет выделить ряд предельных режимов возможных в процессах сушки коллоидных капиллярно-пористых тел. Режимы движения бинарной газовой смеси определяются числами Рейнольдса Яе1, Дарси Ба13 и р1Г. Эффективность внутреннего теплообмена зависит от критериев Пекле Ре1 и Нус-сельта Ми1Л, внешнего - от числа Ми1Г. Режимы масо-обмена в парогазовой среде определяются диффузионными критериями Нуссельта Ми1Вг и Пекле Ре1В.
г
; v1 =
x=
г
l
t
xap
n
xap
h
h
А
При значениях №ш<<1 процесс сушки малоэффективен, концентрация пара в образце не может снижаться вследствие диффузионного механизма выравнивания концентраций пара в материале и в окружающей среде. При числах №ш>>1 концентрация пара практически мгновенно подстраивается под изменяющееся условия теплового взаимодействия материала с окружающей средой.
Длительность переходного диффузионного процесса в парогазовой фазе капиллярно-пористого материала зависит от диффузионного числа Пекле. При Рещ<<1 эта величина весьма невелика. При Реш>>1, поле концентраций пара в паровоздушной смеси по-рового пространства будет перестраиваться весьма медленно.
Учет диффузионного переноса пара в математической модели особенно важен в режиме №ш~1, Реш~1, когда вклад значений концентраций паровой компоненты в развитие полей других теплофизиче-ских переменных наиболее значителен.
Заключение
Сформулированная модель диффузионного переноса пара в парогазовой смеси, заполняющей поры капиллярно-пористого тела, позволяет с использованием численных методов изучить особенности пове-
дения бинарной газовой системы, влияние диффузии паров воды и воздуха на интенсивность процесса сушки материала. Предельные режимы внешнего и внутреннего массообмена могут быть использованы для упрощения математической модели.
Литература
1. Рудобашта С.П. Фундаментальные исследования тепломассообмена при сушке // Современные энергосберегающие тепловые технологии (Сушка и тепловые процессы СЭТТ-2005): Тр. 2 Международ. науч. практич. конф. - М., 2005. - Т.1. - С. 7-17.
2. Дорняк О.Р., Шульман З.П. Математическое моделирова-
ние процесса сушки натуральной и уплотненной древесины // Изв. вузов, Сев.-Кавк. регион. Техн. науки. 2005. Спец. выпуск «Композиционные и порошковые материалы» С. 133-142.
3. Dornyak O.R., Shulman Z.P. Mathematical modeling of two-dimensional field development of temperature and moisture content in wood // Int. J. Applied Mechanics and Engineering. - 2005. - Vol. 10, № 4. - P. 593-604.
4. Нигматулин Р.И. Основы механики гетерогенных сред. -М., 1978.
5. Гринчик Н.Н. Процессы переноса в пористых средах, электролитах и мембранах. - Минск: АНК «Институт те-пло-и массообмена» АН Беларуси, 1991.
6. Чудинов Б.С. Вода в древесине. - Новосибирск, 1984.
23 мая 2006 г.
Воронежская государственная лесотехническая академия
УДК 536.24:674.047
ДИФФУЗИОННЫЙ ПЕРЕНОС В ПАРОГАЗОВОЙ ФАЗЕ ПРИ СУШКЕ ДРЕВЕСИНЫ. Ч. II. ВЫЧИСЛИТЕЛЬНЫЙ ЭКСПЕРИМЕНТ
© 2006 г. О.Р. Дорняк
Численная реализация
Задача диффузионного переноса в двухкомпо-нентной газообразной фазе при тепловых воздействиях на древесину, сформулированная в первой части работы, решается в рамках общей проблемы тепло-массопереноса [1].
Численная реализация задачи проводится на основе нестационарных неявных конечно-разностных схем, построенных методом контрольного объема с использованием процедуры расщепления по физическим процессам. Линеаризация дискретных уравнений проведена методом Ньютона. В задаче движения парогазовой смеси величина давление p1 является искомой. Для расчета изменяющихся полей p1 использован специальный алгоритм SIMPLE [2]. Системы разностных уравнений решаются методом прогонки. Для проверки алгоритма проведена серия тестовых расчетов на основе специально подобранных точных решений исследуемых уравнений с источниковыми членами. Сравнение расчетных значений изменения со
временем избыточного давления в парогазовой среде с опытными данными [3] показало их качественное совпадение.
Изменение концентрации пара в газообразной фазе описано с помощью нестационарного уравнения параболического типа. Уравнение является одномерным, поскольку перенос этой фазы осуществляется преимущественно вдоль волокон, но концентрация х является функцией трех пространственных переменных, поскольку в соответствие с общей постановкой задачи, развитие поля концентраций связано с распределением температуры в парогазовой среде, а также зависит от интенсивности фазовых переходов / Так как характерное время выравнивания температуры в паре намного порядков меньше времени процесса и много меньше времени развития температурных полей в твердом скелете (древесинном веществе) и в жидкой фазе, то следует ожидать, что поля концентраций пара в паровоздушной среде будут близки к однородным в каждом поперечном сечении образца.
Расчетные значения теплофизических параметров пара, воздуха, воды и твердой фазы приведены в [1]. Образец представляет собой брус из древесины сосны с размерами 180x180x6000 мм. Температура окружающей среды Тс = 363 К. Температура мокрого термометра здесь Тм = 323 К. Давление паровоздушной смеси в окружающей среде равно атмосферному рс = 1,013105 Па. Величины характерных параметров приняты следующими: /хар = й = 6 м, ¿хар = 3,6103 с, Рхар = рс, Тхар = 373 К. Начальная влажность образца однородна Ж0 = 25 и 80 %. В первом случае процесс сушки происходит в основном в гигроскопической области, во втором - при значениях влажности вне этого диапазона (расчетное время процесса равно ¿хар). Начальная температура всех фаз одинакова и равна Тс. При / = 0 величина давления паровоздушной смеси в образце и в окружающей среде равны. Значения кон-
0 0,2 0,4 0,6 0,8 t/tx
(Л/Рхар^ИО
12 10 8 6
-2
0,4 0,6
в)
центрации пара и влагосодержания в окружающей среде и образце в начальный момент различны: Хс = 0,086, хо = 0,589. Эти величины рассчитаны по относительной влажности воздуха в окружающей среде и во влажном материале [4].
В построенной математической модели характер диффузионного переноса пара в бинарной системе определяют два безразмерных комплекса - диффузионное число Пекле Реш = / х^ар / / хар В и диффузионный
критерий Нуссельта Ки^ =у^/хар /В (верхний индекс Г относится к внешней границе образца; В - коэффициент взаимной диффузии паровоздушной смеси, м2/с; у1Г - коэффициент массообмена пара, м/с). Интервалы варьирования безразмерных параметров массо-обмена в расчетах составили Ре1В = 3,3-3,3*102, Ки1В = 0,02-2'102.
Ж, %
70
60
50
(///хар)" 10
0,6
г)
Рис. 1. Изменение со временем объемной концентрации пара в газообразной фазе в центральной точке бруса х (а); средней по сечению бруса влажности Ж (б); давления парогазовой смеси рх (в); интенсивности парообразования / (г) при однородной начальной влажности Ж0 = 80 % для Ре1В = 3,3; = 0,02 - 1; 0,2 - 2; 2 - 3; 20 - 4; 200 - 5
4
2
0
0,5
0,4
0,3
0,2
0,1
0 0,2 0,4 0,6 0,! а)
(Р1/Рхар-1)"10-
(///хар)" 107
20
15 / / V \ 2 -
10 / / 3
5 11 1 ^
5
0
4
-5
3,0
2,5 3
2,0
1,5 / 2 / -
1,0 -
0,5 / / 1 ——
0 5
-0,5 4
0
0,2
0,4
0,6
0J
t/t«
0
0,2
0,4
0,6
0,
t/tx;
в) г)
Рис. 2. Изменение со временем объемной концентрации пара в газообразной фазе в центральной точке бруса х (а); средней по сечению бруса влажности Ж (б); давления парогазовой смеси рх (в); интенсивности парообразования / (г); при однородной начальной влажности Ж0 = 25 % для Реш = 3,3; Киш = 0,2 - 1; 2 - 2; 20 - 3; а также Реш = 3,3'102; Киш = 20 - 4 и Реш = 33;
= 20 - 5
Обсуждение результатов Рис. 1-4 иллюстрируют влияние массообменных критериев на интенсивность процесса сушки древесного материала.
На рис. 1 и 2 показано изменение основных характеристик парогазовой смеси при сушке древесного бруса при влажности выше предела гигроскопичности Жпг (рис. 1) и в гигроскопической области (рис. 2).
Как видно из графиков на рис. 1а и 2а, концентрация пара х снижается более интенсивно при больших значениях диффузионного числа Нуссельта. Большим значениям числа Ки1В соответствует более высокая скорость обезвоживания материала (рис. 1б, 2 б). Более низкие значения концентрации пара х приводят к увеличению его производства из-за увеличения разности давления насыщения и парциального давления паровой компоненты (рис. 1г, 2г), вследст-
вие чего развитие переходного процесса в газообразной фазе происходит существенно быстрее при больших числах Ки1В (рис. 1в, 2в). В менее влажном материале длительность процесса выравнивания давления еще более сокращается из-за меньшего сопротивления оттоку парогазовой смеси из образца.
Из рис. 1а - г видно, что при значении Ки1В<<1 процесс сушки существенно замедляется после удаления воды с внешних поверхностей образца и установления на внешних границах равновесной влажности при данной температуре окружающей среды (кривые 1). Его ускорение возможно, например, за счет более интенсивных капиллярных перетоков связанной воды. При числах Ки1В>>1 изменение теплофизических переменных происходит в большем диапазоне. При этом имеет место некоторый предельный режим их изменения (кривые 4, 5). При Ки1В~1 процессы массо-
X
обмена на границе образца вносят наиболее существенный вклад в формирование полей основных теп-лофизических характеристик, влияя тем самым на интенсивность процесса его сушки.
X 1
0,6 -
0,5 v6 ■
0,4 ■ 2
0,3 / 3 \
0,2 ___ \ -
0,1 5
0
0 0,2 0,4 0,6 0,8 x3/1хяр а)
X 0,6 0,5 0,4 0,3 0,2 0,1
расчетного периода обуславливает ту же тенденцию и для парциального давления пара, что приводит вновь к росту производства пара.
При больших числах Реш практически не работает механизм внутреннего диффузионного переноса пара в парогазовой среде (кривая 4 на рис. 2б, в). Даже при достаточно благоприятных условиях внешнего масо-обмена (Киш = 20), процесс сушки идет в этом случае значительно медленнее, чем при малых диффузионных числах Пекле. При Реш = 333 в центре образца в течение расчетного периода происходят процессы конденсации (/ < 0). Средняя влажность материала снижается, как показывают результаты расчетов, в основном, за счет испарения влаги в областях, близких к торцевым поверхностям.
При одинаковых условиях внешнего массообмена (одинаковых значениях Киш) характер распределения концентрации пара в газообразной фазе х вдоль волокон образца качественно одинаков для значений ^о>Жпг и Ж0<Жпг (см. рис. 3). Кривые 4 и 6 на рис. 3а, а также 3-5 рис. 3б, относящиеся к расчету процесса сушки материала при малых и больших величинах Ре1В, демонстрируют вклад внутреннего диффузионного переноса в формирование профиля концентрации пара вдоль направления его преимущественного движения. При больших значениях Ре1В к завершению расчетного периода концентрация пара вдоль центрального волокна практически равна начальной, а вблизи торцевых поверхностей образца, имеют место достаточно большие градиенты величины %.
X
0,6
Рис. 3. Распределение объемной концентрации пара в газообразной фазе вдоль длинной оси симметрии бруса х для момента времени г/гхар = 1 при начальной однородной влажности Ж0 = 80 % (а) 25 % (б). Кривые на рис. 3а получены для Ре1В = 3,3; №1В = 0,02 - 1 ; 0,2 - 2; 2 - 3; 20 - 4; 200 - 5 и Ре1В = 33; Ш1В = 20 - 6; на рис. 3б для Ре1В = 3,3; №1В = 0,2 - 1; 2 - 2; 20 - 3; а также Ре1В = 3,3102, №1В = 20 - 4 и Реш = 33; Ш1В = 20 - 5
Вне гигроскопической области качественная картина изменения со временем давления и интенсивности фазовых переходов в выбранной точке объема совпадают (рис. 1 в, г). Представляет интерес, что в гигроскопической области на фоне роста давления в газовой фазе интенсивность фазовых переходов изменяется со временем немонотонно при больших диффузионных числах Нуссельта (рис. 2в, г). Падение величины / связано с уменьшением давления насыщенного пара вследствие снижения значения активности пара из-за потери влаги пористым материалом. Существенное уменьшение давления р1 в середине
0,5
0,4
0,3
0,2
0,1
Рис. 4. Изменение объемной концентрации пара в газообразной фазе вдоль длинной оси симметрии бруса х при начальной однородной влажности Ж0 = 80 % - 2; 2'; 2м; 25 % - 1, 1', 1м для Ре1В = 3,3; Ки1В = 20. Кривые относятся к моменту времени г/гхар = 0,035 - 1-2; 0,5 - 1'-2'; 1,0 -Г-2"
При одинаковых параметрах внутреннего и внешнего массообмена профиль концентрации пара вдоль длинной стороны образца развивается практически одинаково за расчетный период при Ж0>Жиг и Ж0<Жпг (см. рис. 4). Несколько более значительное выравнивание концентрации пара в образце с влажностью ниже предела гигроскопичности связано с менее интенсивным производством пара в этом случае (кривые 1м, 2м, рис. 4).
0
0
Заключение
Диффузионный перенос пара в парогазовой среде -одной из фаз коллоидного капиллярно-пористого тела может оказывать значительное влияние на кинетику и динамику сушки. Численное моделирование процесса сушки в материалах такого типа может быть действенным инструментом в поиске эффективных режимов сушки, в том числе, с помощью регулирования условий внешнего массообмена. Режимы внутреннего диффузионного переноса пара и массообмена у внешних границ образца определяют диффузионные критерии Пекле и Нуссельта.
Литература
1. Дорняк О.Р., Шульман З.П. Математическое моделирова-
ние процесса сушки натуральной и уплотненной древесины // Изв. вузов. Сев.-Кавк. регион. Техн. науки. 2005. Спец. выпуск. Композиционные и порошковые материалы. - С. 133-142.
2. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. - М., 1984.
3. Шубин Г.С. Сушка и тепловая обработка древесины. - М.,
1990.
4. Крутов В.Н., Исаев С.И., Кожинов И.А. и др. Техническая термодинамика. - М., 1991.
23 мая 2006 г.
Воронежская государственная лесотехническая академия
УДК 628.517.2
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССА ПОГЛОЩЕНИЯ ЗВУКОВОЙ ВОЛНЫ ЗЕЛЕНЫМИ НАСАЖДЕНИЯМИ
© 2006 г. Н.Н. Панюшкин, П.И. Попиков, А.Н. Панюшкин, В.П. Попиков
При прохождении звуковой волны через листву деревьев часть ее энергии поглощается, т.е. преобразуется в механическую энергию колеблющихся листьев.
Каждый лист можно представить себе как пружинный маятник с массой т0 и коэффициентом жесткости к (рис. 1).
1
d 2 х at2
. dx
F,
+ 28 — + ro x = ——cos rot = A0 cos rot, (1)
at
ro 0
где 5 - коэффициент затухания, 8 =
2m 0
Ес = бппЯу = гу ; г = 6ппЯ , где п - коэффициент динамической вязкости воздуха: при Т = 293 К п = 1,72-10 -5 Пуаз; при Т = 300 К п = 1,84-10 -5Пуаз ; Я - эффективный радиус листа, определяется по формуле Я = . £ , где £ - площадь листа.
Решение дифференциального уравнения (1) для установившихся колебаний будет иметь вид [1]
X =
Ä0
■ cos ф ,
:2 го2
Дсо2-го2) + 482 где го0 - циклическая частота свободных незатухающих колебаний листа: го 0 = ; к, т0 - коэффици-
Рис. 1. Расчетная схема колебаний листа: С - центр тяжести листа; О - точка крепления листа; АХ - смещение центра тяжести листа от положения равновесия
Под действием звуковых волн лист будет совершать вынужденные колебания, которые можно описать дифференциальным уравнением второго порядка:
ент жёсткости и масса одного листа; А0, го - амплитуда и циклическая частота колебаний волны.
Интенсивность плоской и сферической волны звука определяется как I = и < Ж >= 2 рдго2А 0 , где
< Ж > - среднее значение плотности энергии за период звуковой волны; и - скорость звука в воздухе [2]:
~ЯГ
и = < Y
где R = 8,31-
Дж
M
- молярная газовая постоянная;
; m0 - масса
одного листа, кг; г - коэффициент сопротивления движению листа, теоретически может быть определен по формуле Стокса:
мольхК
у - постоянная адиабаты, у = 1,4 ; ц - молярная масса, для воздуха ц = 0,029 ^ ; Т - термодинамиче-
моль
ская температура.
r
Для падающей волны, возбуждаемой источником звука, т.е. движущимся по улицам автотранспортом 10 = и < Ж0 > , где < Ж0 > - энергия колебаний звука на пороге слышимости, < Ж0 >=< Ж > + < А Ж > , < Ж >=< Ж0 > - < АЖ >, < Ж > - прошедшая энергия; < АЖ > - поглощенная энергия (часть энергии, которая преобразуется в энергию механических колебаний листа): <АЖ >= пЖл, где п - концентрация листьев (число листьев в 1 м3) кроны дерева; Жл -энергия механических колебаний одного листа:
Wл =
m о ю 2(X о*)2
I и(< w0 > - < Aw >
1 2 л 2 1 2 1 \2
) ^ рю Ао - 4 р лю (Xо)
U<W0 >
1 2 .. 2 -рю2 А 2
= 1 -1 £i 2 р
( X * ^2
Ао
X о*
Ао
Ао =
А0 ^(ю 2-ю2)2 + ;8 2ю2 0 ^(ю 2-ю2)2 + ;8 2ю2
Тогда окончательно имеем формулу для определения коэффициента Ь:
< aw >= Nmm± ю 2(X *)2 =р л ю 2( X *)2
где р л
V
плотность листвы; р л =
N V'-
здесь N - коли-
чество листьев в кроне дерева; V - объем кроны дерева.
Для характеристики ослабления звуковых волн будем исследовать величину уровня интенсивности звука (коэффициент ослабления) Ь, измеряемый в белах (или единицами в десять раз меньшими - децибелами):
L = lg
( 1 Л I о
где I о - интенсивность звука на пороге слышимости, равная 1о-12 Вт/м2.
■3
L = lg
1 -
1 Р.
1
2 р (ю 2-ю 2)2 + ф8 2ю 2
Минимальное значение коэффициента Lm
1 Р л 1
L min = lg
1-
2 р ф 2 ю2
Численный эксперимент проверяли на примере насаждений липы, как наиболее распространённой среди зелёных насаждений в городских условиях [3, 4]. Математическая модель решена программой Ыа^Саё при следующих исходных данных включающих параметры и физико-механические свойства листьев и крон деревьев.
- 6
L := 9о • 1о C := 17о S := о.оо785 m := бо • 1о fmax:= 4оо N := 5ооооо i := 1.. N nl := 5о95 Vk:= 34.4 pl := 148^ m r := 1.72^ 1Ö
3
fmin:= о.1 5
Sk := 42.5
delf :=
(fmax- fmin)
R := S
Ki := log
NC юо := —
1-
f. := fmin+ i • delf i
r:= 6• n-r| • R Pl
Sl := Sk nl
S:=-
2 • m
Sl = 8.342X 1о
ю; := 2 • п • f. 1 i
Г1 := 1.84^ 1о pv := 1.29
5
- 2 "I
2 • pv • 2 / \2 юо - (ю;) + 4 • S • (ю;)
Kmin := log
1-
юо = 1.683X Ш
Pl
3
2 • pv •
[• S2 •(ю1)2
К, "5 -1о
"1 -1о
"8
167о
168о
ю i
а)
= 2о9.333
169о
-о.1
1о
( (Ol
б)
2о
Рис. 2. Графики зависимостей коэффициентов К и К^ от частоты колебаний автотранспортного шума
2
r
о
о
9
На рис. 2а показан график зависимости ослабления от частоты звуковой волны. Из графика следует, что максимальное ослабление соответствует частоте собственных колебаний листа и находится в диапазоне (250...300) Гц.
На рис. 2б показана зависимость максимального ослабления Ктт от резонансной частоты. Из графика следует, что Ктт нелинейно убывает с ростом частоты.
Разработанная математическая модель процесса звукопоглощения листьями деревьев позволяет прогнозировать степень ослабления автотранспортного шума от геометрических и биофизических параметров крон лиственных деревьев, используемых для зеленых насаждений в селитебной зоне городов и поселков.
Одним из способов интенсификации технологических процессов прессования древесины является использование пульсирующих нагрузок. Установлено, что применение направленных вибраций специально подобранной частоты и амплитуды улучшает качество вырабатываемой продукции, способствует автоматизации и механизации трудоемких процессов, повышает производительность труда и улучшает технику безопасности [1].
Принципиальная гидравлическая схема пресса для прессования древесины с пульсирующей нагрузкой представлена на рис. 1.
В качестве пульсатора можно применять вращающейся золотник, аксиально-поршневой насос с одним или двумя работающими плунжерами и другие устройства [2]. Расход жидкости гидропульсатора должен обеспечивать утечки и упругие деформации элементов гидропривода и прессуемой древесины. Гидропульсотор включается в работу после достижения гидроцилиндром рабочего усилия. После чего золотник распределителя, управляющего работой гидропульсатора, переводится в рабочее положение. Работа гидропульсатора осуществляется от гидростанции пресса.
Литература
1. Попиков В.П. Математическая модель процесса звукопоглощения шума зелеными насаждениями в селитебной зоне городов // Математическое моделирование, компьютерная оптимизация технологий, параметров оборудования и систем управления лесного комплекса: сб. науч. тр. / под ред. проф. В. С. Петровского; ВГЛТА. - Воронеж, 2005. - С. 10-14.
2. Трофимова Т.И. Курс физики: учеб. пособие для вузов. -М., 1999.
3. Подольский В.П. Воздействие транспортного шума, вибрации и электромагнитного излучения на окружающую среду в зоне влияния автодорог: учеб. пособие / ВГАСА. -Воронеж, 1996.
4. Репринцев Д.Д., Попиков В.П. Зеленые насаждения, как средство защиты от автомобильного шума / ВГЛТА -Воронеж, 2001. - 8 с. - Деп. В ВНИИТИ 21.05.01, № 1292 - В 2001.
24 мая 2006 г.
ПГ
Рис. 1. Гидравлическая схема пресса с гидропульсатором: ПГ - пульсатор гидравлический; Р1 и Р2 - гидрораспределители; КП1 и КП2 - гидроклапаны предохраните льные; КР - гидроклапан редукционный; РР - регулятор расхода; КО - клапан обратный; МН1 и МН2 - манометры;
Ц - гидроцилиндр
Математическая модель прессования древесины на гидравлическом прессе с гидропульсатором может быть представлена системой дифференциальных уравнений [3]:
Воронежская государственная лесотехническая академия
УДК 674.049.2
МАТЕМАТИЧЕСКАЯ МОДЕЛЬ РАБОЧЕГО ПРОЦЕССА ГИДРАВЛИЧЕСКОГО ПРЕССА С ГИДРОПУЛЬСАТОРОМ
© 2006 г. П.И. Попиков, Р.В. Юдин
dP
1010 (Р +1)
0,45
dt
= (1 - cos rot )sin ф-^F 2,5gP;
V Y
d 2 x
dx P0 nD ц nD ц .
m-+ u—+ cx =-—+-— P sin фt,
dt dt 4 4
где d - диаметр цилиндра; Rg - радиус окружности заделки поршневых шатунов в наклонном диске; ф -угол, образованный осями цилиндрового блока; ю -угловая скорость насоса; Р - давление рабочей жидкости, развиваемое пульсатором, Па; ц - коэффициент расхода для конусных клапанов; у - объёмная сила тяжести; g - ускорение силы тяжести; т - масса подвижных элементов пресса; с - жесткость упругой системы;
С помощью пакета программ «М&квтайса 4.0» были решены дифференциальные уравнения и получены графики, описывающие динамику работы гидравлического пресса с пульсатором (рис. 2, 3).
Р, МПа I
Рн
10 8 6 4 2 0
t, с
А, мм 1,5
1,0
0,5
-0,5 -1,0 -1,5
5
10
15
20
t, с
2 4 6 8
Рис. 2. График пульсации рабочего давления в зависимости от времени
На рис. 2 представлен график пульсации рабочего давления в зависимости от времени, из которого следует, что амплитуда колебаний находится в пределах 8.15 МПа, частота 5.6 колеб/с.
Рис. 3. График зависимости амплитуды пульсации от времени прессования
Из рис. 3 видно, что амплитуда колебаний возрастает по мере уплотнения древесины за 10.15 с примерно в два раза.
Математическая модель динамики гидропривода пресса с пульсирующей нагрузкой позволяет определить параметры гидропульсатора и динамические характеристики колебательной системы, установить оптимальные режимы технологического режима прессования.
Литература
1. Амалицкий В.В., Бондарь В.Г. и др. Теория и конструкция
деревообрабатывающих машин: учеб. пособие. - М., 1983.
2. Башта Т.М. Объемные насосы и гидравлические двигатели гидросистем. - М., 1974.
3. Попиков П.И., Юдин Р.В. Математическая модель процесса прессования древесины на гидравлическом прессе // Математическое моделирование, компьютерная оптимизация технологий, параметров оборудования и систем управления лесного комплекса: сб. науч. тр. / под ред. проф. В. С. Петровского / ВГЛТА. - Воронеж, 2005. -С. 26-30.
Воронежская государственная лесотехническая академия
1 июня 2006 г.
УДК 630.617
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПОЧВООБРАБАТЫВАЮЩЕГО АГРЕГАТА С РЕКУПЕРАТИВНЫМ ГИДРОПРИВОДОМ
© 2006 г. В.И. Посметьев, Е.А. Тарасов, Е.В. Снятков
При работе на лесных вырубках ходовая часть трактора и почвообрабатывающее орудие постоянно подвергаются знакопеременным нагрузкам вследствие наличия большого числа препятствий в виде пней,
корней, камней и неровностей поверхности. Поэтому, для повышения экономической эффективности, целесообразно оснащать их системами рекуперации энергии. Для изучения возможности оснащения рекупера-
тивным гидроприводом почвообрабатывающего агрегата в составе серийных трактора ДТ-75М и культиватора КЛБ-1,7 (рис. 1) разработана математическая модель, описывающая динамическое поведение агрегата в процессе работы на вырубке.
В [1] предложено использовать рекуперативную систему, встроенную в гидросистему трактора. Элементами, непосредственно воспринимающими внешние усилия, являются гидроцилиндры, устанавливаемые между балансирами кареток трактора, гидроцилиндр навесной системы и гидроцилиндр предохранительного механизма культиватора. Для оптимизации параметров элементов рекуперативной гидросистемы
была разработана имитационная компьютерная модель.
В рамках модели почвообрабатывающий агрегат рассматривается как плоский механизм, состоящий из семи твердых тел (рис. 2), для каждого из которых известны координаты центра тяжести (х,, у,), масса т, и центральный момент инерции Ji. Исходя из конструкции агрегата, определен набор точек, в которых тела контактируют друг с другом с помощью шарниров (12-22, 24-32, 13-52, 42-54, 63-71), невесомых нерастяжимых тяг (01-11, 14-62, 15-61) и пружин (23-33, 43-53, 64-72).
Рис. 1. Общий вид и исследуемые элементы почвообрабатывающего агрегата: 1 - трактор; 2 - лесной дисковый культиватор с гидравлическим предохранителем; 3 - звенья механизма навески трактора; 4 - опорный каток; 5 и 6 - внешний и внутренний балансиры каретки; 7 и 8 - оси качания внутреннего и внешнего балансиров; 9 - пружина; 10 - гидроцилиндр навесного механизма; 11 - автоматическая сцепка; 12 - рама культиватора; 13 - дисковая батарея; 14 - поворотная стойка дисковой батареи; 15 - рамка дисковой батареи; 16 - гидроцилиндр предохранителя культиватора
23
22
30
у.
61
62
х6
К7
40
21 Т2
31 Т3
41 Т4
51 Т5
Рис. 2. Представление почвообрабатывающего агрегата в виде совокупности твердых тел
в рамках предлагаемой модели
В основе математической модели лежит система дифференциальных уравнений Лагранжа I рода с неопределенными множителями в виде
p ЭФ mх ¡о + £X =QX1;
Эх,
p ЭФ
тг У г о + XX = Qyi;
.=1 oy 1о
p ЭФ
JФго + . ^ = Qфг, .=1 дФ 10
где Qxi, Qyi - декартовы составляющие равнодействующих сил, приложенных к ■-му телу; Qфi■ - соответствующий момент; - неопределенные множители Лагранжа (5 = 1, 2,..., р); Ф5 - функции связей.
Для составления системы уравнений используется метод [2], основанный на конечно-элементном подходе, согласно которому общая система уравнений составляется из уравнений-шаблонов для соответствующих связей (шарнир, тяга, пружина). Полученная система имеет, укрупненно, следующий вид:
M T" Г х > Г Q-1
T' O_ X V У = 1 и J
(1)
где М - квадратная матрица масс и моментов инерции размерностью 3п х 3п (п = 7 - число тел); Т - прямоугольная матрица размерности 3п х 3пх П - суммарное число степеней свободы, которые «отнимают» у системы все наложенные связи); Т' - транспонированная матрица Т размерности 3п^х3п; О - нулевая матрица размерности 3п^х3п^; Qx - вектор размерности 3 п, где каждый элемент представляет собой сумму всех соответствующих коэффициентов правой части исходных уравнений-шаблонов, выбранных и вычисленных на основании описания массива связей, а также независимые возмущений; и - вектор размерности п%, образующийся из совокупности коэффициентов и■ уравнений-шаблонов.
Для вычисления сил, действующих на ходовую часть трактора (в точках 21, 31, 41, 51) и дисковый рабочий орган со стороны почвы и препятствий, используется разработанная ранее модель [3]. В систему
(1) добавляются также уравнения, описывающие основные элементы гидропривода. Для численного интегрирования полученной системы дифференциальных уравнений используется модифицированный метод Эйлера. Для проведения компьютерных экспериментов используется специально составленная в среде Borland Delphi 7 программа.
В процессе компьютерного эксперимента моделируется движение почвообрабатывающего агрегата по контрольному участку вырубки длиной 1 км. Вероятность встретить препятствие (пень или корень) подчиняется равномерному закону. Для расчета данной вероятности при определенной плотности пней используются результаты работы [4]. В процессе движения модельного агрегата фиксируется зависимость вертикальных перемещений центра тяжести трактора от времени y1(t) и затем строится соответствующая амплитудно-частотная характеристика (АЧХ) A(m). Для оптимизации параметров элементов гидропривода компьютерный эксперимент проводится многократно. При этом критерием оптимизации является некоторый функционал АЧХ F [А(ю)].
Литература
1. Посметьев В.И., Тарасов Е.А., Кухарев В.С. Перспективные рекуперативные системы для гидроприводов лесных почвообрабатывающих агрегатов // Наука и образование на службе лесного комплекса: Сб. матер. МНПК Воронеж, 2005. - С. 132-136.
2. Расчет и проектирование строительных и дорожных машин на ЭВМ / под ред. Е.Ю. Малиновского. - М., 1980.
3. Посметьев В.И., Посметьев В.В. Моделирование взаимодействия дискового рабочего органа лесного почвообрабатывающего орудия с почвой и препятствиями // Математическое моделирование, компьютерная оптимизация технологий, параметров оборудования и систем управления лесного комплекса: межвуз. сб. науч. тр. Воронеж, 2000. - С. 39-44.
4. Посметьев В.И., Пухов Е.В., Посметьев В.В. Обоснование на основе результатов компьютерного моделирования выбора трактора по критерию его проходимости // Лес. Наука. Молодежь ВГЛТА 2004: сб. науч. тр. Воронеж, 2004. - С. 160-165.
Воронежская государственная лесотехническая академия
6 июня 2006 г.
УДК 621.22
ДИНАМИКА ГИДРОПРИВОДА ТРЕЛЁВОЧНОГО ЗАХВАТА ДЛЯ БЕСЧОКЕРНОЙ ТРЕЛЁВКИ СОРТИМЕНТОВ
© 2006 г. П.И. Попиков, С.И. Федяинов
В настоящее время основным направлением в техническом прогрессе лесного комплекса является разработка и внедрение прогрессивных технических и технологических решений, обеспечивающих макси-
мальное повышение производительности труда и создание условий для рационального использования всей биомассы дерева при минимальном воздействии машин на лесную среду и оператора.
Нами была составлена математическая модель [1], описывающая рабочие процессы в гидросистеме сельскохозяйственного трактора МТЗ-82.1, во время трелёвки пачки деревьев клещевым захватом УТБ-0,8 (рис. 1):
d! dt:
h hA
hx
dp = 1 ~dt = Kr
- P
nD 2 h1 + h 5
h1 h 4
Робж h 2 +"
G б n
- P
тр
nD dx q Н n Н----a y P
Н л ж у
где t - время, с; п - коэффициент длины дерева; Р -давление рабочей жидкости в напорной магистрали, Па; Р2 - давление рабочей жидкости в сливной магистрали, Па; Gб - вес бревна, Н; Ртр - сила трения, Н; Кр - коэффициент податливости упругих элементов гидропривода, м3/Па; дн - рабочий объём насоса, м3/об; пН - частота вращения вала насоса, 1/с; ау -коэффициент утечек рабочей жидкости, м3/с-Па.
Р, МПа 6 4 2
Рис. 2. Зависимость давления в гидросистеме горизонтального перемещения захвата от времени: тд = 400 кг, ау = 10-10 м3/(с-Па), Кр = 2 • 10-11 м3/Па
На рис. 3 показана зависимость коэффициента пульсации давления от времени.
k
-пульс
Л
[ay= 1о-
м3/сПа
о
2
4
6
8 Kp, 1о-11 м3/Па
Рис. 3. Зависимость коэффициента пульсаций давления в гидросистеме горизонтального перемещения захвата от коэффициента податливости п = 0,66; Б = 0,075 м; Рсл = 2 • 105 Па; qН = 32 • 10-6 м3; пН = 20 1/с; Нх= 0,25 м; Н2 = 0,93 м; Н2 = 1,02 м; Н3 = 0,45 м; Н4 = 0,39 м; Н5 = 0,45 м; Робж = 1000 Н; тд = 400 кг
Введение в гидросистему демпфирующих элементов способствует увеличению коэффициента податливости Кр и, следовательно, уменьшению коэффициента пульсации &пульс,- При Кр > 1012 динамического повышения давления в гидросистеме практически не наблюдается, что обеспечивает повышение надёжности лесных машин.
4
2
Рис. 1. Расчётная схема трелёвочного захвата
Решая данную модель, получили графические зависимости основных параметров: давления Р, горизонтальной скорости захвата Vx и горизонтальной координаты Х от времени при следующих значениях: п = 0,66; Б = 0,075 м; Рсл = 2 • 105 Па; qН= 32 • 10-6 м3; пН = 20 1/с; Н1= 0,25 м; Н2 = 0,93 м; Н2 = 1,02 м; Н3 = =0,45 м; Н4 = 0,39 м; Н5 = 0,45 м; Робж = 1000 Н.
На рис. 2 представлен график зависимости давления в гидроцилиндре захвата, из которого видно, что максимальное давление в начале движения составляет 4 МПа, а время переходного процесса - 0,6 с.
Литература
1. Попиков П.И., Гончаров П.Э., Федяинов С.И. Математическая модель рабочего процесса гидросистемы колёсного тягача с трелёвочным захватом // Математическое моделирование, компьютерная оптимизация технологий, параметров оборудования и систем управления. - Воронеж, 2005. - Вып. - 10. С. 76-82.
1 июня 2006 г.
Воронежская государственная лесотехническая академия
