Теория метода конечных элементов (FEM)

Мотивация

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

Нужен более универсальный метод - таким и является Finite Element Method.

Общая формулировка метода конечных элементов (FEM)

Постановка математической задачи. Слабая форма

Дифференциальная задача на неизвестную функцию \(\mathbf{u}\) в самом общем случае может быть записана как

\[\begin{split} \mathbf{A}(\mathbf{u})=\left\{\begin{array}{c} A_1(\mathbf{u}) \\ A_2(\mathbf{u}) \\ \vdots \end{array}\right\}=\mathbf{0}, \quad x \in \Omega \end{split}\]

с учётом граничных условий:

\[\begin{split} \mathbf{B}(\mathbf{u})=\left\{\begin{array}{c} B_1(\mathbf{u}) \\ B_2(\mathbf{u}) \\ \vdots \end{array}\right\}=\mathbf{0}, \quad x \in \Gamma = \partial \Omega \end{split}\]

Искомая функция \(\mathbf{u}\) может быть как скаляром, так и вектором. Системы выше могут состоять как из одного дифференциального уравнения, так и из нескольких (система УРЧП), в том числе быть нелинейными по \(\mathbf{u}\).

👉Примеры

Пример 1. Уравнения установившейся теплопроводности в двумерной области

\[\begin{split} \begin{array}{rlr} A(\phi) & =\frac{\partial}{\partial x}\left(k \frac{\partial \phi}{\partial x}\right)+\frac{\partial}{\partial y}\left(k \frac{\partial \phi}{\partial y}\right)+Q=0 \\ B(\phi) & =\phi-\bar{\phi}=0 & \\ \text{или} \quad B(\phi) & =k \frac{\partial \phi}{\partial n}+\bar{q}=0 & \text{на } \Gamma_\phi \end{array} \end{split}\]

где \(\mathbf{u} \equiv \phi\) обозначает температуру, \(k\) — теплопроводность, \(Q\) — источник тепла, \(\bar{\phi}\) и \(\bar{q}\) — заданные значения температуры и теплового потока на границах, а \(n\) — направление, нормальное к границе \(\Gamma\).

В приведённой задаче \(k\) и \(Q\) могут зависеть от координаты, а при нелинейной постановке задачи — также от \(\phi\) или её производных.

---

Пример 2. Уравнение установившейся теплопроводности-конвекции в двумерной области

\[ A(\phi)=\frac{\partial}{\partial x}\left(k \frac{\partial \phi}{\partial x}\right)+\frac{\partial}{\partial y}\left(k \frac{\partial \phi}{\partial y}\right)+u_x \frac{\partial \phi}{\partial x}+u_y \frac{\partial \phi}{\partial y}+Q=0 \]

с граничными условиями, как в первом примере. Здесь \(u_x\) и \(u_y\) — известные функции координат, представляющие собой компоненты скорости несжимаемой жидкости, в которой происходит теплоперенос.

---

Пример 3. Система трёх уравнений первого порядка, эквивалентная Примеру 1

\[\begin{split} \begin{aligned} \mathbf{A}(\mathbf{u}) = \left\{ \begin{array}{c} \frac{\partial q_x}{\partial x} + \frac{\partial q_y}{\partial y} + Q \\ q_x + k \frac{\partial \phi}{\partial x} \\ q_y + k \frac{\partial \phi}{\partial y} \end{array} \right\} = 0 \end{aligned} \end{split}\]

в области \(\Omega\), и

\[\begin{split} \begin{aligned} \mathbf{B}(\mathbf{u}) & = \phi - \bar{\phi} = 0 & & \text{на } \Gamma_\phi \\ & = q_n - \bar{q} = 0 & & \text{на } \Gamma_q \end{aligned} \end{split}\]

где \(q_n\) — это тепловой поток, нормальный к границе.

Вектор неизвестной функции \(\mathbf{u}\) имеет вид:

\[\begin{split} \mathbf{u} = \left\{ \begin{array}{c} \phi \\ q_x \\ q_y \end{array} \right\} \end{split}\]

