Тема 2.4 Решение МКЭ тепловой задачи для цилиндра. Алгоритм расчета

Лекция 12

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

(1)

с граничными условиями:

(2)

(3)

(4)

, (5)

где r и z – радиальная и аксиальная координаты; W – функция распределения внутренних источников тепла, полученных в результате решения электромагнитной задачи (1)..(4); – коэффициент температуропроводности,; - степень черноты материала загрузки; – коэффициент излучения абсолютно черного тела; - коэффициент теплообмена с окружающей средой конвекцией и зависит от геометрических размеров и формы стенки нагреваемого изделия; q – тепловой поток от корпуса взрывателя к корпусу посадочного гнезда.

Начальные условия характеризуются произвольным в общем случае пространственным распределением

(6)

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

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

(7)

где L_h – граница с конвективным теплообменом; L_q – граница, которую пронизывает поток q; .

Исследуемая область аппроксимируется совокупностью элементов с конечным числом узловых точек. Функционал (51) заменяем суммой отдельных вкладов элементов, определяя, таким образом, функциональные соотношения относительно узловых неизвестных.

В качестве элементов использовались симплекс-элементы, т. е. такие, для которых интерполяционный полином имеет первую степень координат.

Вершины треугольников, обозначаемые индексами в направлении против часовой стрелки, образуют локальную систему узлов.

Для произвольного элемента е_i пробная функция выбирается линейной, т. е.

; _i , (8)

где , , -- постоянные, в общем случае отличные для различных элементов. Значения этих постоянных определяются из выражений

(9а)

(10)

(11)

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

Подставляя в (10), получим

, (12)

где

(13)

является матрицей базисных функций, а

(14)

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

Определяем вклады элементов в матрицы [K] жесткости, матрицы [C] демпфирования и в вектор {F}

источников.

(15)

(16)

. (17)

Здесь

(18)

(19)

Вектор {F} источников формируется из внутренних источников тепла w, обусловленных вихревыми токами в изделии, из конвективных потерь, определяемых коэффициентом h теплообмена, и из потока q тепла через стенку. Рассмотрим теплообмен с внешней средой по двум граничащим со средой сторонами. Примем, что имеет место общий случай граничных условий.

Рис. 1. Аппроксимирующий элемент.

Для случая, представленного на рис. 1, получим:

(20)

(21)

(22)

где

, , (23)

– вычисляется по формуле (12). Следует отметить, что в выражениях (21), (22) интегралы вычислены приближенно. Это вполне допустимо, если размеры элемента намного меньше среднего радиуса. Полученные матрицы для элементов объединяются в глобальные матрицы.

Процесс ансамблирования осуществляется с помощью процедуры поэлементного объединения. Такой способ безразличен к разнородности конечных элементов, из которых собрана исследуемая система. В результате ансамблирования преобразованных элементарных матриц , формируются глобальные матрицы ансамбля КЭ , , которые являются симметричными ленточными матрицами размерностью (S*S) с шириной полуленты m_k. Величина S равна S=N^u, N^u – количество узлов. Учитывая особенности этих матриц, в памяти ЭВМ достаточно хранить коэффициентов для каждой матрицы, что существенно снижает потребные ресурсы ЭВМ по памяти и позволяет решать задачи с густой сеткой КЭ. Практически матрицы ансамбля хранятся в виде одномерных массивов размерностью N, а работа с ними производится с помощью вычисляемых индексов.

Полученные матрицы , и с учетом замены временной производной конечно-разностным аналогом, объединяем в систему уравнений (схема Галеркина).

(24)

где Dt – временной шаг, n – номер шага.