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

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

Обычно матрица жёсткости получается разреженной, то есть содержащая большое количество нулей. Для работы с подобным типом матриц созданы специальные библиотеки (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}{2A}}{\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})}{2A}}&S_{02}={\cfrac {(y_{2}-y_{1})(y_{1}-y_{0})+(x_{1}-x_{2})(x_{0}-x_{1})}{2A}}\\S_{10}={\cfrac {(y_{0}-y_{2})(y_{2}-y_{1})+(x_{2}-x_{0})(x_{1}-x_{2})}{2A}}&S_{11}=0&S_{12}={\cfrac {(y_{2}-y_{0})(y_{0}-y_{1})+(x_{0}-x_{2})(x_{1}-x_{0})}{2A}}\\S_{20}={\cfrac {(y_{0}-y_{1})(y_{1}-y_{2})+(x_{1}-x_{0})(x_{2}-x_{1})}{2A}}&S_{21}={\cfrac {(y_{1}-y_{0})(y_{0}-y_{2})+(x_{0}-x_{1})(x_{0}-x_{2})}{2A}}&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_{01}&0&0&S_{02}&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_{11}&0&0&S_{12}&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_{21}&0&0&S_{22}&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}}$

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

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

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

В случае граничных условий первого рода необходимо изменить матрицу ${\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