Этот пример типичен для так называемой смешанной формулировки (mixed formulation). В таких задачах число зависимых переменных можно уменьшить с помощью подходящих алгебраических преобразований, при этом задача остаётся решаемой (например, получение уравнения (пример 1) из (пример 3) путём исключения \(q_x\) и \(q_y\)).

Если такое исключение невозможно (см. уравнение (пример 1)), то формулировка называется неприводимой (irreducible formulation).

Если \(\mathbf{u}\) - решение дифференциальной задачи, заметим, что оно удовлетворяет следующему выражению для любых \(\mathbf{v}\):

\[ \int_{\Omega} \mathbf{v}^{\mathrm{T}} \mathbf{A}(\mathbf{u}) \mathrm{d} \Omega \equiv \int_{\Omega}\left[v_1 A_1(\mathbf{u})+v_2 A_2(\mathbf{u})+\cdots\right] \mathrm{d} \Omega \equiv 0 \]

где

\[\begin{split} \mathbf{v}=\left\{\begin{array}{c} v_1 \\ v_2 \\ \vdots \end{array}\right\} \end{split}\]

вектор из произвольных функций.

С учётом граничных условий, имеем

\[ \int_{\Omega} \mathbf{v}^{\mathrm{T}} \mathbf{A}(\mathbf{u}) \mathrm{d} \Omega+\int_{\Gamma} \overline{\mathbf{v}}^{\mathrm{T}} \mathbf{B}(\mathbf{u}) \mathrm{d} \Gamma=0 \quad \forall \, \mathbf{v} \]

Таким образом возникает идея решать именно эту задачу (вариационную задачу) вместо решения дифференциальной. Т.е. мы ищем \(\mathbf{u}\) исходя из выполнения уравнения выше для любых функций \(\mathbf{v}\) (т.н. тестовая функция)).

Примечание (слабая форма)

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

\[ \int_{\Omega} \mathbf{C}(\mathbf{v})^{\mathrm{T}} \mathbf{D}(\mathbf{u}) \, \mathrm{d} \Omega + \int_{\Gamma} \mathbf{E}(\overline{\mathbf{v}})^{\mathrm{T}} \mathbf{F}(\mathbf{u}) \, \mathrm{d} \Gamma = 0 \]

В этом выражении операторы \(\mathbf{C}\)–\(\mathbf{F}\) обычно содержат производные более низкого порядка, чем операторы \(\mathbf{A}\) или \(\mathbf{B}\). Теперь требуется меньший порядок непрерывности при выборе функции \(\mathbf{u}\), но это достигается за счёт необходимости более высокой непрерывности для тестовых функций \(\mathbf{v}\) и \(\overline{\mathbf{v}}\).

Такая форма уравнения называется слабой формой.

👉Пример

Пример. Слабая форма уравнения теплопроводности из примера 1 — принудительные и естественные граничные условия

Рассмотрим теперь интегральную форму дифференциального уравнения теплопроводности. Мы можем записать его в следующем виде:

\[ \int_{\Omega} v\left[\frac{\partial}{\partial x}\left(k \frac{\partial \phi}{\partial x}\right)+\frac{\partial}{\partial y}\left(k \frac{\partial \phi}{\partial y}\right)+Q\right] \, \mathrm{d}x \, \mathrm{d}y + \int_{\Gamma_q} \bar{v}\left[k \frac{\partial \phi}{\partial n} + \bar{q}\right] \, \mathrm{d} \Gamma = 0 \]

Здесь \(v\) и \(\bar{v}\) — тестовые скалярные функции, причём предполагается, что одно из граничных условий, а именно:

\[ \phi - \bar{\phi} = 0 \]

автоматически удовлетворяется благодаря соответствующему выбору функции \(\phi\) на границе \(\Gamma_\phi\).

Это уравнение можно проинтегрировать по частям, чтобы получить его слабую форму. Для этого используем обобщённые формулы Грина:

