Русская Википедия:Матрица жёсткости

Материал из Онлайн справочника
Перейти к навигацииПерейти к поиску

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

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

Определение

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

<math> A_{ij} = \int_\Omega\nabla\varphi_i\cdot\nabla\varphi_j\, dx.</math>

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

<math> -\nabla^2 u = f</math>

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

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

<math> u \approx u^h = u_1\varphi_1+\cdots+u_n\varphi_n.</math>
<math>u_i</math> — это известные значения функции в узлах, а <math>\varphi</math> — некие базисные функции

то

<math> A^{[k]}_{ij} = \int_{triangle}\nabla\varphi_i\cdot\nabla\varphi_j\, dx.</math>

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

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

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

<math>

\begin{bmatrix} S_{00} & S_{01} & S_{02}\\ S_{10} & S_{11} & S_{12}\\ S_{20} & S_{21} & S_{22} \end{bmatrix} </math>

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

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

<math>\alpha_1 = \cfrac {1} {2A} \big( (x_1y_2 + x_2y_1) \; + \; (y_1 - y_2)x \; + \; (x_2 - x_1)y \big) . </math>
где
<math>A</math> — площадь треугольного элемента.

<math>\alpha_2</math> и <math>\alpha_3</math> получаются из <math>\alpha_1</math> цикличной перестановкой индексов. Удобно искать <math>A</math> как определитель матрицы

<math>

A = \det \begin{bmatrix} 1 & x_1 & y_1 \\ 1 & x_2 & y_2 \\ 1 & x_3 & y_3 \end{bmatrix} </math>

Сами <math>S_{ij} = \int (\nabla \alpha_i) (\nabla \alpha_j) dS \;\;\;\;\; i, j = 0, 1, 2</math>

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

<math>

\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} </math>

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

Файл:Stiffness matrix - adding triangles.png

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

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

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

<math>

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}

</math>

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

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

<math> 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}

</math>

Результирующая матрица <math>S_{global} = C^T S_d C = </math>

<math> =
   \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}

</math>

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

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

Файл:Mesh FEM.png

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

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

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

<math>

S =

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

   \end{bmatrix}

</math>

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

Получим

<math>

\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}

</math>

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

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

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

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

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

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

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

Литература

  • 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