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

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

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

где f_j(t) – приближенные значения искомой функции в узлах, а j_j(x, y, z) – пробные функции. Пробные функции (называемые также интерполяционнымии или функциями формы), как правило, выбираются из числа кусочно-линейных (или полилинейных).

Рассмотрим сначала одномерный случай. Кусочно-линейные пробные функции при этом будут иметь вид:

Представим эти функции графически, получаем

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

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

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

а для прямоугольной

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

Научившись выписывать пробные функции можно выписать систему уравнений, которую предстоит решать. Как и в классическом методе Галеркина это будут уравнения вида

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

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

Задача: в направлении оси X идет пучок частиц первого вида. Все частицы имеют одинаковую скорость. Эти частицы могут сталкиваться с частицами второго вида, при этом частицы первого вида переходят в некоторое другое состояние или превращаются в другие частицы (например, вследствие ионизации, рекомбинации, возбуждения или чего-нибудь еще). Скорость частиц первого вида считаем неизменной. Вопрос: как меняется вдоль оси X концентрация частиц первого вида в начальном состоянии?

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

где величина < sv> характеризует интенсивность взаимодействия частиц разных видов. В силу неизменности скорости первых частиц получаем простое дифференциальное уравнение относительно концентрации

Такое уравнение легко решается аналитически

а в случае, когда подынтегральное выражение является константой, то еще проще

константа находится из граничных условий.

Другой вариант этой же задачи: какова вероятность того, что частица первого вида пролетит расстояние L без каких-либо изменений? Совершенно аналогичным образом получаем ответ (если подынтегральное выражение – константа)

Теперь возвращаемся к методу конечных элементов. Рассмотрим решение дифференциального уравнения

на отрезке [0, 1] с граничным условием f(0)=1. Пример реального уравнения такого типа мы только что рассмотрели. Рассмотрим случай с тремя узлами в точках 0, 0.5 и 1. Значение в одном из них задано граничным условием. С учетом этого аппроксимацию функции можно записать следующим образом

при этом имеем следующие определения для пробных функций.

Невязку исходного уравнения можно записать в виде

Уравнения, которые предстоит решить, выглядят так

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

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

В последнем случае знак положителен на правом конце и отрицателен на левом. Теперь спокойно получаем систему уравнений

после упрощений окончательно получаем

решая эту систему получаем

Аналитическое решение того же уравнения дает f(0.5)»0.607 и f(1)=0.368.

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

для центральных узлов (не граничных) получаем уравнения вида

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

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

Невязкой такого уравнения будет

Следовательно, задача сводится к решению уравнений вида

А теперь для первого интеграла применим интегрирование по частям, получим

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

Теперь можно вспомнить, что

и уравнения преобразуются к виду

интегралы от пробных функций известны. Для линейной интерполяции с равномерным шагом они выглядят так (для центральных узлов)

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

Мы получили (с точностью до постоянного множителя) разностную схему, аналогичную получаемой из метода конечных разностей.

Спектральный метод (метод Фурье)

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

Вариант метода Галеркина, в котором в качестве пробных используются ортогональные функции, называется спектральным методом. Разложение функции по синусам и косинусам традиционно называется нахождением спектра. Отсюда, видимо, и название метода. Хотя, можно использовать и другие наборы ортогональных функций. В отличие от метода конечных элементов, в спектральном методе пробные функции могут принимать отличные от нуля значения во всей расчетной области.

Итак, ортогональными функциями называются функции, удовлетворяющие условию

На отрезке [-p, p] в качестве таких можно использовать функции:

1, cos x, sin x, cos(2x), sin(2x), …, cos(nx), sin(nx), …

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

С помощью нехитрого преобразования можно построить аналогичный набор функций для произвольного отрезка [a, b]:

1, , , …, , , …

Значения интегралов от квадратов функций изменяются при этом пропорционально длине отрезка. У такого набора функций (синусы и косинусы) есть еще одно преимущество (не всегда наблюдаемое у других наборов ортогональных функций): производные от них (с точностью до постоянных множителей) принадлежат тому же набору функций. Единицу иногда рассматривают как частный случай косинуса с нулевым коэффициентом перед аргументом.

Все остальное осталось от классического метода Галеркина. Рассмотрим, например, применение спектрального метода к уравнению

оно моделирует, например, остывание (прогрев) за счет теплопроводности или «рассасывание» за счет диффузии. Рассмотрим решение этого уравнения на отрезке [-L, L] с начальным условием

и граничными условиями: f(-L, t)=0 и f(L, t)=0. Будем искать решение в виде

продифференцировав данное выражение и подставив в исходное выражение, получаем

после сокращения получаем обыкновенное дифференциальное уравнение для функции c(t)

которое легко решается, в результате получаем решение

Рассмотрим случай более сложных начального и граничных условий. Пусть задано f(x, 0)=j(x), f(-L, t)=b₁ и f(L, t)=b₂. Тогда будем искать функцию в виде

Функцию f₀ можно выбрать так, чтобы она удовлетворяла граничным условиям, например

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

в результате решение исходного уравнения получается в виде

коэффициенты A_ci и A_si должны быть найдены из еще неиспользованного начального условия. Подставляя в начальное условие выражение для искомой функции, после переноса f₀ в левую часть, получаем

Отсюда видно, что искомые коэффициенты являются коэффициентами разложения функции j-f₀ в ряд Фурье. На практике от бесконечных рядов оставляют только некоторое ограниченное количество членов. Таким образом, спектральный метод в ряде случаев позволяет существенно упростить решение задачи по сравнению с классическим методом Галеркина.

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

Естественно, у людей возникало желание совместить точность спектрального метода и простоту сеточных методов. Такие методы были созданы. Один из таких методов называется псевдоспектральным. По своей сути он ближе не к методу Галеркина, а к более простому в реализации методу коллокаций.

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

1) нахождение разложения функции в сумму пробных на текущем слое (нахождение спектра). Эту операцию умеют делать быстро.

2) нахождение пространственных производных функции по разложению.

3) нахождение значений функции на новом слое по явной схеме первого порядка.

Рассмотрим, например, применение псевдоспектрального метода к уже рассматривавшейся задаче. Дано уравнение

Рассмотрим решение этого уравнения на отрезке [-L, L] с начальным условием

и граничными условиями: f(-L, t)=0 и f(L, t)=0. Шаг сетки примем равным L. Таким образом, сетка содержит всего 3 узла в направлении оси X. Из них в двух стоят граничные условия. Будем искать решение в единственном оставшемся узле.

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

В момент времени t=0 имеем следующие значения функции в узлах: f(‑ L, 0)=0, f(0, 0)=1, f(L, 0)=0. Разложение будет выглядеть так

теперь оценим значения производных

Следовательно, в единственном узле, где нам нужно найти производную по времени, имеем

В момент времени t=t получаем

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

метод конечных разностей дает значение, отличающееся от теоретического примерно на четверть

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

На этом с методами взвешенных невязок в рамках этого курса покончено.