Если у вас возникли сложности с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой - мы готовы помочь.
Предоплата всего
Подписываем
Тема 2.4 Решение МКЭ тепловой задачи для цилиндра. Алгоритм расчета
Математическая модель линейной задачи теплопроводности с внутренним тепловыделением в цилиндрических координатах имеет вид:
(1)
с граничными условиями:
(2)
(3)
(4)
, (5)
где r и z – радиальная и аксиальная координаты; W – функция распределения внутренних источников тепла, полученных в результате решения электромагнитной задачи (1)..(4); – коэффициент температуропроводности,; - степень черноты материала загрузки; – коэффициент излучения абсолютно черного тела; - коэффициент теплообмена с окружающей средой конвекцией и зависит от геометрических размеров и формы стенки нагреваемого изделия; q – тепловой поток от корпуса взрывателя к корпусу посадочного гнезда.
Начальные условия характеризуются произвольным в общем случае пространственным распределением
(6)
В граничных условиях отражены три вида теплообмена: конвективный, передача тепла теплопроводностью и излучением. Это связано с технологическим процессом нагрева взрывателя. На первом этапе осуществляется индукционный нагрев холодной загрузки при наличии теплообмена конвекцией с окружающей средой. В это время учитывается отток тепла от взрывателя к посадочному гнезду в виде плотности теплового потока.
Как было сказано выше, решение тепловой задачи проведем методом конечных элементов, который дает возможность достаточно точно учитывать все нелинейности, путем изменения всех нелинейных величин с каждым шагом по времени, а также задать сложную геометрию нагреваемого изделия. Следуя МКЭ, дифференциальному уравнению (45) ставится в соответствие вариационная формулировка о минимизации энергетического функционала, характеризующего тепловое состояние тела :
(7)
где Lh – граница с конвективным теплообменом; Lq – граница, которую пронизывает поток 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) с шириной полуленты mk. Величина S равна S=Nu, Nu – количество узлов. Учитывая особенности этих матриц, в памяти ЭВМ достаточно хранить коэффициентов для каждой матрицы, что существенно снижает потребные ресурсы ЭВМ по памяти и позволяет решать задачи с густой сеткой КЭ. Практически матрицы ансамбля хранятся в виде одномерных массивов размерностью N, а работа с ними производится с помощью вычисляемых индексов.
Полученные матрицы , и с учетом замены временной производной конечно-разностным аналогом, объединяем в систему уравнений (схема Галеркина).
(24)
где t – временной шаг, n – номер шага.
Среди различных численных методов решения механических задач (МКГ, МКР) в данном случае выбор делается в пользу метода конечных элементов, поскольку он выгодно отличается от остальных.
Расчет электродинамических усилий, перемещений и концентрации напряжений в элементах конструкций сводится к определению компонентов векторов перемещений точек тела
, (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-координат.
В общем случае процесс непрерывного индукционного нагрева описывается системой уравнений Максвелла для электромагнитного с соответствующими краевыми условиями.
(1)
(2)
(3)
(4)
Здесь {H}, {B}, {D} – векторы напряженности магнитного поля, магнитной и электрической индукции, – вектор плотности приложенного тока, – вектор плотности индуцированного тока,
Вышеприведенные уравнения описывают связь между различными средами, входящими в систему. Индукционный нагрев на средних частотах характеризуется отсутствием свободных зарядов в системе рассматриваемых сред, поэтому из системы уравнений (1)..(4) можно исключить уравнение (4). Кроме того, обоснованны следующие допущения:
Принятые допущения позволяют упростить решение рассматриваемой задачи.
Граница раздела магнитных сред описывается системой:
(7)
Последнее выражение учитывает скачок вектора на границе раздела сред.
При тангенциальные составляющие напряженности на границе раздела непрерывны
(8)
Кроме этого необходимо задать:
- уравнения поверхностей, отделяющих друг от друга среды i и j, fij(x,y,z)=0;
- начальные величины E0(x,y,z), H0(x,y,z) в момент времени t0 в произвольной точке исследуемого объема с границей S;
- касательные составляющие вектора или в произвольной точке поверхности в произвольном временном интервале от t0 до t, или распределения полей и вне исследуемого объема V;
- функциональные зависимости параметров ε, μ, γ от координат пространства или от напряженности соответствующего поля.
При индукционном нагреве на средних частотах влиянием электрической индукции можно пренебречь. Отсутствие в рассматриваемой системе движущихся постоянных магнитов также исключает появление дополнительных источников внутри проводящих материалов. Тогда связь между напряженностью электрического поля и плотностью токов будет иметь вид
, (12)
Решение задачи электромагнитного поля достигается использованием векторного магнитного потенциала {A} и скалярного электрического потенциала V, которые выражаются следующим образом:
(14)
(15)
Чтобы функция была определена, нужно определить значение ее дивергенции. Для этого добавляется условие, которое называется калибровкой Кулона
(16)
В результате получим следующую систему уравнений
(17)
(18)
(19)
Используя соотношение
(20)
при μ=const из (17) получим уравнение
(21)
Уравнение Пуассона (21) дополняется граничными условиями Дирихле и Неймана на различных участках границы:
на S1 (22)
на S2 (23)
Такое упрощение условий задачи объясняется тем, что дальнейший переход к конечно-элементной формулировке намного облегчается для линейной задачи. Реальные нелинейные задачи решаются на базе линейных моделей с помощью итерационных алгоритмов расчета.
Решение краевой задачи расчета магнитного поля в изотропной среде (21)..(23) эквивалентно минимизации энергетического функционала:
(24)
Сущность метода, основанного на МКЭ, заключается в исследовании глобальной функции процесса, в данном случае векторного потенциала , в дискретных частях анализируемой области V, которая должна быть предварительно разбита на конечные смежные подобласти (конечные элементы), что позволяет свести задачу с бесконечным числом степеней свободы к задаче, содержащей конечное число параметров. При этом внутри подобластей искомая функция интерполируется степенными полиномами, сшивается на границах контакта элементов, и при условии малости геометрических размеров последних (число элементов стремится к бесконечности), оказывается решением уравнений в частных производных типа (21)..(23). В качестве интерполирующих полиномов конечных элементов треугольного вида использованы линейные функции формы вида:
(25)
Треугольные элементы позволяют наиболее просто аппроксимировать сложные геометрические границы тел. В настоящее время разработаны другие виды конечных элементов, например четырехугольные с криволинейными сторонами, что обеспечивает при сравнительно небольшом числе элементов гладкую аппроксимацию контуров области. Такая же ситуация имеет место и для объемных областей. Дальнейшие рассуждения проведем для более простой плоской задачи.
В результате векторный потенциал внутри m – го элемента треугольника определяется значениями потенциала в вершинах треугольника, то есть является линейной функцией координат x и y.
(26)
где: - постоянные коэффициенты функций формы Ni , вычисляемые в зависимости от пространственных координат узлов элемента 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 по всей области ферромагнитных макроэлементов, затем вычисляются распределенные параметры поля, что позволяет на следующей стадии расчета корректировать μ внутри каждого конечного элемента в зависимости от значения напряженности магнитного поля в данной области. Итерации повторяются до полной сходимости процесса. Определение магнитной проницаемости производится с помощью введения в программу расчета полинома, аппроксимирующего кривую намагничивания.
PAGE 49