Element Stiffness Matrix

The element stiffness matrix in its matrix formulation is given by

\[\boldsymbol{K}^{e} = \int \limits_{\Omega^{e}} {\boldsymbol{B}^{e}}^{T} \boldsymbol{C}^{e} \boldsymbol{B}^{e} \,\mathrm{d} \Omega.\]

After transforming the integral from real space to reference space the expression becomes

\[\boldsymbol{K}^{e} = \int \limits_{-1}^{1} \int \limits_{-1}^{1} {\boldsymbol{B}^{e}}^{T} \boldsymbol{C}^{e} \boldsymbol{B}^{e} \det\bigl(\boldsymbol{J}^{e}\bigr) \,\mathrm{d} \xi \,\mathrm{d} \eta.\]

In many cases this integral can not be solved analytically. There are, however, a handful of ways to approximate it numerically. The so-called Gauss quadrature is a widely employed method in finite element software for the approximation of this integral.

Number of points, \(n\) Points, \(x_{i}\) Weights, \(w_{i}\)
\(1\) \(0\) \(2\)
\(2\) \(\pm \sqrt{\dfrac{1}{3}}\) \(1\)
\(3\) \(0\) \(\dfrac{8}{9}\)
\(\pm \sqrt{\dfrac{3}{5}}\) \(\dfrac{5}{9}\)

In our case we are going to use two Gauss points with their corresponding weights.

\[\boldsymbol{K}^{\mathrm{e}} = \sum \limits_{i = 1}^{2} \sum \limits_{j = 1}^{2} {\boldsymbol{B}^{\mathrm{e}}}^{T}(\xi_{i}, \eta_{j}) \boldsymbol{C}^{\mathrm{e}} \boldsymbol{B}^{\mathrm{e}}(\xi_{i}, \eta_{j}) \det\bigl(\boldsymbol{J}^{\mathrm{e}}(\xi_{i}, \eta_{j})\bigr) w_{j} w_{i}\]

And now we can implement a routine to give us the element stiffness matrix.

Element.stiffness_matrix()

Return the stiffness matrix of the element.

Returns:The stiffness matrix of the element.
Return type:numpy.ndarray