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

Перейдём к вычислению матрицы жёсткости для трёхмерного случая. В трёхмерном пространстве градиент имеет вид $\nabla \upsilon$ $= \left( \frac{\displaystyle \partial \upsilon}{\displaystyle \partial x}, \frac{\displaystyle \partial \upsilon}{\displaystyle \partial y}, \frac{\displaystyle \partial \upsilon}{\displaystyle \partial z} \right)$, а скалярное произведение градиента с самим собой равно $\nabla \upsilon \cdot \nabla \upsilon$ $= \left( \frac{\displaystyle \partial \upsilon}{\displaystyle \partial x} \right)^2$ $+ \left( \frac{\displaystyle \partial \upsilon}{\displaystyle \partial y} \right)^2$ $+ \left( \frac{\displaystyle \partial \upsilon}{\displaystyle \partial z} \right)^2$. Учитывая, что трёхмерная расчётная область $M$ разбита на симплексы-тетраэдры, исследуемую часть функционала для одного тетраэдра с вершинами

можно записать как

Функция $\upsilon(x, y, z)$ $= \sum_{i=1}^N \upsilon_i(x, y, z)$. Пробная функция на тетраэдре имеет вид $\upsilon_{(i)(i+3)}(x, y, z)$ $= q_i \cdot \phi_i(x, y, z)$ $+ q_{i+1} \cdot \phi_{i+1}(x, y, z)$ $+ q_{i+2} \cdot \phi_{i+2}(x, y, z)$ $+ q_{i+3} \cdot \phi_{i+3}(x, y, z)$. Аналогично двумерному случаю, функции «крышек» для тетраэдра имеют линейный вид

Вычислим частные производные пробной функции

Заметим, что производные не зависят от $x$, $y$ и $z$ и являются константами на тетраэдре. Подставим (6.14$\displaystyle \frac{\displaystyle \partial \upsilon_{(i)(i+3)}(x, y, z)}{\displaystyle \partial x}$ $\displaystyle = q_i \cdot b_i$ $\displaystyle + q_{i+1} \cdot b_{i+1}$ $\displaystyle + q_{i+2} \cdot b_{i+2}$ $\displaystyle + q_{i+3} \cdot b_{i+3}$$\displaystyle \frac{\displaystyle \partial \upsilon_{(i)(i+3)}(x, y, z)}{\displaystyle \partial y}$ $\displaystyle = q_i \cdot c_i$ $\displaystyle + q_{i+1} \cdot c_{i+1}$ $\displaystyle + q_{i+2} \cdot c_{i+2}$ $\displaystyle + q_{i+3} \cdot c_{i+3}$$\displaystyle \frac{\displaystyle \partial \upsilon_{(i)(i+3)}(x, y, z)}{\displaystyle \partial z}$ $\displaystyle = q_i \cdot d_i$ $\displaystyle + q_{i+1} \cdot d_{i+1}$ $\displaystyle + q_{i+2} \cdot d_{i+2}$ $\displaystyle + q_{i+3} \cdot d_{i+3}$) в (6.12$$\int_{\text{тет}} \left[ \left( \frac{\displaystyle \partial \upsilon}{\displaystyle \partial x} \right)^2 + \left( \frac{\displaystyle \partial \upsilon}{\displaystyle \partial y} \right)^2 + \left( \frac{\displaystyle \partial \upsilon}{\displaystyle \partial z} \right)^2 \right] \,dV.$$)

где $V_{\text{тет}}$ — объём тетраэдра, который вычисляется по формуле

где $\Delta$ — определитель матрицы

Раскроем квадраты и перегруппируем члены

Коэффициенты $b_k$, $c_k$ и $d_k$ для каждой вершины $k$ тетраэдра вычисляются через миноры определителя $\Delta$. Для вершины с индексом $k$ коэффициенты имеют вид

где $\Delta_k^{(x)}$, $\Delta_k^{(y)}$ и $\Delta_k^{(z)}$ — миноры, получаемые вычёркиванием $k$-й строки и соответствующего столбца ($x$, $y$ или $z$) из матрицы (6.16$\displaystyle \Delta$ $\displaystyle = \begin{vmatrix} x_i & y_i & z_i & 1\\ x_{i+1} & y_{i+1} & z_{i+1} & 1\\ x_{i+2} & y_{i+2} & z_{i+2} & 1\\ x_{i+3} & y_{i+3} & z_{i+3} & 1 \end{vmatrix}.$).

Введём обозначения для элементов локальной матрицы жёсткости тетраэдра

Таким образом, локальная матрица жёсткости для тетраэдрального элемента имеет вид

Локальная матрица жёсткости является симметричной, то есть $k_{mn}^{\text{лок}}$ $= k_{nm}^{\text{лок}}$. Глобальная матрица жёсткости $\mathbf{K}$ получается путём суммирования вкладов от всех тетраэдральных элементов сетки методом сборки: элементы локальных матриц добавляются к соответствующим элементам глобальной матрицы согласно глобальной нумерации узлов. Размерность глобальной матрицы жёсткости равна $N \times N$, где $N$ — общее количество узлов сетки.