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

Метод конечных элементов - это численный метод решения дифференциальных уравнений. В данной статье будет показан метод подсчёта и приведён код соответствующей программы на C++ (g++ 5.1.0; linux 4.0.x-ARCH). Для работы также потребуется MATLAB 4.0+ с возможностью вызова из терминала.

Суть

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

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

Матрица жёсткости

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

Обычно матрица жёсткости получается разреженной, то есть содержащая большое количество нулей. Для работы с подобным типом матриц созданы специальные библиотеки (mtl4, SparseLib++, SPARSPAK и другие)

Определение

Элементы матрицы жёсткости $A^{k}$ в общем случае равны

$A_{ij}=\int _{\Omega }\nabla \varphi _{i}\cdot \nabla \varphi _{j}\,dx.$

Например, если дано уравнение Пуассона

$-\nabla ^{2}u=f$

в пространстве $\Omega$ и граничные условния — это $u=0.$

Представим функцию как ряд:

$u\approx u^{h}=u_{1}\varphi _{1}+\cdots +u_{n}\varphi _{n}.$

$u_{i}$ — это известные значения функции в узлах, а $\varphi$ — некие базисные функции

то

$A_{ij}^{[k]}=\int _{triangle}\nabla \varphi _{i}\cdot \nabla \varphi _{j}\,dx.$

Создание матрицы

Для одного треугольника

Пусть дан один конечный элемент, для простоты — треугольный. Матрица жёсткости, по сути, задаёт связи между узлами. Так как у элемента три узла (в локальной нумерации — 0, 1 и 2), то матрица будет иметь вид

${\begin{bmatrix}S_{00}&S_{01}&S_{02}\\S_{10}&S_{11}&S_{12}\\S_{20}&S_{21}&S_{22}\end{bmatrix}}$

В дальнейшем матрицу для одного треугольника будем называть локальной, для всей сетки сразу - глобальной.

В общем случае, элементы $S_{ij}$ определяются через линейные функции

$\alpha _{1}={\cfrac {1}{4A}}{\big (}(x_{1}y_{2}+x_{2}y_{1})\;+\;(y_{1}-y_{2})x\;+\;(x_{2}-x_{1})y{\big )}.$

где

$A$ — площадь треугольного элемента.

$\alpha _{2}$ и $\alpha _{3}$ получаются из $\alpha _{1}$ цикличной перестановкой индексов. Удобно искать $A$ как определитель матрицы

$A=\det {\begin{bmatrix}1&x_{1}&y_{1}\\1&x_{2}&y_{2}\\1&x_{3}&y_{3}\end{bmatrix}}$

Сами $S_{ij}=\int (\nabla \alpha _{i})(\nabla \alpha _{j})dS\;\;\;\;\;i,j=0,1,2$

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

${\begin{bmatrix}S_{00}=0&S_{01}={\cfrac {(y_{1}-y_{2})(y_{2}-y_{0})+(x_{2}-x_{1})(x_{0}-x_{2})}{4A}}&S_{02}={\cfrac {(y_{2}-y_{1})(y_{1}-y_{0})+(x_{1}-x_{2})(x_{0}-x_{1})}{4A}}\\S_{10}={\cfrac {(y_{0}-y_{2})(y_{2}-y_{1})+(x_{2}-x_{0})(x_{1}-x_{2})}{4A}}&S_{11}=0&S_{12}={\cfrac {(y_{2}-y_{0})(y_{0}-y_{1})+(x_{0}-x_{2})(x_{1}-x_{0})}{4A}}\\S_{20}={\cfrac {(y_{0}-y_{1})(y_{1}-y_{2})+(x_{1}-x_{0})(x_{2}-x_{1})}{4A}}&S_{21}={\cfrac {(y_{1}-y_{0})(y_{0}-y_{2})+(x_{0}-x_{1})(x_{0}-x_{2})}{4A}}&S_{22}=0\end{bmatrix}}$

Первый вид обобщения на несколько треугольников - сшивание

Для того, чтобы сделать из многих раздельных матриц, полученных выше, одну большую матрицу, описывающую отношения между узлами всей области расчёта, необходимо произвести процедуру объединения матриц. Пусть символ $d$ обозначает разделённые элементы (рис. а), а символ $c$ — объединённые элементы (рис. б).

Обозначим $u_{d}^{T}={\begin{bmatrix}u_{1}&u_{2}&u_{3}&u_{4}&u_{5}&u_{6}\end{bmatrix}}$ — вектор-строку значений функции в вершинах двух треугольников (см.рис). Символ $u^{T}$ обозначает транспонирование матрицы $u.$ То есть, это вектор значений функции в шести узлах треугольников. Очевидно, что при объединении оных получится вектор $u_{c}$, содержащий только четыре компоненты.

Преобразование происходит по схеме

$u_{d}=Cu_{c}\;\Leftrightarrow \;{\begin{bmatrix}u_{1}\\u_{2}\\u_{3}\\u_{4}\\u_{5}\\u_{6}\end{bmatrix}}\;=\;{\begin{bmatrix}1&&&\\&1&&\\&&1&\\&1&&\\&&&1\\1&&&\end{bmatrix}}\cdot {\begin{bmatrix}u_{1}\\u_{2}\\u_{3}\\u_{4}\\\end{bmatrix}}$

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

Запишем теперь матрицу жёсткости для двух треугольников:

$S_{d}={\begin{bmatrix}S^{(1)}&0\\0&S^{(2)}\end{bmatrix}}\;\Leftrightarrow \;{\begin{bmatrix}S_{00}&S_{01}&S_{02}&&&\\S_{10}&S_{11}&S_{12}&&&\\S_{20}&S_{21}&S_{22}&&&\\&&&S_{33}&S_{34}&S_{35}\\&&&S_{43}&S_{44}&S_{45}\\&&&S_{53}&S_{54}&S_{55}\\\end{bmatrix}}$