\[\begin{split} \begin{aligned} \int_{\Omega} v \frac{\partial}{\partial x}\left(k \frac{\partial \phi}{\partial x}\right) \, \mathrm{d}x \, \mathrm{d}y & = -\int_{\Omega} \frac{\partial v}{\partial x} \left(k \frac{\partial \phi}{\partial x}\right) \, \mathrm{d}x \, \mathrm{d}y + \oint_{\Gamma} v \left(k \frac{\partial \phi}{\partial x}\right) n_x \, \mathrm{d}\Gamma \\ \int_{\Omega} v \frac{\partial}{\partial y}\left(k \frac{\partial \phi}{\partial y}\right) \, \mathrm{d}x \, \mathrm{d}y & = -\int_{\Omega} \frac{\partial v}{\partial y} \left(k \frac{\partial \phi}{\partial y}\right) \, \mathrm{d}x \, \mathrm{d}y + \oint_{\Gamma} v \left(k \frac{\partial \phi}{\partial y}\right) n_y \, \mathrm{d}\Gamma \end{aligned} \end{split}\]

Подставив это в исходное выражение, получим следующую слабую форму:

\[\begin{split} \begin{aligned} & -\int_{\Omega} \left( \frac{\partial v}{\partial x} \cdot k \frac{\partial \phi}{\partial x} + \frac{\partial v}{\partial y} \cdot k \frac{\partial \phi}{\partial y} - v Q \right) \, \mathrm{d}x \, \mathrm{d}y \\ & \quad + \oint_{\Gamma} v k\left(\frac{\partial \phi}{\partial x} n_x + \frac{\partial \phi}{\partial y} n_y\right) \, \mathrm{d} \Gamma + \int_{\Gamma_q} \bar{v} \left[k \frac{\partial \phi}{\partial n} + \bar{q} \right] \, \mathrm{d} \Gamma = 0 \end{aligned} \end{split}\]

Заметим, что производная по нормали выражается как:

\[ \frac{\partial \phi}{\partial n} = \frac{\partial \phi}{\partial x} n_x + \frac{\partial \phi}{\partial y} n_y \]

Рассмотрим два важных замечания:

(a) Переменная \(\phi\) исчезает из интегралов по границе \(\Gamma_q\), и соответствующее граничное условие

\[ k \frac{\partial \phi}{\partial n} + \bar{q} = 0 \]

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

(b) Если выбор функции \(\phi\) ограничен так, чтобы удовлетворять принудительным граничным условиям \(\phi = \bar{\phi}\), то мы можем исключить последний интеграл, выбрав такие функции \(v\), что \(v = 0\) на \(\Gamma_\phi\).

Таким образом, полученное выражение представляет собой слабую форму уравнения теплопроводности, эквивалентную исходному дифференциальному уравнению. Эта форма допускает разрывные коэффициенты теплопроводности \(k\) и температурные поля \(\phi\), имеющие разрывные первые производные — что невозможно в дифференциальной форме.

Если без потери общности положить

\[ \bar{v} = -v \quad \text{на } \Gamma, \]

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

\[ \int_{\Omega} \nabla^{\mathrm{T}} v \, k \nabla \phi \, \mathrm{d} \Omega - \int_{\Omega} v Q \, \mathrm{d} \Omega - \int_{\Gamma_q} v \bar{q} \, \mathrm{d} \Gamma - \int_{\Gamma_\phi} v k \frac{\partial \phi}{\partial n} \, \mathrm{d} \Gamma = 0 \]

где оператор градиента записывается как:

\[\begin{split} \nabla = \left\{ \begin{array}{c} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \end{array} \right\} \end{split}\]

Хотим опять свести эту задачу к решению СЛАУ. Для этого разложим нашу функцию по базису \(\phi_i(\vec x)\) (конкретный вид базиса выберем позже):

\[ u(\vec x) \approx u^N(\vec x)=\sum_{i=1}^N u_i \cdot \phi_i(\vec x), \quad u_i \in \mathbb R \]

