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

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

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

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

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

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

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

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

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

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

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

то

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

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

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

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

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

${\displaystyle \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 )}.}$

где

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

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

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

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

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

${\displaystyle {\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}}}$

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

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

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

${\displaystyle 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}}}$

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

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

${\displaystyle 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}}}$

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

${\displaystyle ={\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}}}$

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

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

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

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

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

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

Получим

${\displaystyle {\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}}}$

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

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

Граничное условие гласит, что функция в узлах на границе равна нулю. Для узла ${\displaystyle u_{i,j}}$ необходимо вычеркнуть ${\displaystyle i}$-й столбец и ${\displaystyle j}$-ю строку в матрице ${\displaystyle {\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