2.9 Метод конечных элементов

Метод конечных элементов [3, 12, 13] является основным расчетным методом в современных системах автоматизированного проектирования и инженерного анализа. Он основан на идее замены непрерывной функции (в физической интерпретации – температуры, давления, перемещения и т.д.) дискретной моделью, которая строится на множестве кусочно-непрерывных функций. При этом пространственное представление исследуемого объекта осуществляется с помощью набора конечных элементов, объединяемых в ансамбли. Роль конечных элементов обычно выполняют области канонической геометрии (отрезки, многоугольники, многогранники). Например, на рисунке 20 приведено представление сечения режущей твердосплавной пластины ансамблем из треугольных конечных элементов.

Рисунок 20 – Конечно-элементное представление режущей твердосплавной пластины

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

Существует множество формулировок метода конечных элементов и множество типов конечных элементов. Здесь мы рассмотрим вывод расчетных соотношений для аппроксимации исходной задачи простейшими симплекс-элементами.

На рисунке 21 изображен двумерный симплекс-элемент. Это треугольник с прямолинейными сторонами и тремя узлами, по одному в каждой вершине.

Рисунок 21 – Двумерный симплекс-элемент

В случае линейного распределения искомой функции в рассматриваемой треугольной области приближающий ее интерполяционный полином должен имеет вид:

.

(18)

Тогда, если известны значения искомой функции в узлах

,

то коэффициенты , , соотношения (18)

можно найти из решения матричного уравнения:

,

где  следующая матрица:

.

Объединив в матрицу множители при в интерполяционной формуле (34)

,

получим матрицу так называемых функций формы треугольного элемента , вычисляемую как

.

(19)

Тогда интерполяционный полином (18) для искомой функции может быть получен с помощью следующей матричной операции:

.

(20)

Допустим, что в треугольной области (см. рисунок 21) задан линейный закон распределения температуры в виде (18). Значения координат узлов и температур в этих узлах элемента приведены в таблице 5.

Таблица 5 – Данные о распределении температуры

Номер узла	X, мм	Y, мм	Ф, ^0С
1	1	2	40
2	0	0	34
3	2	0	26

Тогда матрицы и соответственно примут вид:

, .

Обратная матрица вычисляется по известным формулам линейной алгебры:

.

Находим матрицу функций формы по формуле (19):

Найдем интерполяционный полином по формуле (20):

.

Технику связывания нескольких элементов в область проиллюстрируем на примере простой пятиэлементной конфигурации, показанной на рисунке 22 (звездочкой помечен угол элемента, с которого начинается обход узлов против часовой стрелки, в скобках указан номер элемента).

Рисунок 22  Область для иллюстрации алгоритма связывания

Приписывая узлам глобальные номера, составим массив, характеризующий всю пятиугольную область, например, так:

.

(21)

Тогда связывание реализуется в простом цикле с параметром (номер элемента), который представим с помощью следующего псевдокода:

Для стационарных режимов решение двумерного дифференциального уравнения теплопроводности с граничными условиями первого-третьего рода эквивалентно отыскание минимума функционала [5]:

,

где – функция, определяющая поле температур в двумерной области, К;

– внутренний тепловой источник или сток, Вт;

и – коэффициенты теплопроводности в направлении осей координат, Вт/(мК);

– тепловой поток заданной интенсивности, Вт/м²;

– коэффициент конвективного теплообмена, Вт/(м²К);

– температура окружающей среды, К;

– объем тела, м³;

– площадь поверхности, м²;

и – текущие координаты.

При численной реализации метода рассматривают функционал для отдельно взятого конечного элемента и заменяют искомую функцию интерполяционным полиномом типа (18). Кроме того, можно расширить данные элемента. В частности, в массив (21) можно включить дополнительные столбцы, соответствующие: значениям коэффициентов теплопроводности, условиям теплообмена на границах элемента и в направлении нормали к лицевой поверхности, значениям мощности теплового потока в узлах. При этом непрерывность функции во всей области обеспечивается равенством интерполирующих функций на границе между элементами. После указанной модификации внутри цикла по элементам вычисляются функционалы для каждого конечного элемента, находится суммарный функционал:

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

,

решение которой приводит к нахождению неизвестных значений искомой функции (температуры) в узлах.