CC BY
34
А.Н.Поляков
УЧЕТ ФЛУКТУАЦИЙ УСЛОВИЙ ТЕПЛООБМЕНА В ВЕРОЯТНОСТНОЙ ТЕРМОУПРУГОЙ МОДЕЛИ МЕТАЛЛОРЕЖУЩЕГО СТАНКА
Приведен частный случай вероятностной термоупругой модели металлорежущего станка. В качестве системы случайных величин принимались коэффициенты теплоотдачи. При построении модели использована теория линеаризации функций случайных величин. Показано, что главный критерий эффективности практической реализации процедуры линеаризации связан с точностью нахождения соответствующей производной, поэтому выполнен анализ нескольких схем расчета производных. При этом, эффективность решения системы линейных уравнений обеспечивает качество вероятностной термоупругой модели. Учитывая большую размерность решаемой системы, была предложена новая компактная схема хранения разреженных матриц, и выполнен анализ различных численных методов решения систем линейных уравнений.
Развитие термоупругой модели металлорежущего станка, объединяющей две ее составляющие - тепловую и упруго-деформационную модели, неизбежно связано с совершенствованием ее структуры. Одним из направлений совершенствования математических моделей в станках, является переход от детерминированного описания к вероятностному [1].
Проведенный анализ существующих термоупругих моделей станков показал, что наиболее распространенной является модель, в которой числовые характеристики и закон распределения температуры устанавливается либо в ходе проведенных статистических машинных, либо экспериментальных исследований [2]. В данной работе используется аналитическое представление числовых характеристик температуры и температурных деформаций. Положительным свойством этого, является возможность получения более точного и быстрого решения. Негативным, сужающим область моделирования - предположение о независимости всех участвующих в модели случайных величин. В качестве случайных величин в данной работе принимаются лишь флуктуации условий конвективного теплообмена, количественно выражающиеся через коэффициенты теплоотдачи.
Математическое ожидание для температуры принимается с учетом методов линеаризации для функций случайных независимых величин в виде [3]:
М[Т,] = т = 1/Лк • £ р* £р/(1 - е~Лк‘) +Ъ<р-1киоке-Хк (!) к=1 ,=1 к=1
здесь I - подстрочный индекс, соответствующий нумерации узла системы; п - число узловых элементов системы; /, - элемент
вектора тепловой нагрузки; модальная матрица [Ф], столбцами которой являются собственные векторы {pk}, соответствующие собственным значениям 1k, полученным в результате решения обобщенной задачи о собственных значениях:
[ K ]{pk} = Ik [C ][p k ] (2)
где [K],[C] - матрицы теплопроводности и теплоемкости; uok —температура в нормальных координатах в начальный момент времени t=0; {uJ находится из решения системы уравнений: [Ф]{uo} = {To}, где {TJ - вектор температур в физических координатах в начальный момент времени. В уравнении (1) участие случайных величин представлено модальной матрицей и вектором собственных значений, вследствие зависимости матрицы теплопроводности от коэффициентов теплоотдачи. Равенство математического ожидания и узловой температуры в каждом i-ом узле обеспечивается расчетом температурного поля для термодинамической системы с условиями конвективного теплообмена, равными их математическому ожиданию.
Выражение для дисперсии будет иметь
вид:
D[T(a)] = (|T)2. D[a] (3)
где D[a] - дисперсии коэффициентов теплоотдачи.
Математическое ожидание вектора тепловых деформаций определяется по тем же зависимостям, что и для детерминированной модели [4].
Обычно вектор тепловых деформаций определяется из решения задачи статики. Вместе с этим, может быть использован и
модальный подход, применяемый в тепловой модели. Это означает, что система для определения тепловых деформаций в квазистаци-онарной постановке будет иметь вид:
[ Е ]{5} + [К ]{5} = }, (4)
здесь [ Е ] - единичная диагональная матрица {*} 85
и {5} = —- - вектор производных узловых пе-дт
ремещений, [К] - матрица жесткости, {^} -вектор тепловой нагрузки.
Для решения системы (4) как и в работе [5], используется двойное преобразование, с помощью которого осуществляется переход от физических координат {5} к модальным {#} {5} = [¥]№]{*} (5)
здесь [¥] - матрица нормированных векторов Ланцоша, [ф] - модальная матрица, столбцами которой являются соответствующие собственные векторы системы (4).
Таким образом, вектор перемещений {5}, обусловленный действием тепловой нагрузки, для установившегося режима будет иметь вид:
{5} = [У][Ф][Ф]Т {/} / ) к=, ...„, (6)
где {/т} = №4к]-1{^} - вектор нагрузки в обобщенных координатах, полученных при использовании матрицы преобразования [¥], Ма%(Хк )к=1..т - диагональная матрица собственных значений, т - число векторов Ланцоша, участвующих в модели, равное числу собственных векторов и значений, к - текущий индекс собственного значения.
Уравнения (5) и (6) представляют реализацию редукции Ланцоша для упруго-деформационной составляющей термоупругой модели станка. При практической реализации редукции Ланцоша всегда отмечают проблему выбора количества векторов Лан-цоша, и потерю их ортогональности, вследствие роста погрешностей вычислений, неизбежно возникающих при использовании ЭВМ. Но проблема выбора количества векторов Ланцоша возникает лишь при необходимости получения более точного решения для неустановившегося режима. Этой проблемы не возникает, если требуется получить установившееся решение. Этот факт был отмечен еще при подготовке работы [5]. Кроме того, выбор наименьшего количества векторов Ланцоша позволяет минимизировать вычислительные затраты, связанные с реализацией редукции Ланцоша. Поэтому был выбран один вектор Ланцоша, т.к. основным назначением редукции Ланцоша, в данной модели, являлось получение аналитического выражения для вектора тепловых деформаций в форме, удобной для учета влияния на
него флуктуаций условий конвективного теплообмена.
Характерной процедурой для термоупругой модели является решение системы линейных уравнений. Эта процедура используется как в традиционных моделях, не использующих модальный подход, так и в редукции Ланцоша [5]. Несмотря на кажущуюся простоту этой процедуры, в силу вычислительных свойств рассматриваемой системы, ее практическая реализация для систем большой размерности представляет определенную проблему. Во-первых, изменяются требования к организации, хранению и обработки соответствующих матриц термоупругой модели. Так, например, хранение только одной матрицы жесткости с рангом не менее 10000 традиционным способом, в виде двухмерного массива, потребовало бы оснащение компьютера ОЗУ емкостью равной не менее 400 Мб. Во-вторых, прямые методы решения становятся либо недоступны, в силу существующих ограничений вычислительных мощностей, либо малоэффективны. При этом эффективность оценивается тремя показателями: размерностью решаемой системы, точностью и быстродействию вычислений. Свойство сильной разреженности матрицы жесткости и теплопроводности используют при организации, хранении и обработки матриц. На сегодняшний день существуют большое многообразие компактных схем хранения разреженных и симметричных матриц. Вместе с этим, при практической реализации компактных схем хранения матриц, были учтены опыт накопленный автором при разработке конечно-элементного пакета ТЕШОБ, а также стремление поиска наилучшего сочетания между организацией и компактностью хранения матриц и эффективностью их обработки. Более компактная схема хранения матрицы, как правило, сложней реализуется в алгоритмах ее обработки. В данной работе была реализована следующая схема хранения матрицы:
1) все ненулевые элементы матрицы хранятся в виде двухмерного массива, для которого число строк совпадает с рангом матрицы, а число столбцов равно числу ненулевых элементов каждой строки (для задачи статики, при рассмотрении плоско-напряженного состояния и изгиба это число не превышало сорока одного, для задачи теплопроводности еще меньше);
2) вводится дополнительно два целочисленных массива, необходимых при последующей обработки этой матрицы. Первый -одномерный массив, определяет количество ненулевых элементов в каждой строке, а второй - двухмерный, обеспечивает соответствие элементов матриц компактной схемы и традиционной, размера [К х К].
Такой способ организации матрицы обеспечивает преемственность и простоту модификации алгоритмов, применяемых для матриц, хранящихся традиционным способом. Вместе с этим, проведенные вычислительные эксперименты показали, что при реализации прямых методов решения системы линейных уравнений (в данной работе анализировалась эффективность методов Гаусса, Холецкого и метода отражений), при увеличении размерности решаемой системы, снижается эффективность организации матрицы по компактной схеме. Это объясняется тем, что алгоритмы всех прямых методов предусматривают появление промежуточной информации, что неизбежно увеличивает число ненулевых элементов матриц. Поэтому для систем большой размерности, таких как металлорежущий станок или отдельные его узлы, применяют либо методы редукции [6] либо итерационные.
В настоящей работе был проведен анализ пяти итерационных методов решения систем линейных уравнений: метод простых итераций, метод Зейделя, метод верхней релаксации, явный итерационный метод с че-бышевскими параметрами (чебышевский итерационный метод) и метод сопряженных градиентов. Для тепловой модели любой из названных методов позволял получать устойчивые результаты. (, с гарантированной точностью до третьего знака после запятой.) Для упруго-деформационной модели - только метод сопряженных градиентов. В первую очередь это объясняется различием в структуре матриц жесткости и теплопроводности (для тепловой модели). Матрица теплопроводности удовлетворяет требованиям диагонального преобладания и хорошей обусловленности, что не характерно для матрицы жесткости. На рис.1. приведены результаты расчетов шести систем коробчатых конструкций различной размерности (420, 647, 1173, 1855, 2693 и 3687 узлов соответственно) пятью методами на ПК на базе процессора Pentium 233 ММХ: метод Гаусса - кривая 1; метод простых итераций - кривая 2; метод Зейделя - кривая 3; метод последовательной верхней релаксации с параметром релаксации w = 1,4 - кривая 4; метод сопряженных градиентов - кривая 5. Отметим, что расхождение всех результатов расчетов наблюдалось лишь в отдельных узловых точках, начиная с третьего знака после запятой. Анализ полученных результатов показал, что для хорошо обусловленной матрицы теплопроводности метод сопряженных градиентов и метод последовательной верхней релаксации незначительно отличаются по скорости сходимости, причем значительно опережая другие рассмотренные методы. Вместе с
этим, проведенные исследования для упруго-деформационной модели показали, что метод сопряженных градиентов позволил получить наиболее точные решения. Максимальные расхождения решения в сравнении с решением, полученным методом Гаусса, в отдельных узлах не превосходили 1%. При этом время решения систем сопоставимых размерностей для тепловой и упруго-деформационной моделей значительно отличаются. Так время расчета системы, содержащей 3887 обобщенных координат составило 17 с (секунд), что более чем в три раза превышает временные затраты на решение тепловой модели с размерностью 3687 (сопоставимой по размерности, рис.1.). Применение метода Зейделя позволяло получать решения с расхождением до 8%. Однако, по отдельным обобщенным координатам, когда перемещение составляло менее 1 мкм, это расхождение составляло почти 100%. При этом скорости сходимости методов Зейделя и сопряженных градиентов уже значительно отличались. Так для решения системы выше приведенной размерности методом Зейделя потребовалось 648 с - это почти в сорок раз больше чем при реализации метода сопряженных градиентов. Реализация метода последовательной верхней релаксации показала, что хотя метод обладает лучшими свойствами по скорости сходимости, чем метод Зейделя, но необходимость экспериментального определения параметра релаксации ■№, в значительной степени сужает область его практического применения. Аналогичные суждения можно высказать и для итерационного метода с чебы-шевскими параметрами. При удачном их сочетании, в тепловой модели, можно достигнуть даже большей скорости сходимости метода, чем в методе сопряженных градиентов, но этот же факт и определяет неустойчивость данного метода (поэтому результаты расчетов этим методом не нашли отражение на рис.1.). Фактор обусловленности системы дополнительно ограничивает область применимости методов. Как показали проведенные вычислительные эксперименты, чем хуже обусловленность системы, тем ниже скорость сходимости метода. Более того, в отдельных случаях оказывается невозможным использовать метод верхней релаксации и чебышев-ский итерационный метод для получения решения с заранее заданной погрешностью. Для этих методов также следует отметить наличие зависимости между величиной параметров, влияющих на их сходимость, и вектором нагрузки. Поэтому, для систем, отличающихся лишь вектором нагрузки, необходимо выбирать параметры различающиеся по величине. Все эти факты не позволили использовать ни метод Зейделя, ни метод релакса-
ции, ни итерационный метод с чебышевски-ми параметрами в последующих алгоритмах вероятностной термоупругой модели. Отметим, что проверкой сходимости выбранного итерационного метода являлось решение задачи о собственных значений вида:
Л}=4ЛТ ИЛ}, (7)
здесь [Г ] = тт [К гт - трехдиагональная матрица Рэлея, {лк} - к - й собственный вектор модальной матрицы [Ф], 1 - к - ое собственное значение. Из (7) видно, что наличие решения задачи о собственных значениях обеспечивается свойствами матрицы Рэлея [Т. ], полученной в ходе реализации процедуры формирования векторов Ланцоша. Из всех рассмотренных методов лишь метод сопряженных градиентов оказался наиболее приемлемым, т.к. его алгоритм не предполагает дополнительных исследований для получения решения. В качестве иллюстрации практической реализации данного метода отметим, что при расчете станины двустороннего торцешлифовального станка с дискретностью конечно-элементной сетки, обеспечивающей наличие в расчетной схеме 22122 узла, только метод сопряженных градиентов позволил найти решение задачи (7), при этом, время решения системы линейных уравнений, сопутствующей процедуре формирования матрицы Рэлея, составило 320 с.
Получив, согласно уравнениям (4) - (6), математическое ожидание вектора тепловых деформаций, вычисление его дисперсии 0[8 (а)] осуществляется по формуле:
Л[д(а)] = £(дд(а))2 • Я[а], (8)
й да
здесь Б - число поверхностей, по которым наблюдается возмущение условий конвективного теплообмена; Д[а,. ] - дисперсия коэффициентов теплоотдачи по /-ой поверхности; д - оператор дифференцирования.
При использовании аналитического выражения (6), для расчета дисперсии Б[д (а)], как сложной функции, необходимо вычисление производных для функций случайных величин, не имеющих аналитического описания: собственных векторов, модальной матрицы и матрицы векторов Ланцоша. Единственный способ для получения их производных - это использование интерполяционных формул дифференцирования. Для оценки точности вычисления производной
дд (а)
—г— производился вычислительный экспе-да
римент. В ходе этого эксперимента сопос-
тавлялись результаты прямого расчета (для сокращения изложения в последующем будут называться оригиналом) тепловых деформаций д (а0 +Аа) и спрогнозированной величины тепловых деформаций, полученной с использованием разложения в ряд Тейлора, учитывающего лишь линейные члены. При этом эксперимент проводился для шпиндельных узлов различных конструкций для шлифовальных и токарных станков, станины и шпиндельной бабки двустороннего торцешлифовального станка [4]. Для всех конструкций тепловая модель уточнялась, путем сопоставления результатов расчета с данными, ранее проведенных, натурных экспериментов. Для двустороннего торцешлифовального станка имелась возможность провести идентификацию термоупругой модели, т.к. помимо данных о температуре в характерных точках станка, имелись экспериментальные данные по деформациям, измеренным на холостом ходу. Проведенный анализ полученных результатов показал, что при возмущениях условий конвективного теплообмена в размере ± 20%, можно гарантировать достаточно высокую точность вычисления производной. Это подтверждалось тем, что расхождение между спрогнозированной величиной деформации д(а0 +Аа) и оригиналом не превышало 1%. При увеличении возмущения этот показатель ухудшался. С целью увеличения точности прогнозирования, а значит расчета производной и в конечном итоге дисперсии, был выполнен анализ функциональной зависимости вектора тепловой деформации от изменения коэффициентов теплоотдачи. Оказалось, что хотя функция д(а) имеет нелинейный характер, но для относительно небольших возмущений аргумента, до ± 50%, можно использовать линейную аппроксимацию. В данной работе аппроксимация выполнялась тремя способами. Иллюстрацией линейных аппроксимаций служат графики 2
- 4 (рис. 2.), являющиеся результатом вычислительных экспериментов. Кривая 1 иллюстрирует функцию д(а) - оригинал. Прямая 2 показывает результаты расчета спрогнозированной функции д (а), полученной разложением в ряд Тейлора, в котором производная получена как производная аналитического выражения деформаций (6) (для краткости последующего изложения результаты этого расчета будут называться модальный прогноз). Как следует из характера графика 2 - чем больше диапазон флуктуаций коэффициентов теплоотдачи системы, тем значительней расхождение между спрогнозированной величиной д (а0 +Аа) и оригиналом.
Чтобы уменьшить это расхождение, существует несколько способов построения соответствующих линейных аппроксимаций, минимизирующих это расхождение. Первый способ - это аппроксимация кривой 1 прямой 3, полученной по уравнению прямой, проходящей через две точки. Первая точка
- это тепловая деформация, вычисленная для математического ожидания коэффициентов теплоотдачи д (а о). Втор ая точка соответствует деформации, вычисленной для граничных значений коэффициентов теплоотдачи д (а0 +Аа). В этом случае отличие прогнозируемой величины деформации от оригинала будет меньше, чем в первом варианте расчета. Без ограничения общности изложения, для прямой 3 можно выбрать точки д(ао -Аа) и д (а0 +Аа). Второй способ аппроксимации отличается от первого тем, что положение обеих точек, через которые проходит аппроксимирующая прямая, не обязательно являются граничными. В этом случае точность расчета незначительно отличается от первого варианта расчета. При использовании первого и второго способов аппроксимации деформаций, производная принимается равной коэффициенту при аргументе в уравнении соответствующей прямой. Третий способ аппроксимации (прямая 4) отличается от модального прогноза использованием более точного значения спрогнозированного температурного поля. Так, если для модального прогноза модель для прогноза температуры и деформаций идентичны, то в этом случае прогноз температурного поля осуществлялся по модели аналогичной описанной в первом и втором способах аппроксимации деформаций. Для полноты представления о приведенных способах расчета производных, необходимо отметить, что третий способ аппроксимации и модальный прогноз для упруго-деформационной модели для систем большой размерности отличают большие временные затраты, чем для первого и второго способов аппроксимации.
Укрупнено, алгоритм построения термоупругой модели с альтернативными методами расчета производных показан на рис.3. Представленный алгоритм представляет собой последовательное построение двух фактически самостоятельных моделей, объединенных лишь набором передаваемых параметров. Несмотря на независимость построения двух моделей, особенностью данного
алгоритма является идентичность методов, используемых при построении обеих моделей. Весь алгоритм включает восемь блоков. Первые четыре блока 1- 4 отвечают за формирование вероятностной тепловой модели, а блоки 5 - 8 - за формирование упруго-деформационной модели. Блоки 1 и 5 отвечают за расчет математического ожидания вектора температур и вектора тепловых деформаций. Изменяющееся во времени температурное поле станка порождает изменяющееся во времени тепловое деформационное состояние, поэтому, блок 5 содержит помимо традиционной процедуры решения задачи статики и модальный расчет (аналогичный блоку 1). В блоках 2 и 6 выполняется расчет производных соответственно температур и деформаций, представленных в аналитической форме как функции случайных величин. В блоках 3 и 7 выполняется расчет производных альтернативными методами, через аппроксимацию прямой, проходящей через две точки. Завершают построение моделей расчеты Д[Г(а>], Д[£ (а)] в блоках соответственно 4 и 8. Блоки сравнения А? и В? дают право выбора метода расчета производной.
В ы в о д ы
1. Предложенный вариант построения вероятностной термоупругой модели отличается простотой и гарантирует достаточно высокую точность вычисления числовых характеристик. Важной особенностью предложенного подхода является наследственность методов, используемых в детерминированных моделях, что при практической реализации данной модели обеспечивает наименьшую модификацию базового программноматематического комплекса.
2. Установление функциональных связей между входными параметрами тепловой модели и выходными параметрами упруго-деформационной модели расширяет диапазон решения прикладных задач теории вероятностей.
3. Реализация термоупругой модели для систем большой размерности, к которым относятся металлорежущие станки, предъявляет особые требования к выбору соответствующих алгоритмов. Предложенная схема хранения разреженных матриц, обеспечивает достаточно высокую компактность хранения и простоту модификации алгоритмов, рассчитанных на различные варианты организации матриц.
Список использованной литературы
1. Поляков А.Н. Реализация вероятностного подхода в тепловом моделировнии металлорежущих станков // Перспективные материалы, технологии, конструкции: Сб. науч. тр./ Под ред. проф. В.В. Стацуры. - Вып.4. - Красноярск: САА, 1998. - С. 164-168.
2. Пуш А.В. Шпиндельные узлы. Качество и надежность. - М.: Машиностроение, 1992.- 288 с.
3. Хомяков В.С., Досько С.И., Поляков А.Н. Применение теоретического модального анализа к расчету температурных полей в металлорежущих станках // Изв. Вузов. Машиностроение. 1989.- №9.- С.154-158.
4. Поляков А.Н., Никитина И.П. Применение термоупругой модели к анализу тепловых процессов в металлорежущих станках. // Вестник машиностроения. - 1996. - №7.- С.27-30.
5. Поляков А.Н. Применение метода Ланцоша к построению нестационарного температурного поля в металлорежущих станках // Вестник машиностроения. 1997. № 10. С.43 - 46.
6. Поляков А.Н. Применение редукции Г аяна с выбором рационального базиса к построению нестационарного температурного поля в металлорежущих станках // Изв вузов. Машиностроение. - 1997. - №4-6. - С.109-116.
Статья поступила в редакцию 31.08.99г.