При подстановке в вариационную задачу у нас возникает линейное уравнение на \(u_i\), коэффициенты в котором выражаются через интегралы от \(\phi_i(\vec x)\) и тестовой функции \(v(\vec{x})\). Если они возьмутся аналитически, то будет просто линейное уравнение.

Т.к. вариационная задача верна для любых тестовых функций \(\mathbf{v}\), то она верна любого базисного набора функций (строго говоря, работаем в гильбертовом пространстве). А значит, вариационная задача верна на любом подмножестве \(\Omega_i\).

Определение. Сетка или mesh

Сетка или mesh - это представление области \(\Omega\) в виде конечного набора однотипных простых множеств. Ранее мы рассматривали только прямоугольные сетки одинакового размера, но в FEM удобнее использовать неравномерные треугольные/прямоугольные сетки для более точного описания геометрии \(\Omega\).

Задавая сетку, мы переходим от одного уравнения ко многим на каждом «кусочке» (элементе). Осталось задать удобный базис \(\phi_i(\vec x)\), чтобы можно было аналитически взять упомянутые ранее интегралы.

Базисные функции \(\phi_i(\vec x)\) не обязаны быть ортогональными, главное, чтобы коэффициенты разложения определялись единственным образом. Пример базиса в 1D - линейные функции на каждом отрезке:

Basis	Image
$\(\phi_j(x)= \begin{cases}0 & \text { if } x \notin\left[x_{j-1}, x_{j+1}\right] \\ \frac{x-x_{j-1}}{h} & \text { if } x \in\left[x_{j-1}, x_j\right] \\ \frac{x_{j+1}-x}{h} & \text { if } x \in\left[x_j, x_{j+1}\right]\end{cases}\)$

Для любой сеточно-заданной функции \(u^h\) такой набор функций будет базисом на узлах элементов. В общем случае (на произвольном mesh), базисные функции будут похожи на одномерный случай - главное выполнение следующего свойства на узлах элементов:

\[\begin{split} \phi_j\left(x_i\right)= \begin{cases}1 & i=j \\ 0 & i \neq j\end{cases} \end{split}\]

Таким образом, на каждом элементе количество базисных функций \(\phi_i(\vec x)\) должно равняться числу узлов, на которых ищем сеточную функцию \(u^h\).

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

Определение. Конечный элемент или finite element

Конечный элемент или finite element - это тройка \(\left(K, \, P_K, \, \Sigma\right)\) где

• \(K\) - элемент сетки (ex. треугольничек или интервал);

• \(P_K\) - конечномерное линейное пространство полиномов, определенное на \(K\) (две прямые в примере выше);

• \(\Sigma\) - набор степеней свободы (DoF), т.е. значения полинома в вершинах \(K\).

Эти три составляющие определяют единственный полином на элементе \(K\).

Алгоритм решения задачи методом конечных элементов (FEM)

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

Алгоритм FEM

1. Постановка дифференциальной задачи (сильная форма)

   • Сформулировать уравнение или систему уравнений в области \(\Omega\), например:

     \[ \mathbf{A}(\mathbf{u}) = \mathbf{0} \quad \text{в } \Omega. \]

   • Задать граничные условия на границе \(\Gamma = \partial \Omega\), например:

     \[ \mathbf{B}(\mathbf{u}) = \mathbf{0} \quad \text{на } \Gamma. \]

   • Отметить, какие условия относятся к типу Дирихле (принудительные) и какие — к типу Неймана (естественные).

2. Переход к слабой (вариационной) постановке

   • Умножить исходное уравнение на тестовые функции \(\mathbf{v}\) и проинтегровать по области (и, при необходимости, по границе):

     \[ \int_{\Omega} \mathbf{v}^{\mathrm{T}} \,\mathbf{A}(\mathbf{u}) \, d\Omega \;+\; \int_{\Gamma} \overline{\mathbf{v}}^{\mathrm{T}} \,\mathbf{B}(\mathbf{u}) \, d\Gamma \;=\; 0. \]

   • По необходимости выполнить интегрирование по частям, чтобы снизить порядок производных у \(\mathbf{u}\).

   • Учесть граничные условия:

     • Принудительные (Dirichlet) обычно «вшиваются» в выбор функции \(\mathbf{u}\),

     • Естественные (Neumann и др.) автоматически войдут в вариационную формулировку через граничный член.

