CC BY
50
А.Н.Поляков
ИСПОЛЬЗОВАНИЕ ПРИЗМАТИЧЕСКИХ КОНЕЧНЫХ ЭЛЕМЕНТОВ В ТЕПЛОВОМ МОДЕЛИРОВАНИИ СТАНКОВ
В работе приведены конечно-элементные матрицы теплопроводности, теплоемкости и вектора тепловой нагрузки для пяти- и шестигранных призматических конечных элементов. Рассмотрены два пути построения тепловой модели станка: численное интегрирование уравнения нестационарной теплопроводности и модальный подход. По итогам исследований сформулированы рекомендации по совершенствованию теплового моделирования станков.
Большую роль в формировании теплового состояния станка играют тепловые процессы в его базовых деталях. При этом наиболее точно эти процессы находят отражение в объемной тепловой модели. В условиях поиска максимальной инвариантности и автоматизации моделирования, реализация численных и численно-аналитических моделей является предпочтительной. В матричной форме уравнение теплопроводности имеет вид:
[С ]{Т} + [ К ]{Т} = } (1)
где [С],[К] - матрицы теплоемкости и теплопроводности размера пхп, соответственно;
{¥} - вектор тепловой нагрузки;
{Т },{Т}" п - мерные векторы узловых температур и их производных.
Во многих авторитетных автоматизированных системах численного моделирования (ANSYS, COSMOS,NASTRAN, DesignWorks, WinMachine) в качестве базового численного метода моделирования используется метод конечных элементов (МКЭ). Несмотря на многообразие имеющихся типовых конечных элементов, в литературе по МКЭ из объемных конечных элементов подробно рассмотрен лишь тетраэдр. Однако, при представлении больших базовых деталях станка достаточно удобно представлять деталь шести- и пятигранными призмами. Удобство использования этих конечных элементов также объясняется более простыми алгоритмами разбиения области на конечные элементы, во многом схожими с алгоритмами разбиения для плоской задачи.
Учитывая подробность изложения в отечественной и зарубежной литературе основных
уравнений МКЭ и порядка получения основных конечно-элементных матриц для типовых конечных элементов , в данной работе эти соотношения не приводятся.
Шестигранная призма. Рассмотрим построение уравнений для частного случая - прямой призмы с габаритами 2Ь х 2а х 2с, в основании которой задан прямоугольник. Принят следующий порядок нумерации: г,у,к,I - нижнее основание; т,п,о,р - верхнее основание призмы (рис.1 .а.).
В соответствии с теорией МКЭ, построение конечно-элементных матриц для типового конечного элемента начинается с выбора соответствующих функций формы. Для составления функций формы используются изопараметри-ческие координаты, центр которых совпадает с расположением центра тяжести призмы. Функции формы имеют вид:
М, = ^(1 -Х)(1 -п)(1 -С);
О
Му = 8(1 -Х)(1 -я)(1 + с);
Мк = 8(1 + Х)(1 - Л)(1 + С);
О
М1 = 8(1 + Х)(1 -П)(1 - С); (2)
О
Мт = *(1 -Х)(1 + П)(1 -с);
О
Мп = 8(1 -Х)(1 + П)(1 + С);
О
Мо = 8(1 + Х)(1 + П)(1 + С);
О
Мр = 8(1 + Х)(1 + л)(1 -С).
О
Выполнив предусмотренные теорией МКЭ процедуры были получены матрицы теплопро-
водности [K], теплоемкости [С ], а также вектор тепловой нагрузки ^}.
Как известно, матрица теплопроводности [К ] представляется двумя составляющими [ к 1 ] и [ К ]а. Первая составляющая [ К х ] обусловлена теплопередачей в станке только путем теплопроводности. Вторая составляющая [ К ]а связана с наличием конвективного теплообмена на открытых поверхностях станка.
Для составляющей [К х ] диагональные элементы принимают вид:
2Кха - с 2КуЬ - с 2К.а - Ь
ь _ х . ^ ^
к\^ -
к% = А,
1, ор 1,тп
Ч,к/
'■ К,у
К „а • с К уЬ • с 2 • К.а • Ь
к1,пр = к1,то = к1,у/ = к1,ік
9 • Ь 9 • а 9 • с Кха • с КуЬ • с К.а • Ь 9 • Ь 18 • а 9 • с
к = к,
1,по 1 ,тр
: к 1,ук = к1,й
к1 ,/р = к1,ко = к1,уп = к1 ,іт =
к1, /о = к1, кр = к1, ут = к1,і
2Кха • с КуЬ • с Кга • Ь
9 • Ь 9 • а 9 • с
Кха • с 2КуЬ • с К .а • Ь,
1,1п
1,кт
1,ур 1,іо
. х 9 • Ь ^ + 9 • а 9 • с
= Кха •с КуЬ • с К.а •Ь
‘ 18 • Ь 9 • а 1 •9 с
= Кха • с КуЬ • с К.а •Ь
18 • Ь 18 • а 18 • с
= Кха •с КуЬ • с К.а • Ь
9
к
а ,ік
' к а, у/
аЬс
9
[С ] =
срК
27
8 4 2 4 4 2 1 2
4 8 4 2 2 4 2 1
2 4 8 4 1 2 4 2
4 2 4 8 2 1 2 4
4 2 1 2 8 4 2 4
2 4 2 1 4 8 4 2
1 2 4 2 2 4 8 4
2 1 2 4 4 2 4 8
(4)
9 - Ь 9 - а 9 - с
s = г, у, к, I, т, п, о, р.
Все внедиагональные элементы матрицы [К 1 ] сводятся к следующим комбинациям:
к1,1т к1, кп к1, }о к1, гр п ж. л ю
у е 9 • Ь 9 • а 18 • с
Для второй составляющей матрицы теплопроводности [К]а для теплоотдающей грани г - у - к -1, компоненты матрицы принимают вид:
к к к к = 4 - аЬс а
к а,гг к а, уу к а,кк к а,Я 9 а;
к к к к = 2 - аЬс а
к а,гу = к а,И = к а, ук = к а,к1 = ^ а;
Аналогичные соотношения получаются для других граней с соответствующей перестановкой номеров узловых точек их описывающих.
Матрица теплоемкости в этом случае приобретает вид:
Ненулевые компоненты векторов тепловой нагрузки, обусловленные тепловым потоком от источников тепла {/ц} и конвективным теплообменом |/а} имеют вид: /ц, ^ = ц-К;
/а, ss = Т¥а-К; s = г, у, к, I, т, п, о, р, здесь ц - плотность теплового потока, проходящего через теплопринимающую поверхность деталей станка; а - коэффициент теплоотдачи, количественно характеризующий конвективный теплообмен деталей станка с окружающей средой; Т¥ - температура окружающей среды; у - объем конечного элемента.
Пятигранная призма. Геометрическое представление пятигранника высотой Н приведено на рис.1.б. Узлы г,у,к и I,т,п описывают прямоугольные треугольники, соответствующие треугольным граням пятигранника в плоскости иг.
Функции формы для пятигранника можно представить в виде:
[(1 - w) - ь ,Х = г, у, к
I w - Ьх ,Х = I, т, п 1
где ьх = 2^ (а| + ЬХ х + сХ У), здесь а^,Ь^,с^ - коэффициенты, зависящие от координат узловых точек треугольной грани;
а = хуУк- хуУк; Ьг = Уу- Ук; с = хк- ху (остальные коэффициенты получаются круговой заменой индексов);
А - площадь треугольной грани призмы.
Также как и для шестигранной призмы выпишем основные соотношения, формирующие конечно-элементные матрицы теплопроводности [ К ], теплоемкости [С ], и вектора тепловой нагрузки ^}.
Все компоненты составляющей матрицы теплопроводности [К ]х сводятся к четырем соотношениям и принимают вид:
N =■
= (КХЬ„2 + Кус2,)и + КА
ксс =
+—-—; s = г, /, к, I, т, п (6) 12 А 6Н
= Н(КхЬ*Ьг + Кус,сг ) , К-А
к *г =
12А
или s,г = I,т, п; s ^ г
+ Ж; s,г = ,’Л к;
к*г =
(KxЬs2 + Кс2)Н К-А
24А 6Н ’
s = г, у,к; г = I,т,п; (*,г) е ,
где б*г - условное обозначение ребра призмы, не принадлежащего ее треугольной грани.
Для оставшихся сочетаний индексов компоненты имеют вид:
к*г =
(КхЬЬг + Кускс. )Н К-А
24А 12Н
Компоненты составляющей матрицы теплопроводности [ К ]а приведены ниже.
1) При наличии конвективного теплообмена, характеризуемого коэффициентом теплоотдачи а, для прямоугольных граней:
Н - ь
I _ п ^иу
кП( ГС
9
а;
здесь s = (г,у,т,I) и (у,к,п,т) и (г,к,п,I),
ьиу - длина грани в плоскости иу;
18
-а; s = г,у,к;г = I,т,п; (*,г)е б*;
для всех остальных компонентов матрицы [ К ]а: НЬ
к=
па, sr
36
а.
2) При наличии конвективного теплообмена для треугольных граней призмы:
А
ка,** = V'а; s = (г, у, к) и (/, т, п).
6
. = А
ка,sг = тт - а; s,г = г,у,к или s,г = I,т,п.
12
Матрица теплоемкости для пятигранника имеет вид:
[С ]:
срАА
18
1 0,5 0,5 0,5 0,25 0,25
0,5 1 0,5 0,25 0,5 0,25
0,5 0,5 1 0,25 0,25 0,5
0,5 0,25 0,25 1 0,5 0,5
0,25 0,5 0,25 0,5 1 0,5
0,25 0,25 0,5 0,5 0,5 1
(7)
Ненулевые составляющие векторов тепловой нагрузки, обусловленные тепловым потоком от источников тепла {/ц} и конвективным
теплообменом { /а} имеют вид:
1) для четырехугольных граней:
Л,* = 0,25 - ц- 4у- Н;
/a,s = 0,25 - Т¥а- ьиу - Н;
* = г, у, к, I, т, п.
2) для треугольных граней:
/ц,* = 1/3 - ц - А; /а,* = 1/3 - Т„а- А.
Несмотря на простоту формирования соответствующих конечно-элементных матриц, эта процедура требует относительно длительных вычислений из-за решения кратных интегралов (при решении интегралов была использована автоматизированная система символьных вычислений MathCad 7.0). Выбор прямых призм в качестве конечного элемента с прямоугольным основанием для шестигранника и прямоугольным треугольником для пятигранника объяснялся возможностью получения аналитического представления элементов конечно-элементных матриц. В противном случае, получение матриц выполняется численно. Это приводит к существенным расходам вычислительной мощности используемой вычислительной техники.
Для оценки достоверности полученных соотношений для конечно-элементных матриц пяти- и шестигранных призм были разработаны соответствующие алгоритмическое и программное обеспечения. Адекватность расчетов оценивалась сопоставлением расчетных и экспериментальных значений температур для различных типов станков.
В качестве иллюстрации практической реализации разработанной объемной тепловой модели приводится тепловой расчет шпиндельной головки многоцелевого станка мод. МС 12250 М1-2 высокой точности. На рис.2.а и 2.б. приведены расчетные схемы для шпиндельной головки /ШГ/. На рис.2.а. изображена расчетная схема ШГ в сечении передней опоры шпинделя. На рис.2.б. представлена расчетная схема ШГ соответствующая межопорной части шпинделя. Расчетная схема на рис.2.а. описывает корпус ШГ вместе с приливом под подшипники качения. Такая же схема используется для задней опоры шпинделя. Вторая схема включает только корпус ШГ. Области, иллюстрируемые на рисунках
прямоугольниками, в дальнейшем разбиваются в автоматизированном режиме на шестигранные призмы. Области треугольного сечения аппроксимируются пятигранными призмами. Тепловое нагружение осуществляется в области контакта наружных колец подшипников с приливами ШГ , условно изображенное кольцевым тепловым источником с мощностью тепловыделения Q1. Проведенные машинные эксперименты показали, что этого теплового источника совершенно не достаточно для обеспечения адекватного теплового состояния ШГ. Не учет дополнительных тепловых потоков как от опор Q2, так и от нагревающегося в процессе работы станка шпинделя Q3 приводит к существенно заниженному тепловому состоянию ШГ.
Так как не проводилось трибологических исследований для опор, то уточнение уровня температурной стабилизации выполнялось путем уточнения тепловых потоков Q2 и Q3. Кривизна температурных характеристик на верхней и боковой поверхностях ШГ в точках , где были установлены температурные датчики, уточнялась путем варьирования коэффициентов теплоотдачи на этих поверхностях.
В общем случае, для ШГ включающей вращающийся ШУ условия теплообмена с изменением частот вращения должны меняться. Однако, в данном сдучае, флуктуации конвективных потоков от переднего конца шпинделя , при относительно незначительном изменении частот вращения, были не существенны. Поэтому при переходе с одной частоты вращения на другую коэффициенты теплоотдачи принимались одни и те же: для боковых поверхностей ШГ - 30 Вт /(м2 К); для верхней поверхности ШГ - 35 Вт /(м2 К); для остальных поверхностей - 12 н 14 Вт /(м2 К). Учет частоты вращения в большей мере отразился на идентификации тепловых потоков. Мощность тепловыделения Q1, Q2 и Q3 на последний момент времени моделирования работы ШГ в течение 120мин для частоты вращения 2000 мин-1 составила - 9,924; 0,577 и 1,26 Вт, соответственно. Для частоты вращения 1400 мин-1 - 6, 0,344 и 0,774 Вт.
Решение системы (1) выполнялось двумя путями. Первый - численное интегрирование с применением формул численного дифференцирования для вектора температур Т (*). Второй - модальный подход с применением редук-
ции Ланцоша. Первый подход - традиционный, практическая его реализация подробно изложена в отечественной и зарубежной литературе. Второй подход, наибольшее распространение получил в зарубежных системах моделирования, в основном, в динамике технических систем. Модальный подход в теории теплопроводности более известен под названием , отражающим алгоритм решения задачи нестационарной теплопроводности, как метод собственных значений. Решение задачи о собственных значений в модальном подходе , в основном определяющей временные затраты на вычисления, порождает как его достоинства, так и недостатки. Достоинства связаны с большей информативностью решения. Основными недостатками являются большие затраты вычислительной мощности, обусловленные необходимостью хранения разреженных и не разреженных матриц, и чрезмерным увеличением времени вычислений при росте размерности решаемой системы уравнений (1). Для уменьшения этих недостатков используется процедура редуцирования Ланцоша. Применение редукции Ланцоша позволяет перейти от решения задачи о собственных значений большой размерности к значительно меньшей. Размерность задачи о собственных значениях, при использовании редукции Ланцоша, соответствует числу учитываемых векторов Ланцоша.
На рис. 3. приведены четыре семейства расчетных и экспериментальных кривых в системе координат “температура Т - время ^». Кривые 1,2 и 3,4 иллюстрируют показания двух температурных датчиков для двух частот вращения шпинделя - 2000 и 1400 мнн 1 ( кривые 1,3 - показания датчика, установленного на верхней поверхности; кривые 2,4 - показания датчика на боковой поверхности ШГ). Датчики устанавливались на верхней и боковой поверхностях ШГ в области расположения передней опоры шпинделя. Кривые И, ^25; 2^ 2m35; 3^ 3m50; 4^ 4m10 иллюстрируют расчетные температурные характеристики, полученные при численном интегрировании системы (1) (в обозначении присутствует символ 0 и при использовании модального подхода совместно с редукцией Ланцоша (в обозначении использован символ m, цифры после m указывают на количество учитываемых векторов Ланцоша). Из рисунка видно, что расхождения экспериментальных и расчетных температур на
каждой из рассмотренных частот вращения шпинделя настолько несущественны, что в отдельные интервалы времени визуально неразличимы. По этой же причине на рисунке не приведены температурные характеристики одной и той же узловой точки при учете различного числа векторов Ланцоша. Более существенные различия наблюдаются лишь в первые моменты времени моделирования.
Сопоставление результатов расчета традиционным методом и на основе модального подхода выявило еще две отличительные особенности.
1) При традиционном подходе наблюдаются осциляции решения в некоторых точках температурного поля рассматриваемой термодинамической системы; при модальном подходе за счет реализации численно-аналитического подхода этой проблемы не существует. Так на рис.4. приведены температурные характеристики ШГ в наиболее нагретой точке. Кривые 1 и 2 соответствуют расчетным температурным характеристикам, полученным численным интегрированием системы (1), а ^ и 2m - через реализацию модального подхода.
2) Анализ всего температурного поля ШГ, т.е. теплового состояния всех узловых точек расчетной схемы в каждый моделируемый момент времени, показал, что при решении двумя подходами системы (1) в отдельных точках наблюдались существенные количественные расхождения (это видно из рис.4.). Наиболее существенным это расхождение было в местах действия тепловых источников. На других поверхностях, особенно находящихся в условиях конвективного теплообмена, оно сглаживалось и
становилось сопоставимым с точностью измерения температур в ходе натурных испытаний станков, не превышая ± 0,5° С.
Выполненные исследования тепловой модели станка позволили сформировать ряд рекомендаций по совершенствованию теплового моделирования станков:
- Тепловое состояние станка, рассчитанное двумя альтернативными подходами имеет некоторую погрешность. При анализе теплового состояния станков , не относящихся к группе современных сверхточных станков, такая погрешность расчетов может считаться приемлемой.
- Успешное уточнение тепловой модели по экспериментальным данным , означает, что применение призматических конечных элементов позволяет строить адекватные тепловые модели.
- Наличие осциляций температуры в некоторых узловых точках рассматриваемой расчетной схемы при использовании численного интегрирования, показывает достоинства модального подхода, как подхода обеспечивающего наибольшую адекватность моделирования. Тем не менее, максимальная скорость вычислений, бесспорно, на стороне численного интегрирования. Поэтому, при использовании призматических конечных элементов можно предложить двухступенчатое моделирование. На первой ступени моделирования - грубое моделирование -численное интегрирование и модальный подход с небольшим количеством векторов Ланцоша. Приоритетным показателем эффективности моделирования в этом случае выступает производительность моделирования. На второй ступени моделирования - модальный подход с увеличенным числом векторов Ланцоша.
Рисунок 1.а - Шестигранный призматический элемент Рисунок 1.б - Пятигранный призматический элемент
а,
МММ
а
а,
а
а
а
а
а
Рисунок 2.а - Расчетная схема шпиндельной головки Рисунок 2.б - Расчетная схема шпиндельной головки в в сечении передней опоры шпинделя сечении межопорной части шпинделя
Рисунок 3 - Экспериментальные и расчетные температурные характеристики в двух фиксированных точках на шпиндельной головке
О 20 40 60 80 100 і,мин
Рисунок 4 - Температурные характеристики наиболее нагретой узловой точки расчетной схемы шпиндельной головки
