Stiffness matrix for a general bilinear element.

Summary

This notebook describes the computational steps required in the computation of the displacement based finite element stiffness matrix using an isoparametric formulation. In order to show the detailed steps we first describe the computations at a single integration point. These steps are later written in condensed form using available SOLIDSPy subroutines.

Stiffness matrix

The stiffness matrix in the displacement based finite elment method is defined like:

$${K^{QP}} = \int\limits_V {B_{ij}^Q{C_{ijkl}}B_{kl}^PdV} $$

where:

$$C = \frac{{E(1 - \nu )}}{{(1 + \nu )(1 - 2\nu )}}\left[ {\begin{array}{*{20}{c}} 1&{\frac{\nu }{{1 - \nu }}}&0\\ {\frac{\nu }{{1 - \nu }}}&1&0\\ 0&0&{\frac{{1 - 2\nu }}{{2(1 - \nu )}}} \end{array}} \right]$$

is the elastic tensor and $B_{ij}^Q$ is the strain-displacement interpolator associated to the $Q$-th degree of freedom.

Mathematical aspects

The isoparametric formulation can be easily understood if one assumes that the physical and the natural spaces are regarded as the result of deforming the element (in the physical space) into a perfect square (in the natural space). This pseudo-deformation process is ilustrated in the figure below:

and where $t$ is a pseudo-time and ${\vec r}$ and ${\vec x}$ denote the natural and physical space respectively. Mathematically, the relation between these spaces can be written like:

$${x_i} = {x_i}(r)$$$${r_I} = {r_I}(x).$$

Relating differential elements in both spaces we have:

$$d{x_i} = \frac{{\partial {x_i}}}{{\partial {r_J}}}d{r_J} \equiv {J_{iJ}}d{r_J}.$$

where $J_{iJ}$ is the Jacobian of the transformation. We can establish the relation between the gradient operators in both spaces like:

$$\frac{\partial }{{\partial {x_i}}} = \frac{{\partial {r_J}}}{{\partial {x_i}}}\frac{\partial }{{\partial {r_J}}}$$

or in compact notation: $$\nabla _i^x = J_{iJ}^{ - 1}\nabla _J^r$$

and where

$$J_{iJ}^{ - 1} = \frac{{\partial {r_J}}}{{\partial {x_i}}}.$$

Finite element approach

In the isoparametric finite element implementation the space variable (i.e., $x_i$) is approximatted via interpolation using the same shape functions $N_i^K(r)$ as in the approximation of the primary variable $u_i$:

$${x_i} = N_i^Q(r){{\hat x}^Q}.$$

In the above ${\hat x}^Q$ corresponds to the physical coordinates of node $Q$. The contribution to the numerical-Jacobian from the Q-th degree of freedom can now be written like:

$${J_{iJ}} \equiv \frac{{\partial {x_i}}}{{\partial {r_J}}} = \frac{{\partial N_i^Q(r)}}{{\partial {r_J}}}{{\hat x}^Q}.$$

For the computation of the term $B_{ij}$ appearing in the stiffness matrix $K$, and relating nodal displacements to displacement gradients according to a particular definition of the strain tensor, we require the fundamental interpolator $L_{ij}^K$. This will allows us to obtain displacement gradients ${u_{i,j}}$ out of the nodal displacements ${\hat u}^K$.

Using the expression for the displacement vector ${u_i}$ in terms of nodal displacements and the definition of the gradient operator derived pteviously, we have for the displacement gradients:

$${u_{i,j}} = J_{jQ}^{ - 1}\nabla _Q^rN_i^K(r){{\hat u}^K} \equiv L_{ij}^K{{\hat u}^K}$$

from which we identify the fundamental interpolator $L_{ij}^K(r)$ according to :

$$L_{ij}^K(r) = J_{jQ}^{ - 1}N_{i,Q}^K(r).$$

The strain-displacement operator for a given strain definition can now be assembled using combinations of the fundamental interpolator. Also, notice that in the computation of $L_{ij}^K(r)$ we require the inverse of the Jacobian of the transformation. In the actual implementation we evaluate directly the Jacobian using the interpolated version of $x_i$ and perform the invesrion numerically. For our particular definition of the strain tensor the strain-displacements operator takes the form:

$$B_{ij}^Q = \frac{1}{2}\left( {L_{ij}^Q + L_{ji}^Q} \right).$$

Computational aspects

Let us define ${\hat x}^Q$ as the vector of nodal points coordinates and $E$ and $\nu$ as the elastic modulus and Poisson's ratio respectively for the current material.

Computations at a single integration point

Retrieve the Gauss points and weights in vectors $XW$ and $XP$ respectively;

and let us consider the natural coordinates at a single point:

The first step involves the computation of the Jacobian operator and its inverse at the Gauss point:

$$J = \left[ {\begin{array}{*{20}{c}} {\frac{{\partial x}}{{\partial r}}}&{\frac{{\partial y}}{{\partial r}}}\\ {\frac{{\partial x}}{{\partial s}}}&{\frac{{\partial s}}{{\partial s}}} \end{array}} \right] \equiv \left[ {\begin{array}{*{20}{c}} {}&{\frac{{\partial {N^Q}}}{{\partial r}}}&{}\\ {...}&{}&{...}\\ {}&{\frac{{\partial {N^Q}}}{{\partial s}}}&{} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} .&.\\ .&.\\ {{x^Q}}&{{y^Q}}\\ .&.\\ .&. \end{array}} \right\}$$

Denoting the terms of the gradient operator in the physical space like ${L_x^Q}$ and ${L_Y^Q}$ respectivly these can be obtained like:

$$\left\{ {\begin{array}{*{20}{c}} {L_x^Q}\\ {}\\ {L_y^Q} \end{array}} \right\} = {\left[ J \right]^{ - 1}}\left\{ {\begin{array}{*{20}{c}} {}&{\frac{{\partial {N^Q}}}{{\partial r}}}&{}\\ {...}&{}&{...}\\ {}&{\frac{{\partial {N^Q}}}{{\partial s}}}&{} \end{array}} \right\}$$

allowing us to assemble the strain-displacement operator $B$ and to compute the strain tensor at the Gauss point as:

$$\left\{ {\begin{array}{*{20}{c}} {\frac{{\partial u}}{{\partial x}}}\\ {\frac{{\partial v}}{{\partial y}}}\\ {\frac{{\partial u}}{{\partial y}} + \frac{{\partial v}}{{\partial x}}} \end{array}} \right\} = \left[ {\begin{array}{*{20}{c}} {}&{L_x^Q}&0&{}\\ {...}&0&{L_y^Q}&{...}\\ {}&{L_y^Q}&{L_x^Q}&{} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} .\\ .\\ .\\ {{u^Q}}\\ {{v^Q}}\\ .\\ .\\ . \end{array}} \right\}$$

The stiffness matrix at the current integration point is computed according to:

$$k = {B^T}(r,s)CB(r,s)\left\| {J(r,s)} \right\|\alpha $$

Complete loop over the Gauss points. UEL() subroutine.

We now compute the complete stiffness matrix after looping along all the Gauss points. In the actual implementation Python uses the strain-displacement matrix subroutine stdm4n() contained within the femutils module. In SOLIDSPy the resulting block of code is implemented in the form of a single subroutine termed UEL().

Final stiffness matrix