Stiffness matrix 3D

Let us proceed to computing the stiffness matrix for the three-dimensional case. In three-dimensional space the gradient has the form $\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)$, and the dot product of the gradient with itself is $\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$. Since the three-dimensional computational domain $M$ is partitioned into tetrahedral simplices, the part of the functional under study for a single tetrahedron with vertices

can be written as

The function $\upsilon(x, y, z)$​${} = \sum_{i=1}^N \upsilon_i(x, y, z)$. The trial function on the tetrahedron has the form $\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)$. Analogously to the two-dimensional case, the hat functions for the tetrahedron are linear

Let us compute the partial derivatives of the trial function

Note that the derivatives do not depend on $x$, $y$ and $z$ and are constants over the tetrahedron. Substituting (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}$) into (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.$$)

where $V_{\text{tet}}$ is the volume of the tetrahedron, which is computed by the formula

where $\Delta$ is the determinant of the matrix

Let us expand the squares and regroup the terms

The coefficients $b_k$, $c_k$ and $d_k$ for each vertex $k$ of the tetrahedron are computed through the minors of the determinant $\Delta$. For the vertex with index $k$ the coefficients have the form

where $\Delta_k^{(x)}$, $\Delta_k^{(y)}$ and $\Delta_k^{(z)}$ are the minors obtained by deleting the $k$-th row and the corresponding column ($x$, $y$ or $z$) from the matrix (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}.$).

We introduce the notation for the elements of the local stiffness matrix of the tetrahedron

Thus, the local stiffness matrix for a tetrahedral element has the form

The local stiffness matrix is symmetric, that is $k_{mn}^{\text{loc}}$​${} = k_{nm}^{\text{loc}}$. The global stiffness matrix $\mathbf{K}$ is obtained by summing the contributions from all tetrahedral elements of the mesh by the assembly method: the elements of the local matrices are added to the corresponding elements of the global matrix according to the global node numbering. The dimension of the global stiffness matrix is $N \times N$, where $N$ is the total number of mesh nodes.