Применение треугольных конечных элементов в задаче теплопроводности

Постановка задачи. Стационарное уравнение теплопроводности имеет вид

(7)

где u – температура, – коэффициенты теплопроводности (Вт/(м град)), Q – плотность тепловых источников внутри тела (Вт/ ), эта величина считается положительной, если тепло подводится к телу.

Граничные условия могут быть трех типов:

1. На границе тела (или ее части) задана температура, которая может быть функцией координат

.

(8)

2. На границе задан конвективный теплообмен. Конвективный теплообмен – процесс переноса тепла, происходящий в движущихся текучих средах (жидкостях либо газах) (теплообмен с окружающей средой). Конвективный. теплообмен характеризуется величиной , – ккоэффициент теплообмена Вт/ ( град), – температура окружающей среды.

3. На границе задан тепловой поток q (то есть извне подводится или отводится тепло). Тепловой поток q положителен, если тепло подводится.

Граничные условия п. 2 и З записываются с помощью смешанного условия

(9)

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

Решение уравнения (7) с граничными условиями (8) – (9) эквивалентно отысканию минимума функционала:

Разобьем область V на конечные элементы и аппроксимируем решение u конечной суммой где n — число узлов, — функции формы. Записывая условие минимума функционала , получим СЛАУ относительно узловых значений:

Или

,

Элементная матрица системы имеет вид

	(10)
	(11)

Рассмотрим треугольный элемент с тремя узлами i, j, k. Функции формы линейного треугольного элемента (где определяются соотношением (6)). Температура па элементе определяется формулой

Запишем матрицу градиентов от функций формы

матрица свойств

Теперь можно вычислить матрицу теплопроводности элемента (10) (элементную матрицу жесткости). Первое слагаемое принимает вид

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

Вычисляя произведения матриц, имеем

(12)

Второй интеграл должен быть вычислен по поверхности

(13)

Функции формы зависят от х и у, поэтому произведения вида не могут быть вынесены за знак интеграла. Кроме того, значение интеграла зависит от того, на какой поверхности наблюдается конвективный теплообмен. Если, например, конвекции подвержена сторона между узлами i и j, то равно нулю вдоль этой стороны и интеграл сводится к следующему выражению

(14)

Рис. 5: Пример сетки из треугольных элементов. Перечислены номера узлов и номера элементов

Таблица 1: Списки узлов из которых состоит элемент (массив elements)

Если любая из двух других сторон подвержена конвекции, то расположение отличных от нуля членов в (13) будет иным, чем в (14).

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

Программная реализация. Используемые данные: n ­­­-­ число узлов, m - число элементов; сетка (координаты узлов) mesh [ n ] [2]: mesh [ i ] [1] — ко-ордината , mesh [ i ][2] — координата , i, j = ; массив elements [ m ][3] = (nodе1, node2, node З) — ­­­ список элементов (списки узлов из которых состоит элемент), для сетки на рисунке 5 массив elements имеет вид: elements [1, 1] = 1, elements [1, 2] = 5, elements [1, 3] = 4 (см. Табл. 8); элементная матрица размера ; элементный вектор правых частей , ; вектор правых частей ; -вектор решения, граничное условие Дирихле, свойства материала — массив material [ m, 2] (первый индекс — номер элемента, второй индекс — коэффициент теплопроводности ), так material [1,2] есть значение первого элемента, material [20, 1] есть значение двадцатого элемента. Толщина элемента предполагается единичной.

Поиск по сайту: