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

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

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

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

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

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

Для трехмерных областей наиболее употребимы элементы в форме тетраэдра и параллелепипеда, которые также могут иметь прямолинейные и криволинейные границы.

Пример использования МКЭ для расчета одномерного температурного поля в однородном стержне. Пусть имеется стержень длиной L и площадью поперечного сечения S (рис. 1).

Рекомендуемые материалы

Рис. 1. Однородный стержень, находящийся под воздействием теплового потока

Один конец стержня жестко закреплен, и к нему подводится тепловой поток q заданной интенсивности. На свободном конце стержня происходит конвективный теплообмен с внешней средой. Известны коэффициент теплообмена  и температура окружающей среды Т_*. Вдоль боковой поверхности стержень теплоизолирован.

Температурное поле в стержне описывается уравнением теплопроводности, которое в одномерном приближении имеет вид

                                                                                      (1)

Краевые условия определяются уравнениями:

 при х=0;                                                                   (2 а)

 при x=L.                                                      (2 б)

Искомое температурное поле является непрерывной функцией координаты х (рис. 2, а).

Рис. 2. Расчет одномерного температурного поля в однородном стержне методом МКЭ.

В МКЭ стержень разбивается произвольным образом на конечные элементы, которые в данном случае являются отрезками неравной длины. На каждом элементе непрерывная функция Т(х) аппроксимируется некоторой линейной зависимостью, как показано на рис. 2, б (в скобках указаны номера элементов). Аппроксимирующая кусочно-линейная функция определяется через узловые значения T₁-Т₆, которые в общем случае сначала неизвестны и подлежат определению в МКЭ.

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

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

Этап 2. Определение аппроксимирующей функции для каждого элемента (определение функции элемента). На данном этапе значение непрерывной функции  в произвольной точке е-го конечного элемента аппроксимируется полиномом  , где А^(е) - вектор-строка коэффициентов полинома; А_о - свободный член; R=(х, у, z) - вектор координат в рассматриваемой точке.

Задача этапа заключается в определении неизвестного вектора А^(е) и свободного члена А_о. Для этого, используя условие непрерывности функции в узлах, коэффициенты полинома выражают через вектор Ф^(е) узловых значений функции и координаты узлов и, проделав эквивалентные преобразования, получают

                                                        (3)

где N^(е) - матрица-строка, элементы которой называют функциями формы конечного элемента.

Функции формы легко вычисляются в каждой точке конечного элемента через координаты самой точки и координаты узлов элемента.

Этап 3. Объединение конечных элементов в ансамбль. На этом этапе уравнения (4), относящиеся к отдельным элементам, объединяются в ансамбль, т. е. в систему алгебраических уравнений:

                                                               (4)

Система (5) является моделью искомой непрерывной функции.

Этап 4. Определение вектора узловых значений функции. В общем случае вектор Ф в (5) вначале неизвестен. Его определение - наиболее сложная процедура в МКЭ.

Разработано несколько алгоритмов вычисления вектора Ф. Один из алгоритмов основан на минимизации функционала, связанного с физическим смыслом решаемой задачи, он состоит из 4-х этапов.

Найденные значения вектора Ф подставляют в (4), после чего значение функции  легко вычисляется в любой точке заданной области.

Рассмотрим подробно этапы алгоритма МКЭ.

Этап 1: выделение конечных элементов.

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

Разбиение области на элементы обычно начинают от ее границы с целью наиболее точной аппроксимации формы границы, затем производится разбиение внутренних областей. Часто разбиение области на элементы производят в несколько этапов. Сначала область разбивают на достаточно крупные подобласти (подконструкции), границы между которыми проходят там, где изменяются свойства материала, геометрия, приложенная нагрузка и пр. Затем каждая подобласть взбивается на элементы. Резкого изменения размеров конечных элементов на границах подобластей стараются избегать. На рис. 3. приведен пример разбиения двухмерной области произвольной формы на треугольные конечные элементы с криволинейными границами.

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

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

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

Рис. 4. Матрица ленточной структуры.

Целое число L, представляющее собой наибольшую разность между номерами ненулевых элементов в строке, называется шириной полосы. Чем меньше ширина полосы, тем меньший объем ОП требуется для хранения матрицы при реализации МКЭ в САПР и тем меньше затраты машинного времени на решение результирующей системы уравнений. Ширина полосы зависит, в свою очередь, от числа степеней свободы узлов и способа нумерации последних.

Под числом степеней свободы понимают количество неизвестных функций, определяемых в каждом узле. Так, например, для двухмерных задач гидравлики в каждом узле определяют три переменные: давление и составляющие скорости по осям X и Y.

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

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

На рис. 5 приведены два различных способа нумерации узлов произвольной области, разбитой на конечные элементы. При первом способе (рис. 5, а) N=14, при втором (рис. 5, б) N=5. Ширина полосы для представленных способов при одной степени свободы в узле получается равной соответственно 15 и 6; при двух степенях свободы - 30 и 12. Рациональная нумерация в случае рис. 5, б сокращает необходимый объем ОП почти в три раза по сравнению со случаем рис. 5, б.

Рис. 5. Способы нумерации узлов при разбиении двухмерной области на конечные элементы.

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

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

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

Этап 2: определение аппроксимирующей функции элементов.

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

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

