2017 Математика и механика № 47
УДК 532.517.4
Б01 10.17223/19988621/47/9
Р.Р. Турубаев, А.В. Шваб
ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ АЭРОДИНАМИКИ ЗАКРУЧЕННОГО ПОТОКА В ВИХРЕВОЙ КАМЕРЕ КОМБИНИРОВАННОГО ПНЕВМАТИЧЕСКОГО АППАРАТА
Представлено численное моделирование аэродинамики закрученного течения в вихревой камере. Особенностью такой камеры является наличие вращающихся лопаток, расположенных в верхней выходной для несущей среды части аппарата. Дополнительно исследуется динамика несущей среды в этой же камере с расположенным в средней части аппарата вращающимся диском для увеличения окружной скорости несущей среды. Численные исследования показали возможность получения равномерного по высоте профиля радиальной скорости на входе в ротор, что является необходимым условием для эффективной работы сепаратора. Достоверность результатов численных исследований устанавливается тестовыми исследованиями, сравнением с известной аналитической зависимостью, а также сравнением полученных решений в переменных «скорость - давление» и «функция тока - вихрь».
Ключевые слова: вихревая камера, скорость, численное моделирование, давление, аэродинамика, вихрь, закрученное течение, функция тока, угловая точка.
В наши дни потребность в получении тонкодисперсных порошков с заданным гранулометрическим составом значительно возросла. Порошковые материалы наиболее широко используются в активно развивающейся аддитивной технологии, в порошковой металлургии, химической, атомной и других отраслях промышленности. Одним из наиболее эффективных и экологически безопасных способов получения тонкодисперсных порошков с заданными размерами частиц являются пневматические методы переработки сыпучей среды. К такому способу получения тонкодисперсных порошков относится разработанный в Томском университете комбинированный пневматический аппарат, который сдержит два совмещенных процесса - измельчение и сепарацию. Экспериментальные исследования показывают, что определяющим фактором, влияющим на процесс сепарации в таком аппарате, является аэродинамика закрученного потока несущей среды. Таким образом, настоящая работа посвящена численному моделированию поля скорости газовой среды в сепарационной зоне комбинированного центробежного аппарата, а также исследованию влияния режимных и геометрических параметров на аэродинамику закрученного течения, что позволит использовать разработанную математическую модель при оптимизации технологии получения порошков заданного фракционного состава.
Физико-математическая постановка задачи
Сепарационная камера комбинированного пневматического аппарата (рис. 1) представляет собой цилиндрическую область, в нижней части которой располагается труба с подводящим газом и твердыми частицами. В верхней части камеры находится ротор, состоящий из вращающих с постоянной угловой скоростью лопаток. В средней части сепарационной камеры может располагаться вращающийся диск для увеличения окружной составляющей вектора скорости двухфазного потока.
Несущая среда поступает в вихревую камеру через сечение А-А' (рис. 1) и под действием перепада давления обтекает вращающийся диск и вместе с мелкой фракцией проходит через систему вращающихся лопаток и соответственно выходит через сечение Б-Б'. Часть крупной фракции отсеивается на начальном этапе, где крупные частицы, в отличие от мелких, сталкиваются с диском и падают вниз, для дальнейшего измельчения. Остальная доля крупной фракции отсеивается системой вращающихся лопаток Б-В-В'-Б', в силу того, что центробежная сила значительно превышает силу аэродинамического сопротивления частиц [1]. Также рассматривался вариант геометрии с отсутствием центрального дискового элемента.
z /
Г1 Г2
Г1 Г2
Рис. 1. Схема расчетной области Fig. 1. Scheme of the computational domain
z
В силу большого количества лопаток в окружном направлении можно считать, что в зазорах между вращающимися лопатками для газовой фазы имеет место постоянство угловой скорости вращения. Это обстоятельство позволяет считать, что в камере осуществляется осесимметричное течение относительно оси вращения, т.е. отсутствуют изменения скорости относительно окружной координаты, что позволяет рассматривать закрученное течение вязкого газа в осесимметричной постановке задачи.
Движение закрученного потока описывается уравнениями Навье - Стокса и уравнением неразрывности. В силу того что область цилиндрическая, удобнее
всего их рассматривать в цилиндрической системе координат. В безразмерной форме и с учетом осевой симметрии эти уравнения примут следующий вид:
диг + д(игиг) + д(1иА + -иФ =-др + дт дг дг г дг Яе
Гд\ дг2
д2и
+ 1 д дг2 г дг
(1)
дмф + д(МфМг) + д(МфМг ) + 2мфмг = _1_
дт дг дг г Яе
2
дг2
д V
дг2
Ф +1 дмФ мф
г дг
диг + д(игиг ) + д(игиг ) + «г^т. =+
дт дг дг г дг Яе
2
дг2
ЗУ
1 д
дг
2 г дг
(2)
(3)
диг диг иг —- +—- + — = 0 ,
дг дг г
(4)
где иг, иф, и2 - составляющие вектора скорости, Яе = и0г2/V - число Рейнольдса, V - кинематическая вязкость, и0 - скорость потока во входном сечении, г2 - радиус входного сечения. В силу небольших скоростей плотность газа считается постоянной.
Методы численного решения
Полученная система уравнений решается методом физического расщепления полей скорости и давления. После расщепления уравнений (1) - (4) задача разбивается на 2 этапа:
и — и , . * 1 2 * ^ п
-+ (и-У)и =—V2и — Урп ;
дт Яе
ип+1 — и*
= —V(5p).
(5)
(6)
Умножая на градиент последнее уравнение и учитывая, что на (п+1)-м слое выполняется уравнение неразрывности, в результате получаем уравнение Пуассона для поправки к давлению:
V■u ч
-= V2 (5р).
дт
Исследуемая задача является стационарной, поэтому уравнение Пуассона можно представить в виде нестационарного уравнения:
Ш ^2 (5р) —
дт
(7)
После того как стали известны поправка к давлению и значение вектора скорости на промежуточном слое, рассчитываются скорость и давление на (п+1)-м слое:
ип+1 = и* — дт^(5р); (8)
рп+1 = рп +5р.
(9)
Для получения единственного решения ставятся следующие граничные условия. На всех границах для поправки к давлению используется условие Неймана
[4]. На входе безразмерная осевая компонента скорости принимается равной единице, для радиальной составляющей вектора скорости ставится условие Неймана и окружная компонента скорости равна нулю. На оси выполняются условия симметрии. На выходе для всех искомых величин задается условие установления решения. На стенках выполняются условия прилипания, в результате чего образуются следующие критерии:
Ц =-
Q2 =-
ю2 r2
и0 и0
Здесь Ю], ю2 - угловые скорости диска и стенки соответственно. Эти критерии являются обратными числами Россби.
Для дальнейшего решения используется обобщенный неявный метод переменных направлений [2, 3]. Суть метода состоит в расщеплении шага по времени с целью построения многомерной неявной схемы, в которой требуется обращение только к трехдиагональной матрице. Он обладает вторым порядком точности по времени. Для наглядности, продемонстрируем описанный метод на примере уравнения Пуассона:
59 id 2 9 д 2 9^
= a —+ -
dt 1дх dy
д0* a д2 д9* Гд V n Гд 29^|
--- дх2 a 1&2 у + a 1ду2 У
At 2
a д2 д9** д9*
At
dy 2
At
еи+1 = е" +дб .
Конвективные и диффузионные члены аппроксимировались с помощью экспоненциальной схемы, которая снимает ограничение на сеточное число Рей-нольдса и имеет второй порядок точности относительно пространственных координат [4]. Решение задачи проводится на разнесенной сетке (рис. 2).
Рис. 2. Расчетная сетка и контрольные объемы для иф, p (слева), ur (посередине) и uz (справа), где х - ur, о - uz и ■ - иф, p Fig. 2. Computational mesh and control volumes for u,p and p (on the left), ur (in the middle), and uz (on the right) where х, о, and ■ indicate ur, uz, and u^ p, respectively
Для верификации полученного решения, система (1) - (4) также решалась в переменных «вихрь - функция тока». Вихрь и составляющие вектора скорости вводятся следующим образом:
и=1 дг; (1о)
г дг
^ =-^; (11) г дг
ю=^. (12)
дг дг
Подставив (10) и (11) в (12), получим уравнение Пуассона для определения функции тока у:
д 2Ц1 д 2у 1
—— + —— = юг +--— . (13)
дг2 дг2 г дг
Введем итерационный параметр т:
~дт
( д \ д ч дг2 дг2 у
1 5ш
= -юг---—.
г дг
Перекрестным дифференцированием и последующим вычитанием уравнений (1) и (3) получим уравнение переноса вихря:
дю дигю дигю 2мф ди 1 (д2ю д2ю 1 дю ю 1 — +-Г— + —+ —^—^
дт дг дг г дг Яе
+-+
^ 2
дг2 дг2 г дг г2
Разностный аналог граничных условий для вихря на стенке принимает вид формулы Тома [5]:
ю„=+0(лг), гдп
где п - нормаль к поверхности стенки.
Так же как и в предыдущем случае, для получения единственного решения ставятся следующие граничные условия. На входе вихрь задается равным нулю, значение функции тока определяется посредством интегрирования из (11). На оси симметрии ставятся условия симметрии. На выходе ставится условия Неймана для всех переменных.
На первом этапе работы для расчета вихря в угловой точке использовался способ Ричардсона. Суть его заключается в том, что значения вихря в угловых точках зависят от направления расчета. Следовательно, в каждой угловой точке вихрю будет присвоено два различных значения, которые рассчитываются по формуле Тома. Однако, как показали результаты (рис. 3), вблизи угловой точки существует различие в решениях, полученных двумя различными способами.
В связи с этим были проанализированы несколько схем расчета вихря в угловых точках [2]. Сравнение этих методов приведено на рис. 4, откуда видно, что результаты, полученные при использовании способа Кавагути, лучше согласуются с результатами, полученными в естественных переменных. По этой причине, именно этот способ взят за основу для дальнейших расчетов.
Рис. 3. Сравнение двух методов решения в сечении Д-Д' (- - «вихрь - функция тока», • «скорость - давление») Fig. 3. Comparison of two solution methods in the Д-Д' section (- - and • indicate the cases with «vortex - stream function» and «velocity - pressure», respectively)
Рис. 4. Сравнение схем расчета вихря в угловой точке в сечении Д-Д' (- - способ Ричардсона, □ стенка, наклоненная под углом 45°, — способ Кавагути, • «скорость - давление»)
Fig. 4. Comparison of the computational schemes of the vortex at the corner point in the Д-Д' section (- -, □, —, and • indicate the Richardson method, the case with a wall inclined at an angle of 45°, Kawaguchi method, and the «velocity - pressure» case, respectively)
На рис. 5 - 8 показано влияние геометрических параметров на распределение функции тока и окружной составляющей. В частности, на рис. 5 видно, что поток прижимается к верхней стенке камеры, что, безусловно, препятствует более продуктивной работе сепаратора. Посредством простого изменения геометрии, например показанного на рис. 7, можно добиться равномерного распределения изолиний функции тока при тех же режимных параметрах. Также можно отметить, что картина изолиний угловой скорости стала более симметрична (рис. 8), чем была до изменения (рис. 6).
На рис. 9 и 10 присутствуют циркуляционные зоны. Но если в первом случае, при вращающемся диске, они не препятствуют равномерности распределения функции тока, то во втором случае, когда диск не вращается, а стенка под выходным сечением закручивается, они существенно поджимают изолинии. Вследствие этого можно заключить, что изменяя режимные параметры, также достаточно просто можно добиться более равномерного распределения.
На рис. 11 и 12 видно, что при отсутствии диска, за счет изменения режимных параметров, так же как и в предыдущих случаях, можно добиться равномерного распределения функции тока (рис. 12). Либо можно получить такую картину течения, при которой в выходной области будет образовываться завихренная область, препятствующая равномерному распределению поля вектора скорости на выходе из рабочей области (рис. 11).
На рис. 13-14 можно заметить, что в случае, когда скорость вращения диска меньше скорости вращения ротора, картина течения в непосредственной близости с системой вращающихся лопаток изменяется существенно.
3.5
3 2.5 2 1.5
1
0.5
Рис. 5. Распределение функции тока в рассматриваемой области. Re = 10, П1 = 3, H2 = 5 Fig. 5. Stream function distribution in the region considered at Re = 10, П1 = 3, and H2 = 5
3.5
3 2.5 2 1.5
1
0.5
Рис. 7. Распределение функции тока в рассматриваемой области. Re = 10, Hi = 3, H2 = 5 Fig. 7. Stream function distribution in the region considered at Re = 10, H1 = 3, and H2 = 5
Рис. 6. Изолинии угловой скорости в рассматриваемой области. Re = 10, H1 = 3, H2 = 5 Fig. 6. Isolines of the angular velocity in the region considered at Re = 10, H1 = 3, and H2 = 5
0.5 1 1.5 2 2.5
Рис. 8. Изолинии угловой скорости в рассматриваемой области. Re = 10, H1 = 3, H2 = 5 Fig. 8. Isolines of the angular velocity in region considered at Re = 10, H1 = 3, and H2 = 5
3.5 3 2.5 2 1.5 1
0.5
0.5
1
1.5 2
2.5
Рис. 9. Распределение функции тока в рассматриваемой области. Re = 10, П1 = 10, H2 = 0 Fig. 9. Stream function distribution in the region considered at Re = 10, Hi = 10, and H2 = 0
3.5 3 2.5 2 1.5 1
0.5
0.5
1
1.5 2
2.5
Рис. 11. Распределение функции тока в рассматриваемой области. Re = 10, H2 = 10, H3 = 5 Fig. 11. Stream function distribution in the region considered at Re = 10, H2 = 10, and H3 = 5
3.5 3 2.5 2 1.5 1
0.5
0.5
1
1.5 2
2.5
Рис. 10. Распределение функции тока в рассматриваемой области. Re = 10, H1 = 0, H2 = 10 Fig. 10. Stream function distribution in the region considered at Re = 10, H1 = 0, and H2 = 10
3.5 3 2.5 2 1.5 1
0.5
0.5
1
1.5 2
2.5
Рис. 12. Распределение функции тока в рассматриваемой области. Re = 10, H2 = 1, H3 = 1 Fig. 12. Stream function distribution in the region considered at Re = 10, H2 = 1, and H3 = 1
3.5 3 2.5 2 1.5 1
0.5
0.5
1
1.5 2
2.5
Рис. 13. Распределение функции тока в рассматриваемой области. Re = 10, H1 = 8, H2 = 10 Fig. 13. Stream function distribution in the region considered at Re = 10, H1 = 8, and H2 = 10
3.5 3 2.5 2 1.5 1
0.5
0.5
1
1.5 2
2.5
Рис. 15. Распределение функции тока в рассматриваемой области. Re = 10, H1 = 2, H2 = 4 Fig. 15. Stream function distribution in the region considered at Re = 10, H1 = 2, and H2 = 4
3.5 3 2.5 2 1.5 1
0.5
0.5
1
1.5 2
2.5
Рис. 14. Распределение функции тока в рассматриваемой области. Re = 10, H1 = 2, H2 = 10 Fig. 14. Stream function distribution in the region considered at Re = 10, H1 = 2, and H2 = 10
3.5 3 2.5 2 1.5 1
0.5
0.5
1
1.5 2
2.5
Рис. 16. Распределение функции тока в рассматриваемой области. Re = 10, H1 = 4, H2 = 2 Fig. 16. Stream function distribution in the region considered at Re = 10, H1 = 4, and H2 = 2
В противоположном случае, распределение изолиний у выходной области меняется незначительно. Например, как видно из рис. 13, при значениях закрутки диска = 8 и ротора = 10 образуются множественные циркуляционные зоны. Тем не менее в обоих случаях можно подобрать величины закрутки таким образом, чтобы на выходе имелось равномерное распределение профиля радиальной скорости. Один из вариантов такого подбора представлен на рис. 15 и рис. 16.
Рис. 17. Сравнения двух методов решения в сеченииях Д-Д' (а) и Г-Г' (b) (— «вихрь - функция тока», • «скорость - давление») Fig. 17. Comparison of two solution methods in the Д-Д' (a) and Г-Г' (b) sections (— and • indicate the cases with «vortex - stream function»
and «velocity - pressure», respectively)
2.4 2.6 2.8 R
Рис. 18. Сравнение с известной аналитической формулой (— аналитическое решение,
• расчетные данные) Fig. 18. Comparison with the known analytic formula (— and • indicate the theoretical solution and computational data, respectively)
Рис. 19. Тестовое исследование
на сеточную сходимость. (□ 41х51; • 81x101; — 161x201) Fig. 19. Test study for a grid convergence (□, •, and — indicate the grid sizes 41x51, 81x101, and 161x201, respectively)
Достоверность полученных результатов определялась сравнением решений, полученных при помощи двух различных способов («функция тока - вихрь» и «скорость - давление») в разных сечениях (рис. 17) сравнением полученного численного решения с известной аналитической формулой (14) (рис. 18), а также тестовыми исследованиями на сеточную сходимость (рис. 19):
(r2 -r2)lnrf -(rK2 -r42)lnf
Uz = 2ucp-к-. (14)
rk - Г42 + (rK2 + Г) ln f
к
Выводы
Разработана математическая модель для численного расчета аэродинамики закрученного потока в циркуляционном аппарате. Получены распределения полей вектора скорости и давления. Исходя из численных расчетов, вихрь в угловой точке рекомендуется рассчитывать при помощи метода Кавагути. Из анализа результатов можно заключить, что более равномерного профиля скорости в области ротора можно достичь путем сокращения длины лопаток ротора (рис. 8). В случаях, когда угловая скорость вращения ротора больше угловой скорости вращения диска, картина течения вблизи выходной области существенно изменяется. В противоположном случае, она практически остается неизменной. Однако равномерного профиля радиальной составляющей вектора скорости можно достичь как в первом, так и во втором случае и, основываясь на тех же численных результатах, можно заключить, что, как правило, это достигается при относительно небольшой разности скоростей вращения ротора и диска. Полученная математическая модель может быть использована как для модификации существующих пневматических аппаратов, так и для создания новых.
ЛИТЕРАТУРА
1. Патент РФ № 2522674 Способ газовой центробежной классификации и измельчения порошков / Зятиков П.Н., Росляк А.Т., Шваб А.В., Демиденко А.А., Романдин В.И., Брендаков В.Н. Опубл., Б.И. № 20, 20.07.14.\
2. Роуч П. Вычислительная гидромеханика. М.: Мир, 1977. 618 с.
3. Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен: пер. с англ. С.В. Сенина, Е.Ю. Шальмана; под ред. Г.Л. Подвидза. М.: Мир, 1990. Т. 1. 384 с.
4. Патанкар С.В. Численные методы решения задач теплообмена и динамики жидкости: пер. с англ. под ред. В. Д. Виленского. М. : Энергоатомиздат, 1984. 149 с.
5. Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен: пер. с англ. С.В. Сенина, Е.Ю. Шальмана; под ред. Г.Л. Подвидза. М.: Мир, 1990. Т. 2. 726 с.
6. Лойцянский Л.Г. Механика жидкости и газа: учебник для вузов по специальности «Механика». М.: Наука, 1987. 840 с.
Статья поступила 21.10.2016 г.
Turubaev R.R., Shvab A.V. (2017) NUMERICAL STUDY OF SWIRLED FLOW AERODYNAMICS IN THE VORTEX CHAMBER OF THE COMBINED PNEUMATIC MACHINE. Tomsk State University Journal of Mathematics and Mechanics. 47. pp. 87-98
DOI 10.17223/19988621/47/9
This paper presents a numerical simulation of the aerodynamics of a swirled flow in a vortex chamber. A particular feature of this vortex chamber is the presence of rotating blades located at the top of the outlet part for the carrying agent. In addition, the dynamics of the carrying agent in the abovementioned vortex chamber with a rotating disk located in the middle part and designed to increase the circumferential velocity of the carrying agent has been investigated. The numerical investigation showed that it was possible to get a uniform profile along the height of the radial velocity at the rotor inlet. This fact is a necessary condition for the effective operation of the separator. The reliability of numerical simulation results has been confirmed by both the test study, the comparison with the well-known analytical dependence, and the comparison of the results obtained using «velocity-pressure» and «vortex-stream function» variables.
Keywords: vortex chamber, velocity, numerical simulation, pressure, aerodynamics, vortex, swirled flow, stream function, corner point.
TURUBAEV Roman Rinatovich (Student, Tomsk State University, Tomsk, Russian Federation) E-mail: [email protected]
SHVAB Aleksandr Veniaminovich (Doctor of Physics and Mathematics, Professor, Tomsk, Russian Federation) E-mail: [email protected]
REFERENCES
1. Zyatikov P.N., Roslyak A.T., Shvab A.V., Demidenko A.A., Romandin V.I., Brendakov V.N. Sposob gazovoy tsentrobezhnoy klassifikatsii i izmel'cheniya poroshkov [Gas Centrifugal Classification and Grinding of the Powders]. Patent of RU 2522674, 2014.
2. Rouch P. (1977) Vychislitel'naya gidromekhanika [Computational Fluid Dynamics]. Moscow: Mir.
3. Anderson D., Tannehill J., Pletcher R. (1990) Vychislitel'naya gidromekhanika i teploobmen [Computational Fluid Dynamics and Heat Transfer]. Vol. 1. Moscow: Mir.
4. Patankar S.V. (1948) Chislennye metody resheniya zadach teploobmena i dinamiki zhidkosti [Numerical Techniques for Solving Heat Transfer and Fluid Dynamics Problems]. Moscow: Energoatomizdat.
5. Anderson D., Tannehill J., Pletcher R. (1990) Vychislitel'naya gidromekhanika i teploobmen. Tom 2. [Computational Fluid Dynamics and Heat Transfer]. Vol. 2. Moscow: Mir.
6. Loytsyanskiy L.G. (1987) Mekhanika zhidkosti i gaza : uchebnik dlya vuzov po spetsial'nosti "Mekhanika" [Fluid and gas mechanics: tutorial for high schools specialty "Mechanics"]. Moscow: Nauka.