Результирующая матрица $S_{global}=C^{T}S_{d}C=$

$={\begin{bmatrix}S_{00}^{(1)}+S_{55}^{(2)}&S_{01}^{(1)}+S_{53}^{(2)}&S_{02}^{(1)}&S_{54}^{(2)}\\S_{10}^{(1)}+S_{35}^{(2)}&S_{11}^{(1)}+S_{33}^{(2)}&S_{12}^{(1)}&S_{34}^{(2)}\\S_{20}^{(1)}&S_{20}^{(1)}&S_{22}^{(1)}&0\\S_{45}^{(2)}&S_{43}^{(2)}&0&S_{44}^{(2)}\\\end{bmatrix}}$

То есть, на каждом следующем шаге необходимо добавлять новые элементы к уже существующим.

Второй вид обобщения на несколько треугольников - дозаписывание

Пусть есть область, представленная и разбитая на треугольники так, как преставлено на рисунке. Пусть данная сетка содержит $N$ узлов. Создадим глобальную матрицу ${\mathfrak {S}}$ (размера, очевидно, $N\times N$) и заполним её нулями. Начнём строить локальные матрицы $S$ для треугольников, например, для $\Delta 036.$

Введём локальную нумерацию для данного треугольника: пусть его верхняя вершина имеет локальный номер $0$, далее по часовой стрелке $1$ и $2$. Иначе говоря, пусть глобальным номерам $0,3,6$ соответствуют локальные номера $0,1,2$ соответственно.

Составим матрицу для этого треугольника так, как описано выше, получив что-то типа

$S={\begin{bmatrix}S_{00}&S_{01}&S_{02}\\S_{10}&S_{11}&S_{12}\\S_{20}&S_{21}&S_{22}\end{bmatrix}}$

Теперь заменим локальную нумерацию на глобальную. То есть запишем локальное число $S_{00}$ как глобальное число ${\mathfrak {S_{00}}}$, $S_{01}$ - как ${\mathfrak {S_{03}}}$, $S_{02}$ - как ${\mathfrak {S_{06}}}$ и так далее.

Получим

${\mathfrak {S}}={\begin{bmatrix}S_{00}&0&0&S_{03}&0&0&S_{06}&0&0&0\\0&0&0&0&0&0&0&0&0&0\\0&0&0&0&0&0&0&0&0&0\\S_{10}&0&0&S_{13}&0&0&S_{16}&0&0&0\\0&0&0&0&0&0&0&0&0&0\\0&0&0&0&0&0&0&0&0&0\\S_{20}&0&0&S_{23}&0&0&S_{26}&0&0&0\\0&0&0&0&0&0&0&0&0&0\\0&0&0&0&0&0&0&0&0&0\\0&0&0&0&0&0&0&0&0&0\end{bmatrix}}$

С остальными треугольниками поступаем аналогично. Необходимо помнить, что надо "дописать" число в глобальную ячейку, то есть прибавить к уже существующему.

Матрица масс

Матрица масс собирается по тем же правилам, но чуть по другим формулам. Создаётся матрица $S$ размеров три на три, затем говорится, что

$C_{1}={\cfrac {(x_{2}-x_{1})(x_{3}-x_{1})+(y_{2}-y_{1})(y_{3}-y_{1})}{4A}}$

$C_{2}={\cfrac {(x_{3}-x_{2})(x_{1}-x_{2})+(y_{3}-y_{2})(y_{1}-y_{2})}{4A}}$

$C_{3}={\cfrac {(x_{1}-x_{3})(x_{2}-x_{3})+(y_{1}-y_{3})(y_{2}-y_{3})}{4A}}$

• $S_{22}=S_{22}+C_{1}$

$S_{23}=S_{23}-C_{1}$

$S_{32}=S_{32}-C_{1}$

$S_{33}=S_{33}+C_{1}$

• $S_{33}=S_{33}+C_{2}$

$S_{31}=S_{31}-C_{2}$

$S_{13}=S_{13}-C_{2}$

$S_{11}=S_{11}+C_{2}$

• $S_{11}=S_{11}+C_{3}$

$S_{12}=S_{12}-C_{3}$

$S_{21}=S_{21}-C_{3}$

$S_{22}=S_{22}+C_{3}$

где

• $A$ - площадь данного треугольника, которая считается, как в предыдущей главе.
• $C_{2}$ и $C_{3}$ получаются из $C_{1}$ циклической перестановкой, равно как второй и третий блок элементов матрицы из первого.

После чего полученная матрица $S$ записывается в матрицу ${\mathfrak {S}}$ любым известным читателю способом. В коде используется метод дозаписи, приведённый выше.

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

Условия Дирихле

В случае граничных условий первого рода необходимо изменить матрицу ${\mathfrak {S}}.$

Граничное условие гласит, что функция в узлах на границе равна нулю. Для узла $u_{i,j}$ необходимо вычеркнуть $i$-тый столбец и $j$-ую строку в матрице ${\mathfrak {S}},$ а также вычеркнуть сам узел из массива узлов решётки.

Условия Неймана

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

Код

Из-за большого объёма кода он был помещён на подстраницу.

Литература

• Zienkiewicz, K. Morgan — Finite elements and approximation, 1983.
• P.P. Silvester, R.L. Ferrari — Finite elements for electrical engineers, 1986.
• E. Suli — Finite Element Methods for Partial Differential Equations, 2011
• K.J. Bathe, E.L. Wilson — Numerical methods in finite element analysis, 1982