CC BY
101
геометрические характеристики форсунки Ь, 5, S¡l¡, Sm мы получаем объемный расход жидкости, время истечения, угол распыла и дисперсность получаемого аэрозоля.
Показано, что для обеспечения более высокой дисперсности аэрозоля необходимо обеспечить импульсный характер воздействия на жидкость;
предложен безразмерный критерий Ш, характеризующий определяющий механизм диспергирования. Расчеты нашли качественное экспериментальное подтверждение.
Работа выполнена при финансовой поддержке гранта Российского Фонда Фундаментальных Исследований (№ 12-08-90810-мол_рф_нр).
СПИСОК ЛИТЕРАТУРЫ
1. Кудряшова О.Б., Ворожцов Б.И. Математическая модель взрывной генерации жидкокапельных аэрозолей // Известия Томского политехнического университета. - 2011. - Т. 318. -№2. - С. 77-81.
2. Кудряшова О.Б., Ворожцов Б.И., Муравлев Е.В., Ишматов А.Н., Павленко А.А. Ударно-волновая генерация высокодисперсных жидкокапельных аэрозолей // Ползуновский вестник. - 2010. - №4-1. - С. 95-100.
3. Ворожцов Б.И., Кудряшова О.Б., Ишматов А.Н., Ахмадеев И.Р., Сакович ГВ. Взрывная генерация высокодисперсных жидкокапельных аэрозолей и их эволюция // ИФЖ. - 2010. -Т 83. - № 6. - С. 1084-2004.
4. Витман Л.А., Кацнельсон Б.Д., Палеев И.И. Распыливание жидкости форсунками / ред. С.С. Кутателадзе - М.; Л.: Гос-энергоиздат, 1962. - 264 с.
5. Кириллин В.А. и др. Техническая термодинамика. - М.: Энер-гоатомиздат, 1983. - 417 с.
6. Физическая энциклопедия / Д.М. Алексеев, А.М. Балдин, А.М. Бонч-Бруевич, А.С. Боровик-Романов, Б.К. Вайнштейн, С.В. Вонсовский, А.В. Гапонов-Грехов, С.С. Герштейн, И.И. Гуревич, А.А. Гусев, М.А. Ельяшевич, М.Е. Жаботин-ский, Д.Н. Зубарев, Б.Б. Кадомцев, И.С. Шапиро, Д.В. Шир-ков / под общ. ред. А.М. Прохорова. - М.: Советская энциклопедия, 1988-1998.
7. Кедринский В.К. Газодинамика взрыва: эксперимент и модели. - Новосибирск: Издательство СО РАН, 2000. - 435 с.
8. Блох А.Г., Кичкина Е.С. Распыливание жидкого топлива механическими форсунками центробежного типа // Сб. Вопросы аэродинамики и теплопередачи в котельнотопливных процессах. - М.: Госэнергоиздат, 1958. - С. 132-140.
9. Пажи Д.Г., Ламм Э.Л. Распыливающие устройства в химической промышленности. - М.: Химия, 1975. - 200 с.
Поступила 28.06.2012 г.
УДК 674.815-41
МОДЕЛИРОВАНИЕ ПОСТУПАТЕЛЬНО-ВРАЩАТЕЛЬНОГО ДВИЖЕНИЯ ДРЕВЕСНЫХ ЧАСТИЦ
ПРИ ОРИЕНТИРОВАНИИ
С.М. Плотников, М.С. Лурье
Сибирский государственный технологический университет, г. Красноярск E-mail: smplotnikov@rambler.ru
Представлены результаты моделирования ориентирования плоских древесных частиц, используемых при изготовлении плит из крупноразмерной стружки. Получены зависимости времени падения и скорости разворота частиц отих параметров и высоты падения, позволяющие минимизировать угол укладки частиц в стружечном ковре.
Ключевые слова:
Плоская частица, имитационная модель, время падения, импульс вращения, коэффициент аэродинамического сопротивления.
Key words:
Flat particle, simulation model, fall time, momentum of rotation, drag coefficient.
Современная технология изготовления плит из крупноразмерной ориентированной стружки (Oriented Strand Board - OSB) предусматривает их ориентирование при формировании стружечного ковра, который в дальнейшем подвергается горячему прессованию. Благодаря укладке частиц наружных слоев брикета в продольном, а внутренних - в поперечном направлении возможно значительное повышение прочности на изгиб изготовляемых плит. Причем чем меньше угол разброса частиц, тем выше прочность готовой плиты, которая приближается к прочности исходной древеси-
ны. Ориентирование частиц является одной из основных операций в технологическом процессе изготовления OSB, что требует ее подробного изучения с целью совершенствования ориентирующих устройств.
Зная динамику движения древесной частицы, в частности изменение скорости и времени падения и разворота от ориентирующего органа до укладки в стружечный ковер, можно так подобрать действующий на частицу вращающий импульс, что она будет уложена в ковер с наименьшим углом отклонения. Это позволит значительно улучшить ос-
новной механический показатель OSB - прочность на статический изгиб.
Рассмотрим поступательное и вращательное движения частицы раздельно.
Без учета вращательного движения на частицу, свободно падающую в воздушной среде, действуют сила тяжести mg и сила сопротивления воздуха F1.. Дифференциальное уравнение поступательного движения частицы имеет вид:
Л
dt2
т—2 = mg -FB,
(1)
d2 s V2
m^F = mg - CA'p T
(2)
Программа Сош8о1 МиШрЬу81е8 использует обобщенную версию уравнений Навье-Стокса, чтобы учесть переменную вязкость. Уравнения имеют вид:
ди
Р — -У[п(Уи + (V и)Т)] + р( иУ) и+ У р = ^
dt
и Vu = 0,
(4)
где т — масса частицы, кг; 5 - путь, проделанный частицей, м; g=9,81 м/с2 - ускорение свободного падения.
С учетом силы Fъ, которая зависит от скорости движения частицы относительно воздуха, уравнение (1) принимает вид:
где Сх — коэффициент аэродинамического (лобового) сопротивления частицы; ^ —площадь проекции поперечного сечения частицы на плоскость, перпендикулярную направлению ее движения, м2; рв=1,29 кг/м3 — плотность воздушной среды при нормальных условиях; V — скорость поступательного движения при падении частицы, м/с.
Масса частицы зависит от ее размеров и удельной плотности материала:
т=ВЬкрч, (3)
где В - ширина частицы; Ь - длина частицы; к — толщина частицы; рч — плотность материала частицы (400...700), кг/м3.
Площадь ^ при плоскопараллельном падении частицы определяется по формуле: Зё=ВЬ.
Для определения скорости V по уравнению (2) необходимо знать величину коэффициента аэродинамического сопротивления Сх. В [1] значения Сх при малых числах Рейндольса Re имеют большой разброс и недостоверны, т. к. получены на основании малого числа экспериментов. Коэффициент Сх нелинейно зависит от числа Re. Последнее, в свою очередь, зависит от скорости частицы V. Исключение из расчетов значений Сх возможно за счет моделирования аэродинамических процессов, возникающих при падении частицы в воздушной среде.
Расчет падения частицы с учетом ее аэродинамики заключался в определении зависимости силы торможения частицы Fъ от времени ее падения численным моделированием в программе Сош8о1 МиШрЬу81е8 3.5 [2]. Получаемые значения силы F в виде аппроксимирующего выражения вводились в 5-модель падения частицы, созданную в программе МаЛаЬ^тиИпк. При этом скорость и время падения частицы определялись с учетом силы аэродинамического торможения в воздухе. Таким образом, необходимость использования нелинейной зависимости коэффициента Сх от числа Рейнольдса отпадает.
где п - динамическая вязкость, Па-с; и - векторное поле скоростей; p - давление, Па; F - объёмная плотность силы, действующей на среду, в том числе и сила тяжести, Н/м3.
Первое уравнение - это уравнение переноса количества движения, а второе - уравнение неразрывности потока для несжимаемой жидкости.
Поскольку длина частицы в десятки раз больше ее поперечных размеров, то для сокращения машинного времени моделирования решалась плоская задача. При этом мы пренебрегаем аэродинамическим сопротивлением на концах частицы -это приводит к погрешности величиной не более 2...3 %, что вполне допустимо при таких расчетах.
Важной особенностью движения частиц в устройствах ориентирования является ламинарный характер обтекания частиц воздушной средой. Так, для характерных размеров частиц длиной до 150 мм, шириной до 30 мм, толщиной до 1 мм [3] при их свободном падении или вращении со скоростями вращательного движения в пределах от 1 до 8 рад/с число Рейнольдса, характеризующее обтекание частиц потоком, составляет от 10 до 800, что соответствует ламинарному режиму.
Ламинарный режим движения однородных сред и их взаимодействие с обтекаемыми телами хорошо подчиняется численным методам расчета, поэтому для определения коэффициента аэродинамического сопротивления частиц был применен численный эксперимент. Он заключается в следующих шагах:
• производится расчет стационарного процесса движения частицы в потоке воздуха при нормальных условиях методом численного моделирования с помощью программы Comsol Multiphysics 3.5. Для того чтобы при 2D-моделиро-вании учесть длину частицы, при постпроцес-сорной обработке результатов определялась погонная суммарная сила F/l, действующая на обтекаемое тело, в противоположном движению частицы направлении;
• по рассчитанному значению погонной силы определяется коэффициент аэродинамического сопротивления частицы;
• результаты, полученные для частиц близкой формы, усредняются;
• усредненные результаты аппроксимируются аналитическими выражениями, которые будут пригодны для дальнейших расчетов.
Сила взаимодействия потока с обтекаемым им телом описывается формулой [4]:
F V2
V=C Рв В—.
l
2
Тогда коэффициент аэродинамического сопротивления может быть найден:
C =-
(5)
0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11
Рис. 1. Поле скоростей через 0,2 с после начала падения
По результатам моделирования падения частицы получена зависимость силы сопротивления воздуха F от числа Рейнольдса, аппроксимируя которую с помощью пакета Curve Fitting программы MatLab, получим выражение (коэффициент корреляции ^=0,99991):
В качестве примера на рис. 1 показана картина обтекания частицы (без учета ее вращения) через 0,2 с после начала движения.
BL B0 L0
х(0,0007845 Яе2 + 0,0001813Яе + 4,05 -10-5),
где В0=0,002 м - средняя ширина частицы; Х0=0,1 м - средняя длина частицы.
Имитационная модель для данной задачи составлена в пакете Зішиїіпк. Выражение (4) в модели рассчитывает блок Fcn2. Схема модели представлена на рис. 2, некоторые результаты моделирования - на рис. 3.
В результате моделирования установлено, что время падения частицы Т0 зависит от ее толщины и плотности и не зависит от длины и ширины частицы, т. к. при увеличении этих величин во столько же раз возрастает площадь сечения частицы и, соответственно, растет сила аэродинамического торможения.
Рис. 3. Зависимость времени падения частицы отвысоты Но падения при 8=0,01 м; L=0,1 м; Ь=0,0005 м; рч=400 кг/м3:1 - для воздуха; 2 - для вакуума
Для оценки аэродинамического коэффициента к в выражении (6) необходим анализ вращательного движения частицы.
Пусть плоская частица длиной Ь, шириной В и толщиной к вращается вокруг оси ОО1 с угловой скоростью о (рис. 4).
На рисунке: Я=Ь/2; - элементарная сила со-
противления; о - скорость вращения частицы; г -текущая координата рассматриваемого участка частицы.
Сила аэродинамического сопротивления / по длине частицы будет изменяться, поскольку разные ее участки будут двигаться с разной линейной скоростью относительно воздуха. Рассмотрим бесконечно малый участок частицы длиной йг на расстоянии г от оси вращения. Элементарная сила сопротивления действующая на участок частицы длиной йг и шириной В:
V2
ОР = Ср БСг-^,
где К - линейная скорость перемещения краев частицы, м/с.
Зная время падения частицы Т0, можно определить скорость ее вращения [рад/с], необходимую для минимизации угла укладки частицы в стружечный ковер:
, 2п п +а„„
о= к-----—, (6)
* 0
где а, - угол схода частицы с направляющих, зависящий от расстояния между направляющими и длины частицы [рад]; п=0 или 1 - количество оборотов, задаваемое частице; к - коэффициент аэродинамического торможения, учитывающий замедление частицы при вращении.
Число п определяется производительностью ориентирующего устройства. Чем оперативнее требуется ориентировать частицы, тем больший импульс вращения следует им придать, тем больше должна быть начальная скорость о0. При этом частица может совершить дополнительный оборот. В последнем случае возрастает также угол разброса частиц в ковре. На основании (6) можно определить оптимальную скорость направляющих элементов в системе ориентации древесных частиц [5].
Выразив линейную скорость через угловую VK=m, получим:
dF = CxpB Bdr
2
Элементарный тормозящий момент, считая и силу, воздействующую на вторую половину частицы: dM=2dFi=Cxp11B®2tidr. (7)
Для определения суммарного тормозящего момента необходимо проинтегрировать выражение (7) по длине частицы от оси вращения до ее конца:
^ R
M = JR dM = JСхрвB m2r3dr =
0
R R
= рвВт2 JCxr3dr = K JCxr3dr, (8)
0 0
где Х=рвВт2.
Коэффициент аэродинамического сопротивления Сх в (8) зависит от скорости движения частицы относительно воздуха, которая, в свою очередь, определяется текущей координатой рассматриваемого участка и угловой скоростью т. Для общности последующих выводов будем считать, что коэффициент Сх имеет следующий вид:
C = A +-^т = A + -^ = А + ^ г. (9)
VV vmr r
Частицы шириной 4, 15 и 25 мм и толщиной 0,2, 0,5 и 1,0 мм располагались при моделировании широкой стороной вдоль и поперек потока. Пример полученных зависимостей погонной силы от скорости вращения частиц приведен на рис. 5.
После расчета по выражению (5), усреднения результатов и аппроксимации кривых средствами программы MatLab получены выражения, позволяющие для каждой из групп частиц с точностью до 10 % определять коэффициент Сх.
Для частиц, расположенных шириной В вдоль и поперек потока соответственно:
Cx =-3,449 +
407,14 л/Re ,
Cx =1,3293 +
2,5227
VRe '
практического расчета движения древесных частиц в системах ориентирования он неудобен, т. к. скрытно содержит размеры частиц и их скорость, которые приходится аналитически интегрировать в процессе расчета. В нашем частном случае для последующего анализа удобнее выразить Сх как функцию скорости движения частиц относительно воздуха. Для частиц, расположенных шириной В вдоль и поперек потока соответственно:
Cx = -3,449+
11,218
Подставив (9) в (8), получаем:
R R С
8), п
M = K f Ar 3dr + К
О 0 V W r
R TS
= KA Jr3dr +
Cx =1,238 +
0,15
VF •
r dr =
KCRr3
J dr.
0 V® 0
Выражение (10) содержит два интеграла:
Rr 3dr =R и V dr =2 R.
J 4 JJT 7
0 0 > r '
Подставляя (11) в (10) получим:
,, KAR4 2KCR34R M =--------+-
(10)
(11)
4
7yfä
Из (10) получим выражение тормозящего момента:
^ рв BAR4 2 M =—----------m +
2рвBCR3VR m2
Введем обозначения:
а =-
рв BAR4
2рв BCR 3VR
(12)
4 7
Тогда окончательно получим зависимость аэродинамического момента торможения от размеров частицы и ее угловой скорости:
,2
M = am + b
(13)
Выражение (9) получено для аргумента, которым является число Рейнольдса. Хотя это общепринятый критерий движения потоков среды, для
Для определения закона движения частицы запишем второй закон Ньютона для вращающейся частицы:
В=4 мм В=15 мм В=25 мм
Скорость, м/с
Рис. 5. Погонная сила, действующая на частицы с шириной В поперек потока
г da
J-------= -M,
dt
(14)
где J - момент инерции частицы относительно оси ОО1, кгм2.
Для пластины массой m и длиной Ь момент инерции, как известно, составляет [6]:
2
Т 1 т-2 1 Л! m2 mR
J =— mL =—m(2R) =------------
12 12 3
Подставив (15) в (14), получим:
da
dt
mR
a®2 +
a >—
a
2 Л
(15)
(16)
Интегрируем обе части полученного дифференциального уравнения (16):
a =--
mR2
П 2 , ®
II aa +b—j=
1 v®
dt.
(17)
Масса частицы m не входит в коэффициенты уравнения (12), поэтому, согласно (17), она определяет скорость вращения частицы в обратной пропорции.
Уравнение (17) нелинейно, поэтому для его интегрирования применим численные методы, а именно имитационное моделирование с помощью пакета ЗтиКпк программы МаНаЬ.
Общий вид модели для плоских частиц представлен на рис. 6.
Блок подготовки данных вычисляет значение коэффициентов уравнения (17), блок решения дифференциальных уравнений решает (17), интегрируя его правую часть. Блок обработки данных вычисляет угловую скорость частицы во времени. При этом вычисляются два значения углового отклонения частицы: с учетом сил аэродинамического торможения и идеальное - без учета данных сил. Данные значения индицируются на дисплеях Ап-gle_degrees (в градусах) и А^1в_% (в процентах
от идеального значения). Осциллограф Scope показывает процесс изменения угловой скорости при движении частицы во времени.
Моделирование проводилось для плотности воздуха 1,204 кг/м3, его кинематической вязкости 1,51-Ю“5 м2/с и плотности частицы 400 кг/м3.
Толщина В частицы, встречающей поток, не влияет на замедление частицы при вращении. Это следует из выражений (12) и (13). Коэффициенты a и b в выражении (12) прямо пропорциональны ширине частицы В. Но, масса плоской частицы, определяемая по формуле (3), в формуле момента инерции (15) также пропорциональна ширине частицы В.
Учитывая, что R=L/2 и проведя сокращения, получим:
12
Ис„
16
a2 +р
СЬ24Ь a2 ^ 7л/2 4a
dt.
В то же время, ширина частицы Н оказывает влияние на скорость замедления. Чем она больше, тем меньше замедляет вращение частица. Это можно объяснить тем, что увеличение Н приводит к увеличению момента инерции частицы, мало влияя на аэродинамическое сопротивление. Однако следует учитывать, что вследствие увеличения площади трения частицы о воздух с увеличением ее ширины Н аэродинамическая сила торможения также несколько возрастает. Это приводит к увеличению скорости торможения частицы.
Полученные в результате численного эксперимента поля скоростей потоков при обтекании частицы (без учета ее падения) представлены на рис. 7, графики изменения угловой скорости частицы - на рис. 8.
В результате моделирования поступательного движения плоской частицы получено значение времени ее падения в зависимости от ее размеров и высоты падения. Данное время зависит от ее толщины и плотности частицы и не зависит от ее длины и ширины. Определены также отклонения вре-
Блок подготовки данных
Ц>
ш> 2 £ [А>
о
Блок решения диф. уравнения
Блок обработки данных
a
С
4и
1/s
т
1/s a + U-57,296 Aa 9.607
X U-100 Aa% -10.48
■
Scope
Angle_ degrees
Angle_%
Рис. 6. Имитационная модель вращения плоской частицы
R
2
U
b
Рис. 7. Поле скоростей при движении частицы шириной 20 мм, толщиной 1 мм: а) плоская сторона вдоль потока; б) плоская сторона поперек потока
L=75 мм L=100 мм L=150 мм
Время, с
Рис. 8. Зависимость ®=Ц\) для частицы с параметрами: Ь=0,5 мм, 8=20 мм, движущейся узкой гранью навстречу потоку
мени падения от среднего значения при сочетании параметров частицы случайным образом.
При вращении частицы действующая на нее погонная сила аэродинамического сопротивления возрастает нелинейно с увеличением размеров частицы. Угловая скорость частиц относительно начальной скорости снижается по закону, близкому к линейному, причем за время ориентирования для частиц длиной 150 мм угловая скорость снижается на 20 %, для частиц длиной 75 мм - на12 %. Данное снижение линейно возрастает с уменьшением плотности частиц. Таким образом, для всего диапазона диаметров частиц коэффициент аэродинами-
ческого торможения k в выражении (6) принимает значения от1,12 до 1,20 и зависит от размеров частиц, их плотности, начальной скорости вращения и времени разворота, которое пропорционально высоте ориентирующего устройства.
Точность расчетов удалось повысить за счет того, что коэффициент аэродинамического сопротивления в обеих моделях определялся в численном эксперименте.
Использование результатов исследования позволит уменьшить угол разброса частиц в формируемом стружечном ковре и тем самым повысить прочность на изгиб готовых плит ОББ.
СПИСОК ЛИТЕРАТУРЫ
1. Идельчик И.Е. Справочник погидравлическим сопротивлениям. - М.: Машиностроение, 1992. - 672 с.
2. Шмелев В.Е. Femlab 2.3. Руководство пользователя. - М.: Диалог-МИФИ, 1999. - 442 с.
3. Thoemen H., Irle M., Sernek M. Wood-Based Panels: An Introduction for Specialists. - London: Brunei University Press, 2010. -283 р.
4. Лойцянский Л.Г Механика жидкости игаза. 7-е изд. - М.: Дрофа, 2003. - 840 с.
5. Плотников С.М., Пантелеев В.И. Система оптимального ориентирования стружки в производстве древесных плит // Известия Томского политехнического университета. - 2008. -Т.313. - №2. - С. 56-58.
6. Фаворин М.В. Моменты инерции тел. 2-е изд., перераб. идо-полн. - М.: Машиностроение, 1977. - 514 с.
