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

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

(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)

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

Последнее выражение перепишем в виде

(25)

Лекция 13

Тема 3.1 Термонапряжения при индукционном нагреве. Постановка задачи. Алгоритм решения МКЭ

Среди различных численных методов решения механических задач (МКГ, МКР) в данном случае выбор делается в пользу метода конечных элементов, поскольку он выгодно отличается от остальных.

Расчет электродинамических усилий, перемещений и концентрации напряжений в элементах конструкций сводится к определению компонентов векторов перемещений точек тела

, (1)

деформаций

(2)

и напряжений

, (3)

где символ «Т» означает операцию транспонирования матриц. Все эти неизвестные являются функциями координат точек тела.

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

В статической задаче компоненты вектора напряжений должны удовлетворять уравнениям равновесия

(4)

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

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

, (5)

где ,,– компоненты вектора внешних поверхностных нагрузок; l, m, n - направляющие косинусы единичной нормали к поверхности тела.

Компоненты вектора перемещений однозначно связаны с компонентами вектора деформаций соотношениями Коши:

, , , , , . (6)

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

(7)

где – известные функции координат. Заметим, что поверхности и должны образовывать полную поверхность (S) тела, то есть .

Замкнутая система уравнений краевой задачи получается из уравнений (4)-(7), дополненных физическими уравнениями, связывающими векторы напряжений и деформаций . Последние строятся на основе физических и математических моделей конструкционных материалов.

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

, (8)

где первое слагаемое представляет собой вариацию потенциальной энергии деформации, второе – работу внешних поверхностных и объемных сил, то есть

. (9)

Выполнение принципа возможных перемещений равносильно выполнению дифференциальных уравнений равновесия (4), а также краевых условий (5) и (7). Уравнения сплошности в форме (6) предполагаются выполненными. Остается дополнить уравнение (8) физическими уравнениями.

Непосредственно из принципа возможных перемещений можно получить вариационный принцип Лагранжа в виде

, (10)

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

В осесимметричной задаче данной теории рассматривается тело вращения (рис. 1), внешние нагрузки на котором (а также температура) симметричны относительно его оси.

Рис..1 Тело вращения с внешними нагрузками, симметричными относительно его оси

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

, (11)

, (12)

где верхний индекс «Т» означает операцию транспонирования.

Соотношения Коши записываются в форме

, (13)

где u, w – компоненты перемещений, r – текущий радиус.

Для изотропного материала уравнения упругости принимаются в виде

, (14)

где матрица упругости

; (15)

вектор дополнительных деформаций

. (16)

Если эти деформации температурные, то

. (17)

В соотношениях (15)-(17): Е, - упругие постоянные материала; - коэффициент линейного температурного расширения; - перепад температур.

Для некоторого осесимметричного конечного элемента с вершинами i, j, m (рисунок) вектор искомых узловых перемещений имеет следующую структуру

, (18)

а перемещения точек внутри элемента представляются в виде:

, (19)

где матрица функций формы элементов:

и т.д., (20)

где – площадь элемента, а коэффициенты , , и другие ()определяются с помощью зависимостей типа

. (21)

Вектор деформаций выражается через вектор узловых перемещений с помощью зависимости

, (22)

в которой матрица градиентов , как и в плоской задаче, состоит из трех блоков

, (23)

каждый из которых имеет структуру типа

, . (24)

Однако, в отличие от плоской задачи, здесь зависит от координаты r, следовательно, деформации внутри конечного элемента в общем случае не будут постоянными.

Теперь с помощью соотношения (23) выразим напряжения (14) в конечном элементе через узловые перемещения.

Получим

. (25)

Система разрешающих уравнений МКЭ для осесимметричной задачи имеет тот же вид, что и для объемной, то есть

, (26)

где матрица жесткости в конструкции в целом описывается как

; (27)

- число элементов; - объем элемента.

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

, (28)

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

Эти векторы, отнесенные к конечным элементам (n = 1, 2,…, N_Э), находятся из следующих соотношений:

, (29)

, (30)

, (31)

где матрицы распределенных объемных и поверхностных нагрузок имеют соответственно следующую структуру:

; . (32)

В осесимметричной задаче, как и в плоской, матрица жесткости конечного элемента имеет вид

. (33)

Однако матрица градиентов в данном случае зависит от координат, что осложняет интегрирование уравнения (33). В практических расчетах рассматриваемый интеграл можно вычислить приближенно, определив матрицы и для центра тяжести конечного элемента с координатами

