Моделирование нестационарных температурных полей с использованием совместных конечноэлементных схем Текст научной статьи по специальности «Физика»

УДК 539.3

МОДЕЛИРОВАНИЕ НЕСТАЦИОНАРНЫХ ТЕМПЕРАТУРНЫХ ПОЛЕЙ С ИСПОЛЬЗОВАНИЕМ СОВМЕСТНЫХ КОНЕЧНОЭЛЕМЕНТНЫХ СХЕМ

А.А. Мишин

Разработана программа, моделирующая нестационарный процесс теплопроводности в системе тел, для конечноэлементных схем которых справедлив принцип попарного соответствия контактных узлов. Для решения системы линейных уравнений методом Холецкого разработан специальный способ хранения коэффициентов глобальных матриц

Ключевые слова: температурная задача, глобальная матрица

Одним из направлений развития современного транспорта является повышение скоростей движения, что позволяет повысить пропускную способность магистралей, и как результат, получить выигрыш по времени. Но возрастание скорости квадратично связано с возрастанием кинетической энергии, что ужесточает требования к динамике подвижной единицы, к тормозным системам. Принцип действия последних основан на преобразовании кинетической энергии в любой возможный другой вид энергии, чаще всего в тепловую, а для машин с электрической тягой ещё и в электромагнитную. Таким образом, тепловая нагруженность узлов трения тормозов является актуальным вопросом исследования. Повышение скорости движения требует повышение мощности двигателей, в том числе внутреннего сгорания. Так для дизеля, применение турбокомпрессора позволяет повысить литровую мощность двигателя примерно в два раза. Это достигается в основном за счёт сжигания большего количества топлива в том же объёме, что приводит к повышенной тепловой нагруженности цилиндро-поршневой группы. Далеко не последнюю роль в зарождении и дальнейшем развитии трещины в сварных швах являются температурные деформации и напряжения, которые могут быть вычислены по известному температурному полю. Работоспособность сопел ракет, лопаток турбин, печей, котлов связана с их конструкцией, учитывающей потоки подвода и отвода теплоты и распределение температур.

Уравнение теплопроводности в конечноэлементной постановке

Классическая теория теплопроводности исходит из пропорциональности градиента температуры мощности потока:

dQ - . dT

-n = &-

dtdS

dQ

dtdS

dn

либо в скалярном виде

= л gradnT\

где X -коэффициент теплопроводности; и зависимости между изменением температуры массы ёш на величину ёТ путём поглощения массой количества теплоты dQ :

dQ = сёшёТ, либо dQ = срёУёТ ,

где с-удельная теплоёмкость вещества, р-плотность вещества.

Составление уравнения теплового баланса

приводит к взаимодействию этих двух

взаимопереходящих друг в друга состояний теплоты

во времени, что выражается дифференциальным

уравнением нестационарной теплопроводности в

однородной среде:

дТ ____

ср — = div(Agra<ST), (1)

д/

что вместе с начальными и граничными условиями второго и третьего рода, даёт корректно поставленную задачу:

дТ ____

ср— = div(ЯgradT) д/

T (х, y, z)\ t=0 = f (х, y, z) l|grad(T) = q -2|grad(T) = f (T - T0)

(2)

q =

d2 Q dSdt

Мишин Алексей Александрович - БГТУ, аспирант, тел.(4832) 56-86-37, Email: mish7@yandex.ru

Метод конечных элементов (МКЭ)

представляет искомую функцию и в области О как набор кусочно-непрерывных функций и в подобластях ю, называемых конечными элементами. Характерные для задачи физические свойства подобластей ю выражаются интегральными характеристиками в фиксированных точках облсти О, называемых узлами конечных элементов. Точки выбираются из условий крутизны изменения ив О, а также в точках, где заданы граничные условия, интегрально приведенные к узлам. При

математическом моделировании задачи

нестационарной теплопроводности функции и возможно соответствие температурного поля Т , а граничным условиям - тепловых мощностей.

Кусочно-непрерывные функции и

аппроксимируют ю, то есть в узлах значения и и и равны. Следовательно, можно составить количество равенств для и не превышающих число узлов, им же будет ограничено и число независимо определяемых коэффициентов для и . При решении пространственной задачи в декартовых координатах самым простым видом функции и можно назвать

полином. Решив задачу относительно его

коэффициентов, получим выражение для и относительно её узловых значений, множители при которых в МКЭ называют функциями форм. Таким образом, значение и в произвольной точке ю будет суммой произведений функций формы узлов на их узловые значения. Применительно к температурной задаче:

Т(X, у, 2) = ЫТ (X, у, 2)Туз,

