10606

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

Лекция

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

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

Русский

2013-03-29

635.5 KB

32 чел.

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

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

                                               (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


 

А также другие работы, которые могут Вас заинтересовать

72484. Основные сведения о конструкции и технологии заклепочных соединений, классификация, области применения 469 KB
  Заклепочное соединение является неразъемным. В большинстве случаев его применяют для соединения листов и фасонных прокатных профилей. Соединение образуют расклепыванием стержня заклепки, вставленной в отверстие деталей – рис.3.1, где 1 – обжимка; 2 – прижим при машинной клепке...
72485. Східні слов’яни. Зародження української державності. Київська Русь 78 KB
  На початку XII ст. Київська Русь розпадається на окремі князівства. Син Володимира Мономаха Мстислав продовжував справу свого батька, тримав укупі руські землі та зберігав владу Київського престолу. Але він був останнім київським правителем, якому це було під силу.
72486. Українські землі у складі Великого князівства Литовського та Речі Посполитої 100 KB
  Економічний розвиток та соціальний устрій українських земель у ХІV – першій половині ХVІ ст. Литва здобула більшу частину білоруських українських частину російських земель. Політика Литви в українських землях Литовські князі щоб забезпечити управління величезними завойованими територіями...
72487. Ожоги, отморожения, отравления. Оказание ПМП 82.5 KB
  Первая медицинская помощь при термических ожогах. Первая медицинская помощь при химических ожогах. Первая медицинская помощь при отморожениях. Первая помощь в скрытом периоде: медленное согревание в ваннес с постепенным увеличением температуры.
72488. Электротравма. Утопление. Тепловой и солнечный удар. Инородные тела. Укусы животных, насекомых, змей. Оказание ПМП 79 KB
  Первая медицинская помощь при электротравме: устранить воздействие тока на пострадавшего выключить электроустановку откинуть электропровод и т.; осмотр пострадавшего; при потере сознания немедленно проводить искусственное дыхание при остановке сердца сочетать искусственное дыхание...
72489. Сердечно – легочная реанимация 63.5 KB
  Обеспечить проходимости дыхательных путей: освободить полость рта и глотки от инородных масс кровь слизь рвотные массы зубные протезы жвачка рукой обернутой салфеткой платком предварительно повернув голову спасаемого набок.
72490. Введение лекарственных средств 79 KB
  Перед инъекцией кожу протирают стерильным шариком смоченным 70 раствором спирта левой рукой собирают ее в складку треугольной формы основанием вниз правой рукой берут шприц и придерживая иглу и поршень быстрым движением вкладывают иглу в основании треугольника под углом 45 на глубину 12 см.
72491. ПМП при острых заболеваниях брюшной полости, почечной колике 92 KB
  Возникает рвота кишечным содержимым имеющим неприятный каловый запах. Местный перитонит возникает при отграничении участка воспаления спайками петлями кишечника; локализуется вблизи источника воспаления желчный пузырь червеобразный отросток. Возникает в результате воспалительных процессов в матке или придатках.
72492. Раны. Кровотечения 101.5 KB
  Рана – механическое нарушение целостности кожи, слизистых оболочек с повреждением глубинных тканей. Основные признаки раны: боль – наиболее выражена в местах с наибольшим количеством нервных окончаний (кончиках пальцев, надкостница, плевра). кровотечение – абсолютный признак раны.