3. Построение сетки (разбиение области на конечные элементы)

   • Разбить область \(\Omega\) на конечное число «элементов» (треугольников, четырёхугольников, …).

   • Определить «тип элемента»: выбор пространства полиномов \(P_K\), степеней свободы и т.д.

4. Выбор базиса и дискретизация решения

   • На каждом элементе задать базисные функции \(\phi_i(\mathbf{x})\), обладающие нужными свойствами (например, принимают значение \(1\) в своём узле и \(0\) в соседних).

   • Приближённо представить искомое решение в виде:

     \[ \mathbf{u}^N(\mathbf{x}) \;=\; \sum_{i=1}^{N} u_i \,\phi_i(\mathbf{x}), \]

     где \(u_i\) — искомые числовые коэффициенты.

5. Составление системы уравнений (СЛАУ)

   • Подставить \(\mathbf{u}^N(\mathbf{x})\) в слабую форму.

   • Выбрать тестовые функции (обычно те же \(\phi_j\), т.н. метод Галёркина).

   • Для каждого базисного \(\phi_j\) получится уравнение:

     \[ \sum_{i=1}^N u_i \int_{\Omega} \bigl[\dots\bigr] \, d\Omega \;+\; \sum_{i=1}^N u_i \int_{\Gamma} \bigl[\dots\bigr] \, d\Gamma =\text{(члены от источников и граничных значений)}. \]

   • Собрать все уравнения в глобальную матрицу (матрица жёсткости и т.д.) и в глобальный вектор (правая часть).

6. Решение системы уравнений

   • Получить СЛАУ вида:

     \[ \mathbf{A}\,\mathbf{u} = \mathbf{f}. \]

   • Решить эту систему (линейными или нелинейными методами, если задача нелинейная).

   • Найти коэффициенты \(u_i\).

7. Постобработка и оценка результата

   • Восстановить приближённое решение \(\mathbf{u}^N(\mathbf{x})\) из найденных \(u_i\).

   • Построить графики/контурные карты/векторные поля (в зависимости от задачи).

   • При необходимости выполнить оценку погрешности и адаптивную перестройку сетки, повторив процесс.

Такова общая схема решения (FEM) для стационарных задач в слабой форме.

Пример. Как-то так выглядит решение одномерного уравнения Пуассона с \(f=-1\) при разных сетках \(h\). Ясно видим отличительную особенность метода FEM - вне сетки решение получается нефизичным, но в узлах сетки сходимость крайне быстрая.

Пример. То же самое уравнение, но с другим базисом. На этот раз используются параболы. За деталями см. дополнительные материалы.

Аппроксимация и устойчивость

Аппроксимация

Если базисные функции - полиномы р степени, а решение нашей задачи \(\mathrm{p}+1\) раз дифференциируемо, то

\[ \left\|u-u_h\right\|_{L^2(\Omega)} \leqslant c \cdot h^{p+1}\cdot \|u\|_{H^{p+1}(\Omega)}, \]

где

\[ \|f\|_{H^p}=\left(\int_{\Omega}|f|^2+|\nabla f|^2+\ldots+\left|\nabla^p f\right|^2\right)^{\frac{1}{p}} \]

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

В случае с уравнением Пуассона получаем, что \(\mathrm{p}=1\) и

\[ \left\|u-u_h\right\|_{L^2(\Omega)} \leqslant c \cdot h^2 \cdot \|u\|_{H^2(\Omega)}, \]

Устойчивость

Устойчивость решения задачи определяется устойчивостью решения системы линейных уравнений. При плохо выбранном базисе мы получим плохо обусловленную матрицу.