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

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

Суть

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

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

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

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

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

Определение

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

$A_{ij}=\underset{\mathrm{\Omega }}{\int }\mathrm{\nabla }{\phi }_i\cdot \mathrm{\nabla }{\phi }_{j}dx.$

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

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

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

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

$u\approx u^h=u_1{\phi }_1+\cdots +u_n{\phi }_n.$

$u_i$ — это известные значения функции в узлах, а $\phi$ — некие базисные функции

то

$A_{ij}^{[k]}=\underset{triangle}{\int }\mathrm{\nabla }{\phi }_i\cdot \mathrm{\nabla }{\phi }_{j}dx.$

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

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

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

$\left[\begin{array}{ccc}{S}_{00}& S_{01}& S_{02}\\ S_{10}& S_{11}& S_{12}\\ S_{20}& S_{21}& S_{22}\end{array}\right]$

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

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

${\alpha }_1=\frac{1}{4A}((x_1{y}_2+x_2{y}_1)+(y_1-y_2)x+(x_2-x_1)y).$

где

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

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

$A=\det \left[\begin{array}{ccc}1& x_1& y_1\\ 1& x_2& y_2\\ 1& x_3& y_3\end{array}\right]$

Сами $S_{ij}=\int (\mathrm{\nabla }{\alpha }_i)(\mathrm{\nabla }{\alpha }_j)dSi,j=0,1,2$

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

$\left[\begin{array}{ccc}{S}_{00}=0& S_{01}=\frac{(y_1-y_2)(y_2-y_0)+(x_2-x_1)(x_0-x_2)}{4A}& S_{02}=\frac{(y_2-y_1)(y_1-y_0)+(x_1-x_2)(x_0-x_1)}{4A}\\ S_{10}=\frac{(y_0-y_2)(y_2-y_1)+(x_2-x_0)(x_1-x_2)}{4A}& S_{11}=0& S_{12}=\frac{(y_2-y_0)(y_0-y_1)+(x_0-x_2)(x_1-x_0)}{4A}\\ S_{20}=\frac{(y_0-y_1)(y_1-y_2)+(x_1-x_0)(x_2-x_1)}{4A}& S_{21}=\frac{(y_1-y_0)(y_0-y_2)+(x_0-x_1)(x_0-x_2)}{4A}& S_{22}=0\end{array}\right]$

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

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

Обозначим $u_d^T=\left[\begin{array}{cccccc}{u}_1& u_2& u_3& u_4& u_5& u_6\end{array}\right]$ — вектор-строку значений функции в вершинах двух треугольников (см.рис). Символ $u^T$ обозначает транспонирование матрицы $u.$ То есть, это вектор значений функции в шести узлах треугольников. Очевидно, что при объединении оных получится вектор $u_c$, содержащий только четыре компоненты.

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

$u_d=C{u}_c\iff \left[\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\\ u_5\\ u_6\end{array}\right]=\left[\begin{array}{cccc}1& & & \\ & 1& & \\ & & 1& \\ & 1& & \\ & & & 1\\ 1& & & \end{array}\right]\cdot \left[\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\end{array}\right]$

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

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

$S_d=\left[\begin{array}{cc}{S}^{(1)}& 0\\ 0& S^{(2)}\end{array}\right]\iff \left[\begin{array}{cccccc}{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{array}\right]$

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

$=\left[\begin{array}{cccc}{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{array}\right]$

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

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

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

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

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

$S=\left[\begin{array}{ccc}{S}_{00}& S_{01}& S_{02}\\ S_{10}& S_{11}& S_{12}\\ S_{20}& S_{21}& S_{22}\end{array}\right]$

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

Получим

$𝔖=\left[\begin{array}{cccccccccc}{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{array}\right]$

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

Матрица масс

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

$C_1=\frac{(x_2-x_1)(x_3-x_1)+(y_2-y_1)(y_3-y_1)}{4A}$

$C_2=\frac{(x_3-x_2)(x_1-x_2)+(y_3-y_2)(y_1-y_2)}{4A}$

$C_3=\frac{(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$ записывается в матрицу $𝔖$ любым известным читателю способом. В коде используется метод дозаписи, приведённый выше.

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

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

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

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

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

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

Код

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

Литература

• 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