15.1 Общий ход решения задачи на основе метода конечных элементов

Решение задач анализа при функциональном моделировании сводится к определению вектора V, который удовлетворяет исходным уравнениям задачи и граничным условиям. Часто искомую функцию задают в виде ряда

(1)

где – базисные функции от независимых переменных,

– неизвестные постоянные коэффициенты.

Методы определения постоянных коэффициентов весьма разнообразны. Так, в методе Ритца их находят из условия минимума потенциальной энергии системы, в методе коллокаций – из условия удовлетворения дифференциальных уравнений задачи в отдельных точках и т. д. Выбор базисных функций в конечном счете определяет успех решения задачи: если эти функции подобраны удачно, решение получается простым, в противном случае приходится удерживать большое число членов ряда (9.1), что существенно затрудняет расчет и не всегда приводит к желаемым результатам. Для простых случаев разработаны рекомендации по выбору базисных функций, в сложных случаях проблема назначения этих функций может оказаться не проще решения неходкой задачи.

Идея метода конечных элементов (МКЭ) заключается в том, что вместо поиска единого аналитического представления функции v используют ее кусочно-линейную аппроксимацию, т.е. всю область решения разбивают на подобласти конечных размеров, достаточно малых, чтобы обеспечить требуемую точность линеаризации.

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

. (2)

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

Рис. 1. Одномерные конечные элементы

Узловые значения функций вы можете искать различными способами, исходя из физического смысла задачи. Задача может быть задана в форме дифференциальных уравнений, либо в виде некоторого функционала . В первом случае удобно воспользоваться методом Бубнова – Галеркина, минимизируя ошибку приближенного решения, во втором – каким-нибудь вариационным принципом, обеспечивающим минимизацию функционала. Используя кусочно-линейное представление функции и записывая функционал в виде , вы обеспечите реализацию минимума энергии, приравнивая нулю частные производные:

В результате будет получена система уравнений

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

Унификация вычислительных процедур накладывает некоторые особенности на реализацию МКЭ, для уяснения которых вновь обратимся к решению одномерной задачи (рис. 1).

В общем случае алгоритм МКЭ состоит из четырех этапов:

1. выделение конечных элементов (разбиение заданной области на конечные элементы);

2. определение аппроксимирующей функции для отдельного конечного элемента (определение функции элемента);

3. объединение конечных элементов в ансамбль;

4. определение вектора узловых значений функции.

Первый этап. Первый этап расчета с использованием МКЭ состоит в разбиении области на конечные элементы. Такое разбиение начинают обычно от границы области, стараясь наиболее точно повторить ее конфигурацию, затем производят разбиение внутренних областей. Сначала выделяют достаточно крупные подобласти, которые отличаются по свойствам материала, геометрии, напряженному состоянию и пр. Затем каждую подобласть разбивают на конечные элементы принятой формы, чаще всего треугольные, при этом размеры конечных элементов могут быть приняты различными (рис. 2) в зависимости от требуемой точности описания. Резкого изменения размеров конечных элементов на границах подобластей стараются избегать.

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

Информация о разбиении области на конечные элементы и нумерации узлов является исходной для последующих этапов расчета. Обычно требуется указывать тип конечного элемента, его порядковый номер, номер узла сети, координаты этого узла, значения физических параметров в пределах конечного элемента и др. Такая топологическая информация обычно содержит примерно в шесть раз больше чисел, чем количество узлов сетки разбиения. Для объектов средней сложности объем такой информации измеряется сотнями тысяч. Понятно, что ручной ввод этой информации, помимо чрезвычайной утомительности, может сопровождаться ошибками.

Рис. 2. Пример разбиения подобласти на конечные элементы

Процедуры по вводу информации оптимальной нумерации узлов удается формализовать, что позволило разработать специальные программы, называемые препроцессорами, которые либо существуют в виде автономных программ, либо непосредственно включаются в состав программных комплексов, реализующих МКЭ.

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

Второй этап. Рассмотрим (рис. 1, б) один из выделенных элементов (одномерный симплекс-элемент). Для определения коэффициентов а и b полинома (2) запишем граничные условия:

при

при

Подставляя эти значения в равенство (2), получим систему двух уравнений, решив которую относительно а и b, имеем

.

Подставляя эти коэффициенты в формулу (2), можно записать:

(3)

Члены уравнения (3), заключенные в скобки, называют функциями формы одномерного симплекс-элемента:

. (4)

С учетом обозначений (4) уравнение (3) принимает вид:

, (5)

или в матричной форме

, (6)

где – матрица-строка,

– вектор-столбец.

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

Третий этап – объединение конечных элементов в ансамбль. Основу этого этапа составляет замена номеров узлов i и j в уравнении (5) на номера, присвоенные узлам в процессе разбиения рассматриваемой области. Для примера, приведенного на рис. 1, будем иметь:

Полученная система уравнений представляет собой математическую модель исследуемой задачи. Запишем эту систему в расширенной форме:

(7)

что соответствует векторной записи модели

(7а)

Четвертый этап – определение вектора узловых значений функции. Этот этап не является формальным, и мы не можем ограничиться ссылкой на рис. 5.1 абстрактного примера. Нужно знать конкретное дифференциальное уравнение или конкретный функционал исследуемого явления.

Допустим, что в нашем примере этот функционал имеет простейший вид:

(8)

Разобьем пределы интегрирования в соответствии с делением области на конечные элементы, тогда

Первый член подынтегрального выражения определим, возведя в квадрат равенство (5), второй – путем его дифференцирования с учетом обозначений (4). Выполняя интегрирование и простейшие преобразования, можем записать:

(9)

где

Для того, чтобы определить минимум функционала (9), приравняем нулю частные производные, предварительно развернув выражения сумм для всех шести конечных элементов и заменяя индексы i, j номерами соответствующих узлов:

Таким образом, получена следующая система уравнений:

или

. (10)

Матрицу С называют матрицей жесткости, а вектор – вектором нагрузок. Эти названия пришли из строительной механики. Они отвечают сущности рассматриваемых там задач, поскольку МКЭ первоначально был разработан и в основном сформировался применительно к расчетам строительных конструкций. Затем метод конечных элементов был распространен на другие предметные области, а термины строительной механики сохранились и стали общепринятыми, хотя и отражали иную физическую сущность.

Дальнейший расчет сводится к обращению матрицы С и определению вектора . Зная компоненты вектора и используя систему уравнений (7), можно вычислить значение искомой функции в любой точке. Таким образом, задача полностью решена.

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

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