где ЫТ (х, у, 2) - вектор-строка функций форм,

Туз - вектор-столбец узловых температур.

Тогда градиент температуры в пространственных декартовых координатах равен:

'дЫ дЫ Ы

дх ду дх

Учитывая, что правый и левый члены уравнения (1) - объёмная плотность мощности теплоты, а члены граничных условий (2) поверхностная плотность мощности теплоты, запишем интегральное уравнение баланса мощностей для тела в исходном виде:

ЯгасТ (х, у, і) — Гтуз

ГГГ дТ д

^ ср — + —

О д дх

Хх

дТ

дх

д

ду

У

дТ

ду

дТ Л

ді )

= Ц Ч - Г(т - Т0 К £

Поскольку О состоит из непересекающихся ю, то значение интеграла для О будет суммой интегралов для ю, таким образом, для конечного элемента, имеющего приток и отток теплоты получим:

Ш ср-+ ~{

ді дх V

дТ

дх

д ( дТЛ д ( дТЛ

+-----1 Х I +------------1 Х шю —

ду V у ду ) ді V 2 ді )

— Ц Я - Г (Т - Г)*.

Подставив в уравнение выражение для градиента, получим баланс мощностей для конечного элемента через узловые значения температуры, её производные по времени и суммарным потоком тепла на поверхности, но прежде необходимо эту мощность разнести по узлам конечного элемента. Для этого умножим подынтегральные выражения справа на N(х, у, г):

гг/ дТ д ( дТ Л д ( ЮюСР ді + дх Х дх ) + ду [

дТ д ( _ дТ ) д ( дТ ) д ( дТ , ,

| + | Ху I + І Хі I I —

ду) ді V ді

= Л(? - т (Т - Т0 ^

После интегрирования и преобразований получим:

дТ„

cрNN у + уХх XX + Ху УУ + Хі 22

__У3_

ді

= + Т0 Шзг ,

)+ К

где

NNV — |ЦNNTdV , NS — ЦШЗ , NNS — ЦNNT ЫБ V Б Б

XX — УУ —

V

дх дх

22 — ^.

V ді ді

Для изотропного упрощается:

V ду ду

материала выражение

cрNNv -

дТуз

ді

- + \AXYZ + hNNS

Т

г уз

Зс г уз

где

— + Т0ШЗ,

Я О ~'с

— XX + X + XX ; Х — Хх —Х у — Хі ■

Пусть матрица С — cрNNV - глобальная

а М — ?ХУ2 + hNNS -

матрица теплоёмкости,

глобальная матрица теплопроводности, а

Р — qNS + Т0 hNS - глобальный вектор узловых

Я с

мощностей, тогда уравнение нестационарной теплопроводности для метода конечных элементов запишется в виде:

дТ „

С—^ + МТз — Р .

ді у

Применив неявный метод Эйлера как абсолютно устойчивый при любом шаге Ді Ф 0 для аппроксимации первой производной, получим систему линейных уравнений для приращения температуры на і -том шаге:

С

Ді

+М

,дт; —-м • Т; + р

Т,+1 — Т‘ +ДТі

уз уз уз

Представление глобальной матрицы как многосвязного списка

Два свойства глобальных матриц дают возможность поиска специальных методов хранения, адаптирующих матрицу к специальным операциям и сокращающих объём выделяемой памяти: первое - симметричность относительно главной диагонали, второе - сильная разреженность. Поскольку решение температурной задачи сводится к решению системы линейных уравнений, то выбор метода решения может значительно повлиять на выбор организации хранения и доступа к коэффициентам. Ещё одно свойство глобальных матриц МКЭ предопределило выбор метода решения СЛАУ - это диагональное преобладание и отсутствие нулевых коэффициентов на главной диагонали, все эти свойства матрицы как нельзя лучше подходят для прямого метода Холецкого. В процедуре разложения матрицы системы и последующего решения прямым и обратным ходом требуется последовательный доступ к элементам как строки, так и столбца. В случае классического представления матрицы в языке С/С++ как массива указателей на строки, доступ к у -тому элементу строки осуществляется смещением по массиву указателей на / -тую строку и последующему

+

Х

х

смещению по строке на у -тый столбец. Таким образом, последовательный доступ к элементам строки осуществляется движением по адресному пространству с шагом, равным размеру типа данных массива, а шаг доступа к элементам столбца будет зависеть от способа хранения элементов в строке, если матрица обычный двумерный массив - шаг будет равен размера строки в байтах, то есть доступ к элементам строки и столбца можно считать равноправно-быстрым. Но хранить глобальную матрицу как двумерный заполненный массив, значит хранить и нули. Иметь преимущества доступа двумерного массива и сократить число хранимых нулей, к тому же использовать свойство симметрии матрицы, позволяет профильный метод хранения [1]. Полностью исключить хранение нулевых элементов позволяет компактное хранение матрицы [2], при котором вместе с элементом хранится номер строки и столбца, но доступ к элементам столбца возможен либо перебором элементов строки, либо хранением матрицы без учёта симметрии, в этом случае движение от диагонального элемента влево будет означать движение по столбцу верх, а вправо - вниз. Есть ещё одно обстоятельство, которое влияет на объём занимаемой матрицей памяти, это появление новых элементов при разложении. Использование метода Холецкого гарантирует отсутствие ненулевых элементов за пределами специального профиля, построенного по самым длинным строкам профиля матрицы. Если строка матрицы - массив, то добавление нового элемента приводит к перевыделению памяти, и копированию элементов в новую строку. Представление строки как списка позволяет добавить новый элемент, изменив лишь связи взаимодействующих с ним членов.

На основе принципа списочного представления данных в настоящей работе была разработана архитектура списочного хранения матрицы. Одна из её разновидностей, рис.1, полностью учитывает сильную разреженность матрицы, её симметрию, реализует равноправный последовательный доступ к элементам строки и столбца, который при использовании метода Холецкого, делает минимальными различия в скорости вычислений в сравнении с использованием индексированного доступа.

Элемент списочной матрицы, показанный на рисунке квадратом, содержит значение

коэффициента матрицы, номера строки и столбца (обозначено кругом), а также четыре указателя на соседние строковые и столбцовые элементы, (обозначены трапециями и стрелками). Матрицу из таких взаимосвязаных элементов назовём списочной. Списочная матрица содержит массив указателей на диагональные элементы, которые на рисунке показаны квадратами с заштрихованными кругами. Таким образом, доступ к произвольной строке остаётся индексированным, а значит быстрым. Заполнение матрицы ведётся поиском нужной позиции и установлением связей по месту.

Для организации итераций признаком конца строки и столбца является инициализированный нулём соответствующий указатель (на рисунке -

заштрихованные трапеции). Данная схема

представления данных может быть удобной при

Рис. 1. Схема организации коэффициентов глобальной матрицы как многосвязного

формировании суперэлементов, поскольку затраты на удаление и добавление новых элементов связаны преимущественно с их поиском, а результирующую заполненную матрицу суперэлемента всегда можно представить двумерным индексированным массивом.

Программная реализация

Для программной реализации решения температурной задачи методом конечных элементов описаным выше способом удобным будет язык, поддерживающий создание и работу с основными типами данных, средствами работы с памятью, объектами, шаблонами, файлами, библиотеками, в том числе и графическими, для которого существовали бы компиляторы под основные операционные системы, имеющие хорошие вычислительные свойства. Этим пожеланиям в полной мере удовлетворяет язык C++, компилятор под OS Windows MS Visual C++ 6.0 и графическая библиотека OpenGL.

Программа состоит из семи основных библиотек и пользовательского интерфейса. Структура программы построена по принципу тематического обособления и однонаправленного использования элементов, начиная от базового (рис.2). Первой в иерархии стоит библиотека Programm, в которой описаны и реализованы шаблоны классов список и n-мерный вектор, шаблон иерархии классов матриц размерности n х m, объявлены и инициализированы макроопределения и константы.

Второй в иерархии стоит математическая библиотека Mathematics. В ней находится класс кватернион, иерархии классов систем координат, трёхмерных векторов, матриц 3 х 3 . Также здесь создан класс «списочная матрица» и разработаны функции для работы с его объектами, такие как добавление и удаление элементов, сжатие профиля, копирование профиля из другой матрицы,

умножение на вектор, решение системы линейных уравнений.

Третья библиотека Geometry содержит классы геометрических объектов, как-то: отрезок, дуга,

Interfасе

Pr oblens

suai

Physical

G r i n d

Geonetr'

Mathenc

es

Program

Рис. 2. Структурная схема программы

(стрелкой показано обращение к библиотеке)

плоскость, обобщённая четырёхугольная область. Библиотека широко использует векторную алгебру библиотеки Mathematics.

Четвёртая библиотека Grind предназначена для разбивки геометрических объектов на конечные элементы, выделения особых областей разбивки, их узлов. Здесь же реализована иерархия классов конечных элементов. Предусмотрена функция считывания конечноэлементной разбивки из файла заданного формата, подготовленного сторонней программой.

Пятая библиотека Physical, используя конечноэлементное представление объекта, формирует по шаблонам из Programm вектора функций форм и производных от них. Полученные на их основе локальные матрицы XX, YY, ZZ, NN и другие, могут либо объединяться в соответствующую глобальную, либо сразу формировать нужную матрицу, учитывая данные о материале элементов. Библиотека содержит иерархию классов узлов, таких как температурный, упругий, контактный и их комбинации.

Шестая библиотека Visual содержит функции инициализации графики OpenGL, функции работы со сценой, классы изображения на экране необходимых объектов предыдущих библиотек. Используются функции рисования из библиотек opengl32.dll, glu32.dll.

Седьмая библиотека Problems является основной библиотекой, к которой подключаются все ресурсы предыдущих библиотек. Она содержит классы задач и характерные для них методы. Так класс, моделирующий нестационарный

температурный процесс, реализует итерации по времени, где на каждом шаге решается система линейных уравнений относительно температуры с новыми граничными условиями.

Пользовательский интерфейс содержит главное окно для вывода графической информации, меню, кнопки, панель инструментов для работы с

программой. Интерфейс написан под OS Windows в многодокументном исполнении, с архитектурой документ/представление. Класс документа содержит объект основного класса библиотеки Problems, через который ведётся создание и управление всеми объектами из семи библиотек, включая саму Problems.

Решение тестовых задач

Для проверки правильности моделирования тепловых процессов использовалось аналитическое решение задачи распределения температуры T ,0C в полупространстве, к границе которого подводится теплота с постоянной плотностью мощности Вт

q

м ■ К

qz

expl

+ erfl

2-Jat

2л[а

где

Я

а = —; ср

Я - коэффициент теплопроводности,

Дж ;

Вт м ■ К

с - удельная теплоёмкость,

кг ■ К

р - плотность,

кг

-.3

2 - координата точки, расположенной на внутренней нормали к плоскости полупространства; Т0 - температура при 2 ^ ж, 0 С ,

/ - время, с.

Конечноэлементная модель - прямоугольный параллелепипед со сторонами 70 х 70 х 280 , с дискретизацией: по двум направлениям - 10

конечных элементов (11 узлов) и вглубь полупространства - 40 элементов (41 узел). Итого схема содержит 4961 узел и 4000 конечных элемента (рис.3).

Для расчёта задавались (размерности в системе

СИ): ч = 106, 2 = 26, с = 490, р = 7800, / = 80,

Д/ = 5, Тз = 20 . При выбранном тепловом

нагружении полупространства температурный градиент в любой точке направлен перпендикулярно плоскости полупространства, коллинеарен орту оси 2, следовательно нормальный градиент температуры для боковых плоскостей схемы будет равен нулю, что соответствует условию отсутствия теплообмена через эти поверхности и необходимости его учёта. Плоскость полупространства - плоскость симметрии температурного поля, что также означает отсутствие теплообмена и нулевой градиент. Плоскость, параллельная границе полупространства, расположена на таком удалении, 0,28м, которое позволяет за общее время 80с, узлам, в ней

z

z

Z

лежащим, прогреться до температуры, не превышающей 1% от максимальной. Наглядные результаты расчёта для Г = 80с представлены на рис.3. Сопоставление изменения температуры по глубине 2 для моментов времени 5, 20, 40, 80 с, расположение соответствующих кривых снизу вверх, приведено на рис.4, а изменения температуры во времени для узла на плоскости полупространства - на рис.5.

Рис. 3. Изолинии температур тестовой задачи при Г = 80с

—

1

1

1 Гь

г!

1 \

1 \

1 \

1 Ч

*

1 г

1 \

1 \

1 \

|1 у

|

¡С

V

л\ — V. Ын-

0 0 — 01 0 —— 02 0 —1 03 0 04 0 05 0 Об 0 07 0 08 0 9 1 0 1 0 12 0 13 0 —04 0 5 0 б 0 17 0.18

Рис. 4. Изменение температуры по глубине для Г = 5;20;40;80с (соответствующие кривые расположены снизу вверх)

т,с г

—

х"

¡Г

■г.

Объект исследования

Для практического применения программы была выбрана задача тепловой нагруженности башмака с колодками дискового тормоза скоростного вагона ввиду его малой освещённости в литературных источниках и актуальности развития высокоскоростных пассажирских перевозок.

Расчётная схема башмака с колодками является сборочной единицей с контактными поверхностями, обладающими термическим сопротивлением (рис.6). Применение тонкого слоя контактных элементов, с коэффициентом теплопроводности, учитывающим давление и шероховатость на контактных поверхностях [3] позволяет вести расчёт сборки как единого тела описанным выше способом. К схеме возможно присоединение и тормозного диска, если решён вопрос взаимодействия конечноэлементных сеток различной топологии, как у башмака, колодок и диска. В качестве альтернативы, можно позаимствовать подход из хорошо разработанной методики расчёта колодочного торможения, где введён коэффициент распределения суммарной тепловой мощности торможения между колодкой и колесом. Сама мощность торможения находится для узла как произведение скорости проскальзывания, значения нормальной силы и коэффициента трения.

Рис. 6. Распределение контактных давлений между колодкой и диском с учётом касательных сил при суммарной силе нажатия на колодку 23кН

Вследствие нахождения поверхности трения в постоянном контакте, примем допущение, что теплообмен с диском отражён в коэффициенте перераспределения теплового потока. Все остальные поверхности обдуваются воздухом со скоростью, равной скорости поезда, отводят теплоту за счёт конвекции, граничные условия по конвективному теплоотводу можно позаимствовать из работы [2], и лучеиспусканием с введением коэффициента черноты тела.

Числовые значения параметров брались следующие (размерности в системе СИ): 31% от мощности торможения приходится на колодку с башмаком, коэффициент трения между колодкой и диском 0,25 , торможение со скорости 55,56 до полной остановки за промежуток времени 80 , шаг интегрирования по времени 5 , коэффициент черноты тела 0,9 , коэффициент загрязнения 0,9 ; коэффициент теплопроводности колодки - 50,3 , башмака - 50,6, удельная теплоёмкость материала

Рис. 5. Изменение температуры точки на поверхности подвода теплоты от времени

колодки 363,6, башмака - 465, плотность

материала колодки - 5713 , башмака - 7660 . Для слоя контактных элементов коэффициент теплопроводности, вычисленный на основе среднего значения контактного давления, принятого 10МПа , получен равным 6,217, остальные теплофизические характеристики как и у материала башмака. Распределение контактных давлений на поверхности трения от нормальной силы нажатия на колодку 23кН с учётом касательных сил показано на рис.6. Конечно-элементная схема содержит 18789 узла и 14068 конечных элемента.

Изменение температуры на поверхности трения по времени показано на рис.7 на примере одного из наиболее нагреваемых узлов.

На 45-тои секунде температура узла максимальна. Изолинии температур объекта для этого момента времени представлены на рис.8, и после прекращения торможения - на рис.9. Полученные температуры следует оценивать как завышенные, поскольку эпюра контактных сил для расчта бралась постоянной во времени. На самом же деле на первых секундах торможения происходит интенсивный износ материала в области пиковых контактных сил и эпюра становится более равномерной, как и само распределение температур. Также не учитывалось изменение свойств материала в зависимости от температуры в связи с отсутствием таковой информации для материала колодки, которая испытывает наибольшее нагревание.

Выводы

Описанный способ организации глобальной матрицы в сочетании с прямым методом Холецкого показал хорошие результаты при решении системы линейных уравнений относительно температуры. Дальнейшее его применение для нахождения перемещений, возникающих в результате температурного расширения материала, показало неэффективность метода, связанную не столько с большей размерностью самой матрицы, сколько с большим в три раза количеством ненулевых элементов в строке по причине взаимозависимости деформаций. Возможно более эффективным в этом случае будет применение итерационных методов.

Литература

1. Новиков С.П. Напряжённо-деформированное состояние в области контакта массивных деталей и оболочек. Дис. канд.техн.наук. Брянск: БГТУ, 2003. -197 с.

Брянский государственный технический университет

/

/

/

/

Рис. 7. Изменение температуры в точке с пиковым значением контактной силы

Рис. 8. Изолинии температур при t = 45c

Рис. 9. Изолинии температур в момент полной остановки. Внизу показано поперечное сечение по левому максимуму температур

2. Титарёв Д.В. Обоснование и разработка рациональной конструкции диска тормоза пассажирского вагона. Дис.канд.техн.наук. Брянск, 2008.-115 с.

3 Сакало В.И., Коссов В.С. Контактные задачи железнодорожного транспорта. - М.: Машиностроение, 2004. - 496с.

SIMULATION TIME-DEPENDENT THERMAL FIELDS BY USING CO-PRECIPITATION FINITE ELEMENT MESHES A.A. Mishin

The programm, which simulation time-dependent thermal process for solids, which contact nodes in finite elements meshes is pairwise interaction, is developed. The method of special organized matrices for solve of a system of linear equation by Holessky's method is developed too

Key words: thermal model testing, general matrices