10606

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

Лекция

Математика и математический анализ

Решение МКЭ тепловой задачи для цилиндра. Алгоритм расчета Математическая модель линейной задачи теплопроводности с внутренним тепловыделением в цилиндрических координатах имеет вид: 1 с граничными условиями:

Русский

2013-03-29

635.5 KB

33 чел.

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

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

                                               (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 – номер шага.

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

                                                                                                 (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). Кроме того, обоснованны следующие допущения:

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

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

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

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