, . (34)

С учетом этого, а также соотношения , из (30) найдем

. (35)

Аналогичным образом могут быть найдены и приближенные значения векторов узловых сил (26)-(28). Опыт показывает, что при достаточно мелкой сетке конечных элементов рассмотренный прием обеспечивает приемлемую точность вычислений. Отметим, что точный расчет интегралов (26)-(30) уравнений может быть выполнен с помощью L -координат.

Лекция 14

Тема 4.1 Электромагнитная задача. Постановка проблемы. Алгоритм решения МКЭ

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

(1)

(2)

(3)

(4)

Здесь { H}, { B }, { D } – векторы напряженности магнитного поля, магнитной и электрической индукции, – вектор плотности приложенного тока, – вектор плотности индуцированного тока,

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

Поле принимается квазистационарным. Под этим понимается отсутствие запаздывания электромагнитной волны в воздухе (но не в металле). Это допущение позволяет пренебречь токами смещения по сравнению с токами в проводниках.

Не учитываются потери на гистерезис при нагреве ферромагнитных тел в силу их незначительности по сравнению с потерями от вихревых токов.

Принятые допущения позволяют упростить решение рассматриваемой задачи.

Граница раздела магнитных сред описывается системой:

(7)

Последнее выражение учитывает скачок вектора на границе раздела сред.

При тангенциальные составляющие напряженности на границе раздела непрерывны

(8)

Кроме этого необходимо задать:

- уравнения поверхностей, отделяющих друг от друга среды i и j, f_ij(x,y,z)=0;

- начальные величины E₀(x,y,z), H₀(x,y,z) в момент времени t ₀ в произвольной точке исследуемого объема с границей S;

- касательные составляющие вектора или в произвольной точке поверхности в произвольном временном интервале от t₀ до t, или распределения полей и вне исследуемого объема V;

- функциональные зависимости параметров ε, μ, γ от координат пространства или от напряженности соответствующего поля.

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

, (12)

Решение задачи электромагнитного поля достигается использованием векторного магнитного потенциала { A } и скалярного электрического потенциала V, которые выражаются следующим образом:

(14)

(15)

Чтобы функция была определена, нужно определить значение ее дивергенции. Для этого добавляется условие, которое называется калибровкой Кулона

(16)

В результате получим следующую систему уравнений

(17)

(18)

(19)

Используя соотношение

(20)

при μ =const из (17) получим уравнение

(21)

Уравнение Пуассона (21) дополняется граничными условиями Дирихле и Неймана на различных участках границы:

на S₁ (22)

на S₂ (23)

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

Решение краевой задачи расчета магнитного поля в изотропной среде (21)..(23) эквивалентно минимизации энергетического функционала:

(24)

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

(25)

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

В результате векторный потенциал внутри m – го элемента треугольника определяется значениями потенциала в вершинах треугольника, то есть является линейной функцией координат x и y.

(26)

где: - постоянные коэффициенты функций формы N_i, вычисляемые в зависимости от пространственных координат узлов элемента m; - комплексные амплитуды вектора в узлах конечного элемента.

В дискретной модели функционал (24) определяется суммой вкладов всех КЭ, входящих в ансамбль

, (27)

а условие его минимума приобретает вид

(28)

где Ne – полное число всех элементов.

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

Матрица жесткости определяется следующим выражением:

(29)

Матрица вихревых токов рассчитывается следующим образом

. (30)

Матрица внешних источников тока вычисляется согласно выражению

(31)

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

Согласно выражению (28) элементные матрицы (29) должны объединяться в глобальные матрицы, характеризующие поведение дискретной системы в целом.

(32)

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

(33)

Решение данной задачи осуществляется итерационным методом. Краевые условия вида Дирихле учитываются путем принудительного исключения столбцов и строк глобальных матриц (32), относящихся к узлам дискретной системы, лежащих на удаленных границах S области V. Условия симметрии удовлетворяются при ансамблировании элементов автоматически. Распределенные параметры магнитного поля вычисляются по выражению (18):

; (34)

(35)

. (36)

Напряженность электрического поля

. (37)

Мощность внутренних источников тепла, характеризующих нагрев проводящих тел индукционной системы, вычисляется для каждого элемента по закону Джоуля-Ленца

, (38)

где - величина, сопряженная к .

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

Магнитные силы в токопроводящем проводнике определяются численным интегрированием

(39)

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

(40)

(41)

(42)

(43)

(44)

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