Одномерный симплекс-элемент представляет собой отрезок (рис. 6).

Рис. 6. Одномерный симплекс-элемент.

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

                                                                     (5)

Коэффициенты  и  определяются через узловые значения функции  и  в соответствии с условием непрерывности:

                                                               (6)

Подставив (6) в (5), получим систему уравнений:

решая, которую определим  и :

т.е.

Подставив вычисленные значения коэффициентов аппроксимирующего полинома в (5), получим

Проведем эквивалентные преобразования правой части:

(7)

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

;                                         (8)

С учетом обозначений (8) уравнение (7) принимает вид

                                                              (9)

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

                                                                             (10),

где  - матрица-строка;  - вектор-столбец.

Функция формы обладает следующим свойством: функция формы с номером i равна 1 в узле с соответствующим номером и равна 0 во всех других узлах. Не представляет труда убедиться в наличии этого свойства у функций формы (8).

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

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

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

Можно написать следующее соответствие между произвольными номерами i, j, фигурирующими в уравнении (10), и глобальными номерами узлов рассматриваемой дискретной модели для

элемента 1:                                                    (11 а)

элемента 2:                                                   (11 б)

элемента 3:                                                   (11 в)

элемента 4:                                                   (11 г)

элемента 5:                                                   (11 д)

Подставив значения номеров узлов (11) в (9), получим:

                                                       (12)

где верхние индексы в скобках относятся к номеру элемента.

В выражениях для функций формы элемента (9) значения произвольных номеров i, j также следует изменить в соответствии с (11). Тогда значения , , например, определяются по формулам:

Очевидно, что  и  не равны друг другу даже в случае равенства длин элементов L⁽²⁾ и L⁽³⁾. При известных значениях узловых величин Т₁-Т₆ уравнения (12) позволяют определить значение температуры в любой точке стержня.

Расширенная форма описания моделей имеет некоторые преимущества при реализации следующих этапов алгоритма МКЭ.

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

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

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

Рис. 7. Пример минимизации функционала при нахождении распределения температуры в стержне.

При указанном методе минимизируется функционал:

                             (13)

где V - объем тела; S - площадь границы.

В функционал F входят оба граничных условия (2). При минимизации функционала используется множество функций элементов дискретизированной области. Для простоты вычислений будем считать, что стержень разбит всего на два элемента (в практических случаях этого недостаточно) (рис. 12). Тогда

;                                                 (14)

Функционал (13) удобно представить в виде

                            (15)

где S₁ и S₂ - площади сечений стержня, на которых заданы граничные условия (2 а) и (2 б) соответственно.

Для вычисления объемного интеграла в (15) его необходимо разбить на два слагаемых в соответствии с принятой конечной элементной моделью:

            (16)

Производные в (16) вычисляются с учетом (14) и (8)

                                                                              (17)

Подставив (17) в (16) и считая, что  получим

Второе и третье слагаемые в (15) вычисляются просто так как подынтегральным функциям соответствуют узловые значения Т₁ и Т₃:

;

где S₁ и S₂ - площади поверхностей, на которых заданы q и  (для рассматриваемого примера  и ).

Значение функционала F вычисляется простым суммированием последних трех выражений:

           (18)

где  и

Для минимизации функционала F необходимо выполнение условий:

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

                                 (19)

В общем виде (19) можно представить так: , что соответствует (7).

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

Зная характеристики материала, из системы (19) можно определить узловые значения T₁, Т₂, Т₃. Из (18) и (19) нетрудно заметить, что однотипные конечные элементы вносят в эти выражения слагаемые одного вида. Поэтому при реализации МКЭ в САПР вклад элемента определенного типа в матрицу жесткости вычисляется только один раз, а затем используется во всех необходимых случаях.

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

Метод Галеркина основан на минимизации ошибки  приближенного решения u исходного дифференциального уравнения , где L - дифференциальный оператор.

Для минимизации  в заданной области G требуется выполнение равенства  для каждой из функций N_i .

Сочетание метода Галеркина с МКЭ приводит к системе уравнений

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

Применив метод Галеркина к (1), получим

                                                                              (20)

Подставим в (33) формулу дифференцирования произведений:

                 (21)

Интерполяционная функция Т не сохраняет постоянства по длине стержня, поэтому интеграл в (21) можно представить суммой соответствующих интегралов для отдельных элементов. Так, второй интеграл в (21) можно представить в виде

Вычислим в этом уравнении интегралы, относящиеся к отдельным элементам:

                                          (22)

                                   (23)

С учетом (22) и (23)

       (24)

Первый интеграл в (21) на основании теоремы Остроградского-Гаусса преобразуется к виду

                                    (25)

где ; n - внешняя нормаль к рассматриваемой поверхности.

С учетом краевого условия (2 а) в точке для первого элемента интеграл (25) принимает вид

                                   (26)

Подставив значения ,  в (26), получим

                        (27)

С учетом краевого условия (2 б) в точке  второго элемента интеграл (25) запишется так:

                   (28)

Информация в лекции "2 Начало формирования русского централизованного государства" поможет Вам.

Просуммировав выражения вида (24) для первого и второго элементов и выражения (27) и (28) и приравняв сумму нулю, получим систему уравнений

                               (29)

где    

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

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