CC BY
33
УДК 621.521
Ю. А. Паранин, А. В. Бурмистров, О. Ю. Паранина
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТЕПЛОВЫХ ПОЛЕЙ
РАБОЧИХ ЭЛЕМЕНТОВ СПИРАЛЬНЫХ МАШИН СУХОГО СЖАТИЯ
Ключевые слова: спиральный компрессор, спиральный вакуумный насос, теплообмен, тепловые поля, математическое моделирование, стационарная задача, нестационарная задача.
Описаны методики расчета тепловых полей рабочих элементов спиральных машин для решения стационарной и нестационарной задач, в том числе методы решения, алгоритм расчета граничных данных. Представлены результаты численного эксперимента, хорошо совпадающие с измеренными температурами.
Key words: scroll compressor, volute vacuum pump, heat exchange, thermal fields, mathematical simulation, stationary problem, time-
dependent problem.
Methodology for the analysis of thermal fields of the working elements of scroll machines used for solving the stationary and time-dependent problems has been described, including the solving methods and the computation algorithm for calculation of boundary data. Results of numerical experiment showing good agreement with the measured temperatures have been presented.
Введение
При математическом моделировании рабочих процессов спиральных машин определяющее значение имеет учет теплообмена в рабочих полостях и изменение монтажных зазоров в результате тепловых деформаций рабочих элементов (подвижной и неподвижной спиралей), в основном, за счет нагрева горячим сжимаемым газом. Так, например, не учет теплообмена и изменения зазоров от тепловых деформаций рабочих элементов воздушного спирального компрессора сухого сжатия при расчете внутреннего адиабатного КПД может привести к расхождению с экспериментальными данными на 12-32% [1]. Для расчета теплообмена необходимо знать тепловые поля поверхностей спиралей, образующих рабочие полости (вогнутая профилированная поверхность ребра неподвижной спирали, выпуклая профилированная поверхность ребра подвижной спирали, поверхность основания неподвижной спирали, поверхность основания подвижной спирали), а для расчета тепловых деформаций -тепловые поля всех поверхностей рабочих элементов. Спираль представляет собой сложное геометрическое тело, состоящее из основания и тонкого профилированного ребра, поэтому экспериментальное определение температурных полей рабочих элементов, в необходимом объеме, не всегда возможно. Для подвижной спирали такая задача практически неосуществима. В данной статье представлены методики расчета тепловых полей для эвольвентных спиралей, основанные на методах математического моделирования и предназначенные для решения стационарной и нестационарной задач. Спиральные машины (компрессоры и вакуумные насосы) в основном находят применение в технологических процессах, требующих их продолжительной работы. Однако в некоторых отраслях промышленные спиральные машины работают кратковременно, например, воздушные спиральные компрессоры сухого сжатия для стоматологии. Аналогичный режим работы имеет
место при откачке вакуумной камеры от атмосферного давления до предельного остаточного с помощью спирального вакуумного насоса.
Во время продолжительной работы эксплуатация машины проходит на установившемся режиме в прогретом состоянии (стационарная постановка задачи определения тепловых полей рабочих элементов). При кратковременной работе -на неустановившемся режиме, во время ее прогрева (нестационарная постановка задачи определения тепловых полей рабочих элементов).
Описание метода решения стационарной задачи
Разработанная методика предназначена для расчета стационарных трехмерных температурных полей подвижной и неподвижной спиралей, соответствующих средней по циклу температуре сжимаемого газа.
Стационарное распределение температуры в спирали описывается трехмерным уравнением теплопроводности (уравнением Лапласа):
АTi(х) = 0, х еО.;, i = 1,2. (1)
Здесь аi — трехмерная область, занимаемая спиралью, х = (х1, х2 х3) — точка трехмерного пространства.
A T (x) =
д 2T (x) + д 2T (x) + д 2T (x) dx2 dx2 dx 32
трехмерный
оператор Лапласа. Значение i = 1 соответствует неподвижной спирали, i = 2 — подвижной спирали, Т — температура спирали.
Граничные условия. Граница Г области ^ разбивается на две части Г = Г1 + Г2, Г1 — та часть границы, на которой задается температура (см. ниже описание граничных условий для неподвижной спирали). На Г1 ставится граничное условие Дирихле:
Т(х) = Т, (х), х еГ,, (2)
где Т, — заданная функция.
На Г2 ставится граничное условие третьего рода, описывающие теплообмен с внешней средой, имеющей заданную температуру, и
противоположной спиралью:
к ^Т^ + * х)(Т (х) - 7х (х)Т (х) - (3)
-у2( х)Гз-. (х)) = д( х), х еЦ где / = 1,2 , 7 (х) = а (х) + в( х),
(х) = а( х )/(а( х) + в( х)), у2( х) = в( х )/(а( х) + в( х)),
к - коэффициент теплопроводности материала спирали, а - коэффициент теплообмена с внешней средой (воздухом, охлаждающей водой, валом, корпусом), Т0 - температура внешней среды, в = X/ Ь - коэффициент теплообмена с
противоположной спиралью, X - коэффициент теплопроводности воздуха; Ь - расстояние между соответствующими участками подвижной и неподвижной спиралями; д - плотность теплового потока (за счет тепловыделения подшипников и от трения торцевой уплотнительной ленты).
Уравнение (1) с граничными условиями (2), (3) аппроксимировалось системой линейных алгебраических уравнений по методу конечных элементов [2]. Предварительно на подвижной и неподвижной спиралях строилась трехмерная призматическая сетка с треугольными основаниями. Сетка адаптировалась к поверхностям ребра спирали и основания (рис. 1).
60
х (мм)
Рис. 1 - Общий вид трехмерной сетки на рабочем элементе
При аппроксимации уравнения
использовались изопараметрические элементы с узлами, расположенными в вершинах призмы. При вычислении локальных матриц жесткости использовались формулы численного
интегрирования, позволяющие точно находить элементы этих матриц. Полученная в результате аппроксимации система уравнений решалась при помощи итерационного метода.
На первой итерации в принималось равным нулю. Температурные поля подвижной и неподвижной спиралей определяются при этом независимо. На всех последующих итерациях в
граничные условия для определения температуры спирали входит температура на поверхности противоположной спирали.
Каждый шаг итерационного метода требует решения системы линейных уравнений, аппроксимирующей уравнение Пуассона со смешанными граничными условиями.
Использовалась стандартная программа,
реализующая метод типа Гаусса для разреженных матриц.
Описанная методика реализована в виде комплекса программ, представляющего оконное приложение, и позволяющего варьировать исходные данные и параметры численного метода, отображать результаты расчетов в числовом и графическом виде.
Описание метода решения нестационарной задачи
В общем случае уравнения
теплопроводности рабочих элементов должны быть трехмерными. Однако математическое описание полей температур сложных геометрических тел в функции времени чрезвычайно сложно, поэтому приняты следующие допущения. Профилированные ребра спирали имеют большую длину, для возможности образования несколько камер сжатия. Вследствие этого, по длине ребра температура рабочего газа, с которым контактирует спираль, существенно меняется, а, следовательно, меняется температура спирали. Ближе к центру спираль имеет более высокую температуру, в то же время, высота ребра существенно меньше ее длины, поэтому температура по высоте ребра спирали изменяется незначительно. Таким образом, для моделирования тепловых полей ребра спирали, можно ограничиться одной пространственной переменной, которая отсчитывается вдоль направления закрутки спирали.
Основание спирали, на котором размещено профилированное ребро, также соприкасается с неоднородно нагретым газом, в результате температура основания максимальна в центре и понижается к периферии. В силу небольшой толщины основания, по сравнению с его диаметральным размером, изменения температуры по толщине не должны быть большими. Поэтому для расчета тепловых полей в основании достаточно учесть изменение температуры только по двум пространственным координатам, пренебрегая изменением температуры по толщине основания.
Вышеуказанные допущения справедливы для спиральных машин, у которых подвижная и неподвижная спирали при их эксплуатации находятся в схожих тепловых условиях.
Тепловой поток между профилированным ребром и основанием рассматривался как тепловой поток, ограниченный площадью их
соприкосновения и обусловленный разностью их температур. Как подвижная, так и неподвижная спирали разбиты на два конструктивных элемента -профилированное ребро и основание. Поэтому система уравнений содержит подсистемы нестационарных одномерных уравнений
теплообмена для профилированных ребер и нестационарных двумерных уравнений теплообмена для оснований.
Выделим в спиралях элементы ребра (рис.2а) и основания (рис.2б). Считая а — коэффициент температуропроводности, с — теплоемкость и р — плотность постоянными, напишем для них уравнения теплового баланса.
На рис.2 , qр, , ^р, ^^^ — пл°тн°сть
теплового потока и площадь элемента ребра в направлениях г , р , р. Здесь ,— высота и толщина ребра; р , р — полярные координаты.
дТ
= а Ц- - / Т + / Т )а1 + / (Т )а2,
б
Рис. 2 - К выводу дифференциального уравнения теплопроводности: а - для элемента ребра, б - для элемента основания
Согласно рис.2 изменение температуры ребра йТз в объеме йУз за время йт :
срйТ^йУ, = (+ йар + dQр)йт , (4)
где Qz, Qр, (2р — тепловые потоки в направлениях ^ р,р.
Рассчитав притоки теплоты через поверхности элемента в направлениях г, р,р, и учитывая принятые допущения, и выполнив необходимые преобразования уравнения (4) получим дифференциальное уравнение
теплопроводности для элемента — ребро спирали:
дт
д1
(5)
где / (Т1В = а
Т - Т
К (к8-дв)
1(Т )а1 ==а (Т§1 -Т ) ,
рсбз
рсб3
1 (Т )а2 = а (Т2 -Т ) ,
&
Яр15 = 1-—; 5р25 = 1+^; бз =-£■ ■
2р
2р
К
а„а2 — коэффициент теплоотдачи между газовым потоком и ребром; т§1, Т§2 — температура газового
потока; I — координата, отсчитываемая по длине ребра.
Уравнение представляет собой выражение теплопроводности с источниковыми членами ^Т8), которые описывают теплообмен между ребром и основанием а также между газовым потоком и ребром.
На рис. 2б qz, qr, qр, qSв, ^,2 3Вг, 3Вр, Ззв — плотность теплового потока и площадь элемента основания в направлениях г, г, р и в местах соприкосновения с ребром; qfr — плотность теплового потока от трения уплотнительной ленты; КВ — толщина основания; р, г — полярные координаты.
Согласно рис.2б изменение температуры основания йТв в объеме йУв за время йт :
срйТвйУв = (dQz + dQг + dQр)йт ,
(6)
где Qz, Qr, Q9 — тепловые потоки в направлениях
г, г, р .
Выполнив преобразования уравнения (6), получим дифференциальное уравнение
теплопроводности для основания спирали:
^ = а V 2Тв + 1 (Тв), дт
™ тт 1 д2 1 д д2
где V =--+--+-,
г2 др2 г дг дг2
1 (Тв ) = 1 (Тв \в + 1(Тв - 1(Тв - 1(Тв и — — функция источника.
(7)
1 (Тв )в = 2аТ^Гв
дЗз,
Кз -Кв дрг дг ■ К
1Т), = р
рс
ав
в
(Тв - Т),
1Т Ц а (Тв -Ъ
а
* (Тв 2).
ва Нврсу'а '*идргдг Нврс™ '*"дрдг ав/ ,авв1,авв2 — коэффициент теплоотдачи между
охлаждающей средой и основанием и между газовым потоком и основанием; ^, Т. 2 —
температура охлаждающей среды и газового потока.
Функция источника учитывает теплообмен между спиралью и основанием, от трения торцевой уплотнительной ленты, между охлаждающей средой и основанием, а также между газовым потоком и основанием.
Функции источника ^Т8) и ^Тв) могут учитывать теплообмен и с другими элементами, в зависимости от конструкции спиральных машин.
В результате дискретизации уравнения теплопроводности для ребра (5) была получена неявная разностная схема II порядка, построенная на шеститочечном шаблоне (рис. 3 а).
б
Рис. 3 - Шаблон конечно-разностной схемы: а -для ребра, б - для основания
Для основания (7) - неявная разностная схема II порядка, построенная на десятиточечном шаблоне (рис. 3б).
На рис.3 к-номер временного шага; 1, > номера пространственных шагов. Полученные конечно-разностные уравнения для ребра решались методом прогонки, для основания методом переменных направлений [3]. Описанная методика реализована в виде программы, построенной по модульному принципу и содержит следующие основные модули: головной модуль, модуль расчета коэффициентов в методе прогонки, модуль формирования граничных данных, модули прямой прогонки уравнений теплопроводности для ребер и оснований неподвижной и подвижной спиралей, модуль обратной прогонки, модуль формирования массива выходной информации и записи в файл.
Расчет граничных данных
Основной сложностью при задании граничных данных на поверхностях рабочего элемента, для вышеуказанных методов, является поверхность со стороны профилированного ребра (лицевая поверхность). На других поверхностях задание граничных данных не вызывает
затруднений. Далее приводится алгоритм расчета граничных данных для лицевой поверхности спирали.
Для задания средних по циклу значений температуры То лицевая поверхность спирали разбивается на шесть зон (рис. 4).
Первая зона образуется при движении полости вдоль внутренней поверхности ребра спирали, имеет ширину, равную толщине ребра 8 , ее границы задаются эвольвентами со следующими параметрами (здесь и далее углы перечисляются в следующем порядке: угол поворота системы координат, начальный угол закрутки, конечный угол закрутки, радиус начальной окружности всюду равен г0).
Рис. 4 - К расчету средней температуры газа на поверхности спирали: 1-3 - зоны, образующиеся при движении полостей вдоль спирали, 4 - ребро спирали, 5 - зона нагнетания, 6- зона всасывания
Внутренняя граница:
901 = 9о - 2S1 - S2 '
9rn =9о + 2S1 + s2 -ж, 9ki =9k + 2Si + ô2 - 2n, ô2 = 2n — 3S1,
Si = S/ r о
Внешняя граница:
9о n1 = Pо, 9nn1 =9n + п, 9kn1 = 9k ■ Вторая зона образуется при движении полости вдоль внешней поверхности ребра спирали, имеет ширину, равную S , ее границы задаются эвольвентами со следующими параметрами. Внутренняя граница:
902 =9о —ô 9n2 =9n +ô
9k 2 = 9k +ô-n Внешняя граница: 9о n 2 = 9о — S 9nn 2 =9n + 2Ô1,
9n 2 =9k + 2S1 —n
Третья зона есть пересечение следов от движения полостей вдоль внутренней поверхности ребра спирали и ее внешней поверхности. Она имеет ширину t - 3 S , где t - шаг ребра спирали, ее
границы задаются эвольвентами со следующими параметрами. Внутренняя граница:
Фоз = Фо - 28
Фп3 = Фп + 281 - я,
Фк 3 = Фк + 281 -п
Внешняя граница:
Фоп3 = Фо - 281 -82>
ФппЗ = Фп + 281 + 82
Фкп3 = Фк + 281 + 82
Четвертая зона - проекция профилированного ребра на основание, имеет ширину 8, границы задаются эвольвентами со следующими параметрами. Внутренняя граница: Ф0 = 0 - угол поворота системы координат, фп -начальный угол закрутки, Фк - конечный угол закрутки (задаваемые параметры).
Внешняя граница:
Фоп =Фо
Фш =Фп + 8
Фкп =Фк +81
Пятая зона - зона нагнетания - ограничена отверстием нагнетания и соответствующими участками границ зон 1 - 4. Шестая зона - зона всасывания - ограничена окружностью с центром в начале координат, радиусом, равным радиусу основания спирали и соответствующими участками границ зон 1 - 4.
Во всех указанных областях средние по циклу температура и давление сжимаемого газа определяются по методике расчета рабочего процесса. Для воздушного спирального компрессора сухого сжатия используется методика, изложенная в работах [1, 4, 5].
В зоне нагнетания давление и температура полагается равным давлению и температуре нагнетания соответственно, в зоне всасывания -давлению всасывания и температуре всасывания.
В первой зоне температура и давления полагаются зависящим лишь от угла закрутки ребра спирали, т. е. остаются постоянным по ширине зоны. При этом текущему углу закрутки приписывается объем полости, полученный по программе расчета рабочего процесса, значение температуры и давления определяются затем по PV и TV-диаграммам.
Аналогично вычисляются температура и давление во второй зоне. В третьей и четвертой зонах температура и давление также зависят лишь от угла закрутки. Они полагается равными средним значениям в соответствующих точках, примыкающих зон.
На боковых поверхностях ребра спирали температура и давление полагаются постоянными по высоте и равными температуре и давлению в соответствующей зоне.
В результате выполнения описанного алгоритма давление и температура оказываются определенными во всех точках лицевой поверхности спирали.
Коэффициент теплообмена а в зонах 1 - 4, 6 вычисляется по критериальному уравнению, например, расчет а для воздушного спирального компрессора сухого сжатия проводится по уравнению, полученному в работе[1]. В зоне нагнетания величина а полагается равной ее среднему значению на прилегающих частях границ областей 1-4. Отметим, что для расчета а используется средняя ширина зоны всасывания, которая предполагается задаваемым
конструктивным параметром.
Тепловой поток д отличен от нуля лишь в зоне трения, т. е в области, ометаемой за цикл поверхностью уплотнительной ленты на лицевой поверхности основания спирали (рис. 5).
80
60
40 20 0 -20 -40 -60
-60 -40 -20 0 20 40 60 80
Рис. 5 - К определению области трения: 1 -область трения уплотнительной ленты; 2 - ребро спирали
Ширина области трения равна г -1.5 8 . Границы задаются эвольвентами со следующими параметрами.
Внутренняя граница:
Ф о гг = Ф о -1.25 81, 9шг = Фг о +1.25 81 -п/2, Фк,г = Фк +1.2581 -п Внешняя граница: Фопгг =Фо + °-258, Фппгг =Ф1 о - о.258 + 3п/2, Фк„ = Фк - °-2581 + п где ф1 о - угол начала уплотнительной ленты.
Результаты численного эксперимента
Отладка программных модулей, описанных методов, проводилась на решении тестовых задач [6]. Сравнение теоретических и экспериментальных температур рабочих элементов позволили сделать вывод о хорошей качественной и количественной сходимости результатов. На рис.6 представлены некоторые результаты расчета и эксперимента стационарной температуры внутренней поверхности ребра неподвижной спирали, на середине его высоты, по углу закрутки эвольвенты, а так же аппроксимирующие зависимости температур первого порядка для воздушного спирального компрессора сухого сжатия, работающего на
режиме с охлаждением (п- отношение давлений; п - число оборотов вала; УВ - расход воды, охлаждающей с тыльной стороны неподвижную спираль). В данной машине использовались чугунные спирали. На рисунке сплошные линии соответствуют расчетным данным, штриховые -экспериментальным.
Т,К
425
405
385
365
345
325
О'У' t / /X
¿1 У / г / У
/ ' V Л VI
у\ I J
/ / /
' X ч / /
180
360
540
720
900 1080
ф.град
Рис. 6 - Зависимость распределения температуры внутренней поверхности ребра неподвижной спирали, на середине его высоты, от угла закрутки эвольвенты на режиме работы компрессора: п=7, п=3000 об/мин, УВ=4 л/мин
Сопоставление результатов расчета с экспериментальными данными на участке ф от 90° до 900° показало расхождение по абсолютной величине от 12 до 6 К или в относительном выражении от 14,5 до 3,5% [1]. Полученные результаты позволили использовать разработанные
методики при математическом моделировании рабочих процессов спиральных машин.
Статья подготовлена в ФГБОУ ВПО «КНИТУ» при финансовой поддержке проекта «Создание высокотехнологичного производства безмасляных спиральных вакуумных насосов для индустрии наносистем и наноматериалов» открытого публичного конкурса по отбору организаций на право получения субсидий на реализацию комплексных проектов по созданию высокотехнологичного производства согласно постановления Правительства Российской Федерации от 9 апреля 2010 года N 218 «О мерах государственной поддержки развития кооперации российских высших учебных заведений и организаций, реализующих комплексные проекты по созданию высокотехнологичного производства».
Литература
1. Ю.А. Паранин. Дисс.канд.техн.наук, Казанский национальный технологический ун-т, Казань, 2011. - 254 с.
2. Ф. Сьярле, Метод конечных элементов для эллиптических задач. Мир, Москва, 1980.512 с.
3. А.А. Самарский, Теория разностных схем. Наука, Москва, 1977. 656 с.
4. Ю.А. Паранин, Р.Р. Якупов, А.В. Бурмистров, Вестник Казанского технологического университета, 16,19, 267-271(2013)
5. Ю.А. Паранин, Р.Р. Якупов, А.В. Бурмистров, Вестник Казанского технологического университета, 17, 1, 248252 (2014).
6. Ю.А. Паранин, В.Б. Явкин, В тр. XIV международной научн. техн. конф. по компрессорной технике. Т.1, Казань, 2007. С.172-182.
© Ю. А. Паранин - к.т.н., доц. каф. компрессорных машин и установок КНИТУ, yurij.paranin@yandex.ru; А. В. Бурмистров -д.т.н., проф. каф. вакуумной техники электрофизических установок КНИТУ, burm@kstu.ru; О. Ю. Паранина - студ. КНИТУ.
