Матрица жёсткости
Матрица жёсткости (матрица Дирихле) — матрица особого вида, использующаяся в методе конечных элементов для решения дифференциальных уравнений в частных производных. Она применяется при решениях задач электродинамики и механики.
Обычно матрица жёсткости получается разреженной, то есть содержащая большое количество нулей. Для работы с подобным типом матриц созданы специальные библиотеки (mtl4, SparseLib++, SPARSPAK и другие)
Определение
Элементы матрицы жёсткости A k {displaystyle A^{k}} в общем случае равны
A i j = ∫ Ω ∇ φ i ⋅ ∇ φ j d x . {displaystyle A_{ij}=int _{Omega } abla varphi _{i}cdot abla varphi _{j},dx.}Например, если дано уравнение Пуассона
− ∇ 2 u = f {displaystyle - abla ^{2}u=f}в пространстве Ω {displaystyle Omega } и граничные условния — это u = 0. {displaystyle u=0.}
Представим функцию как ряд:
u ≈ u h = u 1 φ 1 + ⋯ + u n φ n . {displaystyle uapprox u^{h}=u_{1}varphi _{1}+cdots +u_{n}varphi _{n}.} u i {displaystyle u_{i}} — это известные значения функции в узлах, а φ {displaystyle varphi } — некие базисные функциито
A i j [ k ] = ∫ t r i a n g l e ∇ φ i ⋅ ∇ φ j d x . {displaystyle A_{ij}^{[k]}=int _{triangle} abla varphi _{i}cdot abla varphi _{j},dx.}Создание матрицы
Для одного треугольника
Пусть дан один конечный элемент, для простоты — треугольный. Матрица жёсткости, по сути, задаёт связи между узлами. Так как у элемента три узла (в локальной нумерации — 0, 1 и 2), то матрица будет иметь вид
[ S 00 S 01 S 02 S 10 S 11 S 12 S 20 S 21 S 22 ] {displaystyle {egin{bmatrix}S_{00}&S_{01}&S_{02}S_{10}&S_{11}&S_{12}S_{20}&S_{21}&S_{22}end{bmatrix}}}В дальнейшем матрицу для одного треугольника будем называть локальной, для всей сетки сразу - глобальной.
В общем случае, элементы S i j {displaystyle S_{ij}} определяются через линейные функции
α 1 = 1 4 A ( ( x 1 y 2 + x 2 y 1 ) + ( y 1 − y 2 ) x + ( x 2 − x 1 ) y ) . {displaystyle alpha _{1}={cfrac {1}{4A}}{ig (}(x_{1}y_{2}+x_{2}y_{1});+;(y_{1}-y_{2})x;+;(x_{2}-x_{1})y{ig )}.} где A {displaystyle A} — площадь треугольного элемента.α 2 {displaystyle alpha _{2}} и α 3 {displaystyle alpha _{3}} получаются из α 1 {displaystyle alpha _{1}} цикличной перестановкой индексов. Удобно искать A {displaystyle A} как определитель матрицы
A = det [ 1 x 1 y 1 1 x 2 y 2 1 x 3 y 3 ] {displaystyle A=det {egin{bmatrix}1&x_{1}&y_{1}1&x_{2}&y_{2}1&x_{3}&y_{3}end{bmatrix}}}Сами S i j = ∫ ( ∇ α i ) ( ∇ α j ) d S i , j = 0 , 1 , 2 {displaystyle S_{ij}=int ( abla alpha _{i})( abla alpha _{j})dS;;;;;i,j=0,1,2}
В описываемом случае для каждого треугольника составляется такая матрица:
[ S 00 = 0 S 01 = ( y 1 − y 2 ) ( y 2 − y 0 ) + ( x 2 − x 1 ) ( x 0 − x 2 ) 4 A S 02 = ( y 2 − y 1 ) ( y 1 − y 0 ) + ( x 1 − x 2 ) ( x 0 − x 1 ) 4 A S 10 = ( y 0 − y 2 ) ( y 2 − y 1 ) + ( x 2 − x 0 ) ( x 1 − x 2 ) 4 A S 11 = 0 S 12 = ( y 2 − y 0 ) ( y 0 − y 1 ) + ( x 0 − x 2 ) ( x 1 − x 0 ) 4 A S 20 = ( y 0 − y 1 ) ( y 1 − y 2 ) + ( x 1 − x 0 ) ( x 2 − x 1 ) 4 A S 21 = ( y 1 − y 0 ) ( y 0 − y 2 ) + ( x 0 − x 1 ) ( x 0 − x 2 ) 4 A S 22 = 0 ] {displaystyle {egin{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}=0end{bmatrix}}}Первый вид обобщения на несколько треугольников
Для того, чтобы сделать из многих раздельных матриц, полученных выше, одну большую матрицу, описывающую отношения между узлами всей области расчёта, необходимо произвести процедуру объединения матриц. Пусть символ d {displaystyle d} обозначает разделённые элементы (рис. а), а символ c {displaystyle c} — объединённые элементы (рис. б).
Обозначим u d T = [ u 1 u 2 u 3 u 4 u 5 u 6 ] {displaystyle u_{d}^{T}={egin{bmatrix}u_{1}&u_{2}&u_{3}&u_{4}&u_{5}&u_{6}end{bmatrix}}} — вектор-строку значений функции в вершинах двух треугольников (см.рис). Символ u T {displaystyle u^{T}} обозначает транспонирование матрицы u . {displaystyle u.} То есть, это вектор значений функции в шести узлах треугольников. Очевидно, что при объединении оных получится вектор u c {displaystyle u_{c}} , содержащий только четыре компоненты.
Преобразование происходит по схеме
u d = C u c ⇔ [ u 1 u 2 u 3 u 4 u 5 u 6 ] = [ 1 1 1 1 1 1 ] ⋅ [ u 1 u 2 u 3 u 4 ] {displaystyle u_{d}=Cu_{c};Leftrightarrow ;{egin{bmatrix}u_{1}u_{2}u_{3}u_{4}u_{5}u_{6}end{bmatrix}};=;{egin{bmatrix}1&&&&1&&&&1&&1&&&&&11&&&end{bmatrix}}cdot {egin{bmatrix}u_{1}u_{2}u_{3}u_{4}end{bmatrix}}}Нумерация, конечно же, произвольная: необходимо равенство функции в соответствующих вершинах. Матрицу C {displaystyle C} называют матрицей преобразования, а само уравнение называют связанной системой.
Запишем теперь матрицу жёсткости для двух треугольников:
S d = [ S ( 1 ) 0 0 S ( 2 ) ] ⇔ [ 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 ] {displaystyle S_{d}={egin{bmatrix}S^{(1)}&0 &S^{(2)}end{bmatrix}};Leftrightarrow ;{egin{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 g l o b a l = C T S d C = {displaystyle S_{global}=C^{T}S_{d}C=}
= [ 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 ) ] {displaystyle ={egin{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)}&0S_{45}^{(2)}&S_{43}^{(2)}&0&S_{44}^{(2)}end{bmatrix}}}То есть, на каждом следующем шаге необходимо добавлять новые элементы к уже существующим.
Второй вид обобщения на несколько треугольников
Пусть есть область, представленная и разбитая на треугольники так, как представлено на рисунке. Пусть данная сетка содержит N {displaystyle N} узлов. Создадим глобальную матрицу S {displaystyle {mathfrak {S}}} (размера, очевидно, N × N {displaystyle N imes N} ) и заполним её нулями. Начнём строить локальные матрицы S {displaystyle S} для треугольников, например, для Δ 036. {displaystyle Delta 036.}
Введём локальную нумерацию для данного треугольника: пусть его верхняя вершина имеет локальный номер 0 {displaystyle 0} , далее по часовой стрелке 1 {displaystyle 1} и 2 {displaystyle 2} . Иначе говоря, пусть глобальным номерам 0 , 3 , 6 {displaystyle 0,3,6} соответствуют локальные номера 0 , 1 , 2 {displaystyle 0,1,2} соответственно.
Составим матрицу для этого треугольника так, как описано выше, получив что-то типа
S = [ S 00 S 01 S 02 S 10 S 11 S 12 S 20 S 21 S 22 ] {displaystyle S={egin{bmatrix}S_{00}&S_{01}&S_{02}S_{10}&S_{11}&S_{12}S_{20}&S_{21}&S_{22}end{bmatrix}}}Теперь заменим локальную нумерацию на глобальную. То есть запишем локальное число S 00 {displaystyle S_{00}} как глобальное число S 00 {displaystyle {mathfrak {S_{00}}}} , S 01 {displaystyle S_{01}} - как S 03 {displaystyle {mathfrak {S_{03}}}} , S 02 {displaystyle S_{02}} - как S 06 {displaystyle {mathfrak {S_{06}}}} и так далее.
Получим
S = [ 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 ] {displaystyle {mathfrak {S}}={egin{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&0S_{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&0S_{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&0end{bmatrix}}}С остальными треугольниками поступаем аналогично. Необходимо помнить, что надо "дописать" число в глобальную ячейку, то есть прибавить к уже существующему.
Учёт граничных условий
Условия Дирихле
В случае граничных условий первого рода необходимо изменить матрицу S {displaystyle {mathfrak {S}}} .
Граничное условие гласит, что функция в узлах на границе равна нулю. Для узла u i , j {displaystyle u_{i,j}} необходимо вычеркнуть i {displaystyle i} -й столбей и j {displaystyle j} -ю строку в матрице S {displaystyle {mathfrak {S}}} , а также вычеркнуть сам узел из массива узлов решётки.
Условия Нейнмана
В случае граничных условий второго рода глобальная матрица не меняется.