Шабров Н.Н. - Метод конечных элементов в расчётах деталей тепловых двигателей (1061803), страница 16
Текст из файла (страница 16)
32) где !!!'„— трехмерные функции формы; я =- 1, 2, ..., !т!; Л' двухмерные функции формы; Т = 1, 2„..., 4 или у = 1, 2, ..., 8. Для задания функции Я в объеме конечного элемента и функций 6, Т , д на поверхности элемента можно воспользоваться соотношениями типа (5.28) или (4.116). Замечание, высказанное относительно порядков интегрирования в (5.12) — (5.14), справедливо и в данном случае для выражений (5.31) — (5.33). Глава 6 НЕСТАЦИОНАРНЫЕ ЗАДАЧИ ТЕОРИИ ТЕПЛОПРОВОДНОСТИ 6.!. УРАВНЕНИЯ ТЕПЛОВОГО БАЛАНСА МКЭ В НЕСТАЦИОНАРНОЙ ЗАДАЧЕ ТЕОРИИ ТЕПЛОПРОВОДНОСТИ В прикладных исследованиях очень часто требуется решать задачи нестационарной теплопроводности.
Например, анализ напряженного состояния деталей тепловых двигателей на переходных режимах требует предварительного решения задачи нестационарной теплопроводности. Уравнение для линейной трехмерной задачи нестационарной теплопроводности имеет вид где с — удельная теплоемкость материала; р — плотность материала. Решение уравнения (6.1) должно удовлетворять краевым условиям (3.2) — (3.4) и начальному условию в виде Т(х, д, г, 1) ~,, = Т,(х, !!, г), (6.2) где Т0 (х, у, г) — температура в начальный момент времени. 105 Для того чтобы задачу нестационарной теплопроводности можно было решать методом конечных элементов, нужно для уравнения (6.1) построить соответствующий функционал. В работах 18, 191 описывается один из простых приемов получения такого функционала, когда член с частной производной по времени в уравнении (6.1) рассматривается в каждый фиксированный момент времени как функция источников или стоков теплоты.
При таком подходе функционал для уравнения (6.1) будет иметь вид (3.5) при условии, что функция Я в (3.5) заменяется разностью д( дТ (6.3) Рассмотрим отдельно интеграл в (3.5), соответствующий функции стока или источников теплоты Я. С учетом (6.3) получим у() —— — 1(Я вЂ” ср д ) Тдш (6.4) )(.~ ) = — Я вЂ” ср — Т~'~ (1о. у(Е) (6.5) Далее на основании представления (3.10) можно написать дТ(') дТ 1~1(е) (е) ~') (6.6) д( где 1) =- д д1 Тогда с учетом (6.6) и (3.10) интеграл (6.5) будет иметь следующий вид: у,'~) = 1.( авч ) Х" срМ" (1()а" Ь вЂ” 1) а" ~ Я" Я(Ь.
(6.7) г(е) г(с) Введем обозначение С') = ~ Х~' М" ср((а, ) (е) (6.8) где С(') — матрица теплоемкости конечного элемента. При этом выражение для у~ запишется как (е) (6.9) где Г~~' дается выражением (3.16). Условие стационарности функционала (3.8) в форме (3.19) с учетом представления (6.9) приводит к следующему матричному линейному обыкновенному дифференциальному уравнению первого порядка С1.1+ Р1.1 = Г. (6.10) 106 Вклад каждого конечного элемента в (6.4) определится соот- ношением Здесь С вЂ” глобальная матрица теплоемкости всей системы конечных элементов, для которой имеет место равенство С = ~ а" СЧ"1а". е (6.
11) Матрица теплопровадности Р и вектор тепловых сил Г в (6.10) даются соотношениями (3.21) и (3.22). Из (6.8) видна, что выражение для матрицы теплоемкасти конечного элемента совпадает с точностью да скалярного множителя с выражением для матрицы поглощения (3.15). Отсюда следует, что при вычислении матрицы теплоемкости конечного элемента можно целиком воспользоваться результатами, получе1шыми при выводе выражений для матрицы поглощения соответствующего конечного элемента.
Подставим (6.13), (6.14) и (6.15) в уравнение (6.10) и получим С(1)1 — 1!1)+Р(Л';111+ Л1;1!1) =- Л';Г;+М Г;. (6.16) 107 6Л. ПРИМЕНЕНИЕ МЕТОДА КОНЕЧНЫХ РАЗНОСТЕЙ ДЛЯ РЕШЕНИЯ ДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ НЕСТАЦИОНАРНОЙ ТЕПЛОПРОВОДНОСТИ Матричное линейное дифференциальное уравнение (6.10) относительно неизвестного вектора 13 является одномерным уравнением по временной координате 1. Оказывается удобным для представления приближенного решения уравнения (6.10) воспользоваться одномерной схемой метода конечных элементов по временной координате.
Для этого разобьем отрезок временной аси, на котором ищется решение уравнения (6.10), последовательностью интервалов Лт, вообще говоря, различной величины. На каждом интервале определим линейный одномерный конечный элемент по времени, который имеет следующие функции формы: А~1 (г) = Рт — т)IР«); А'~ (т) = «$(Ю (О ~ т ( ~т) (б 12) Тогда интерполяционное соотношение для вектора 13 на интервале Лт запишется 13 (т) — Л'; 13; + У;Ц, (6.13) где 111 и 1!; — значения решения уравнений (6.10) в точках интерполяции или узлах конечного элемента. Если закон изменения по времени вектора нагрузки в (6.10) близок к линейному, то справедлива интерпаляциопное соотношение Г(т) = Г;У;+ Г,У;, (6.14) где Г; н Г; — значения вектора Г в узлах конечного элемента.
Дифференцирование в (6.13) па времени даст —,', и = 1) = —,', (11, — Ц). (6.15) Уравнение (6.16) является исходным для получения различных конечно-разностных схем. Рассмотрим три из них. 1. Центральная разностная схема первого порядка получается из (6.16) при т = О,БЛт и имеет вид 2. Левая разностная схема первого порядка получается из (6.16) при т = 0 и имеет вид (6. 18) 3. Правая разностная схема первого порядка получается из (6.16) при т =- Лт и имеет вид (6.19) Аи'~ =- Ви" > + а ! (6.20) где А и  — матрицы; б — вектор, конкретный вид которых завися от используемой конечно-разностной схемы, и~(0) =- То. 6.3.
ПРИМЕНЕНИЕ МЕТОДА ГАЛЕРКИНА ДЛЯ РЕШЕНИЯ ДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ НЕСТАЦИОНАРНОЙ ТЕПЛОПРОВОДНОСТИ Весьма эффективным способом решения нестационарных задач является сочетание одномерной схемы метода конечных элементов и метода Галеркина. Метод Галеркина принадлежит классу методов взвешенных невязок 113] и представляет собой альтернативный подход при построении систем уравнений МКЭ. Главное преимущество такого подхода заключается в том, что для получения системы уравнений МКЭ необязательно существование !08 Сумма уравнений (6.18) и (6.19) приводит к уравнению (6.17).
Рассмотрим любую из приведенных конечно-разностных схем для первого конечного элемента по времени. В этом случае значение неизвестного в 1-м узле 11; совпадает с заданным начальным значением и в этом смысле является известным. Тогда значение неизвестного в 1-и узле Г1; определяется в результате решения системы линейных алгебраических уравнений, Для второго конечного элемента значение неизвестного в 1-м узле Ц совпадает с найденным значением 13; для первого конечного элемента и в этом смысле также является известным. Таким образом, решение нестационарной задачи теплопроводности сводится к решению рекуррентной последовательности стационарных задач на каждом временном слое (6.21) Рис.
б.1. Сравнение решений дифференциального уравнения нестаципнарнай теплопроводнасти: — — метод конечных разностей; — — — — метод Галеркнна Рассмотрим уравнение (6.21) в сочетании с одномерной схемой МКЭ. Тогда соотношение (6.21) запишется так: Ьт ~ (СЬ + Р11 — Г) 611 с(т —.
О, (6.22) о где на основании (6,13) имеет место выражение 611 = Ж,61.1; + Лг~б Ц. (6.23) Вариация ЬЬ; всегда равна нулю, поскольку для всех конечных элементов значение 1.1; является заданным. !09 вариационного принципа. В связи с этим появляется возможность применять МКЭ для решения задач, не имеющих соответствующего функционала.
Применяя метод Галеркина (13] непосредственно к матричному дифференциальному уравнению (6.10), получим ) (С11+ РЬ вЂ” Г) И) 4т = О. о Соотношение (6.21) выражает условие ортогональности невязки уравнения (6.10) и произвольной вариации решения 1.1. Интегрирование в (6.21) осуществляется по соответствующему отрезку времени. Подставляя в (6.22) выражения (6.23) и (6.13) — (6.15), получим ~т ~ Х~ ~ — С Щ вЂ” 13;) + Р(У;Ц + Л~;131) — (Лг;Г;+ У1Г1)] 613! дт=О. (6.24) Интегрирование в (6.24) с учетом зависимостей (6.12) при условии произвольности вариации 6].1, приводит к следующей системе алгебраических уравнений: (л', +з') '=(л', — з') +з('-')(") Видно„что структура системы уравнений (6.25) представляет рекуррентное матричное соотношение типа (6.20). Следует отметить, что решение задачи нестациопарной теплопроводности с использованием схемы (6.25) дает лучшие результаты для первых временных шагов, чем решение с использованием центральной разностной схемы (6.17).
Сказанное иллюстрируется на примере решения задачи о разогреве призматического стержня (рис. 6.1). Глава 7 ОРГАНИЗАЦИЯ ДАННЫХ ПРИ РЕШЕНИИ ЗАДАЧ МНЭ 7.1. АВТОМАТИЧЕСКАЯ ГЕНЕРАЦИЯ ГЛОБАЛЬНЫХ КООРДИНАТ И НОМЕРОВ УЗЛОВЫХ ТОЧЕК КОНЕЧНЫХ ЭЛЕМЕНТОВ Затрата физических и творческих усилий для разработки большой мощной программы со сложной структурой логических связей разнообразных программных модулей оправдана в том случае, если эта программа предназначена для многократного использования в серийных расчетах какого-либо класса задач. Решение прикладных задач теории теплопроводпости и термоупругости на основе МКЭ, имеющих практическое значение, требует рассмотрения систем с большим числом конечных элементов.
Это обстоятельство вызывает необходимость формирования и обработки огромных числовых массивов данных, которые, разумеется, невозможно одновременно разместить в оперативной памяти ЭВМ. В связи с этим наряду с точностью и высокой скоростью счета программа при ее использовании в серийных расчетах должна удовлетворять некоторым сервисным функциям, среди которых важное место отводится автоматической генерации массивов глобальных координат и номеров узловых точек конечных элементов.
Очевидно, что автоматическая генерация координат и номеров узловых точек значительно снижает возможность появления оши- 110 бок по сравнению с ручной подготовкой. Особенно это относится к построению трехмерной сетки конечных элементов. В современной литературе рассматриваются различные под-- ходы к решению проблемы автоматической генерации сетки конечных элементов. Один из них заключается в использовании техники изопараметрических элементов и обсуждается в работе (19 ). Проиллюстрируем процесс автоматической генерации сетки конечных элементов на примере плоской двухмерной задачи. 1 7 1,1 1~ 1б Ю 6) Х 1 7 '5 Ф Х 11 Ч ~то 11 Рис. 7.1.
Вид макроэлементов для формирования сетки линейных (а) и квадратичных (б) четырехугольных ко- нечных элементов Вернемся к рассмотрению изопараметрического четырехугольного конечного элемента второго порядка. Форма такого элемента определяется интерполяционными соотношениями (4.34) и (4.35). Очевидно, эти соотношения можно использовать для вычисления глобальных координат точек конечного элемента по известным локальным координатам. В качестве локальных координат можно взять координаты точек равномерной сетки, построенной в локальной системе координат конечного элемента (рис.
?.1). Таким образом, вводится локальная сетка. Соотношениями (4.34) и (4.35) локальная сетка переводится в глобальную, точками которой являются узловые точки получившихся базовых конечных элементов. Определим четырехугольный изопараметрический макроэлемент второго порядка, когда наряду с интерполяционными соотношениями (4.34), (4.35) заданы числа разбиений п~ и пп на 111 базовые конечные элементы по каждой из локальных координат. Процедуры генерации сетки базовых четырехугольных конечных элементов первого и второго порядков отличаются друг от друга, Рассмотрим каждую из них в отдельности.