Introduction
This article is part of my series on basic FEA theory. The series started with two Python code examples (Full Finite Element Solver in 100 Lines of Python and Full Finite Element Solver in 200 Lines of Python), and articles on (1) virtual work, (2) isoparametric shape functions, and (3) strain calculations.
This article focuses on how to calculate the stress from the nodal displacements, and how to calculate the element stiffness matrix.
Calculate Stresses from Strains (and Nodal Displacements)
For a linear elastic material the stress can be calculated from the strain components using Hooke’s law: \[ \sigma_{ij} = 2 \mu \varepsilon_{ij} + \lambda \varepsilon_{kk} \delta_{ij},\] which for plane strain can be written \[\begin{cases} \sigma_{11} &= 2\mu \varepsilon_{11} + \lambda(\varepsilon_{11} + \varepsilon_{22})\\ \sigma_{22} &= 2\mu \varepsilon_{22} + \lambda(\varepsilon_{11} + \varepsilon_{22})\\ \sigma_{12} &= 2\mu \varepsilon_{12}.\end{cases}\] Most of the time we are more interested in using the Young’s modulus and Poisson’s ratio than the 2 Lamé constants (\(\mu,\lambda\)). We can easily convert the equation using the following table of elastic constants. For plane strain the stress can then be obtained from: \[ \boldsymbol{\sigma}^{(m)} = \mathbf{C}^{(m)} \boldsymbol{\varepsilon}^{(m)} = \mathbf{C}^{(m)} \mathbf{B}^{(m)} \mathbf{U}^{(m)},\] where:
\[ \mathbf{C}^{(m)}_{ij} = \displaystyle\frac{E}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1-\nu & \nu & 0 \\ \nu & 1-\nu & 0\\ 0 & 0 & 0.5-\nu \end{bmatrix}.\]
Element Stiffness Matrix
The first term in the Principle of Virtual Work can be written:
\[ \displaystyle\sum_m \int_{V^{(m)}} \delta\boldsymbol{\varepsilon}^\top \boldsymbol{\sigma} dV.\]
This equation is difficult to calculate since it needs to be integrated over an irregular shape.
To simplify the calculation we first change variables from \((x,y)\) to \((\xi,\eta)\), so the integral becomes of the form: \[\int_A f(x,y) dx dy = \int_A f(\xi,\eta) J d\xi d\eta,\] where we know the area factor \(J = \det[d\mathbf{N} \mathbf{x}]\) from the shape functions. To evaluate the integral in the Principle of Virtual Work we will use Gaussian Quadrature.
Recall that a simple 1D integral can be numerically approximated as follows: \[ \int_{-1}^1 f(x) dx \approx \displaystyle\sum_{i=1}^n w_i f(x_i).\] If \(n=1\) then \(x_1=0\) and \(w_1=2\), and if \(n=2\) then \(x_i=\pm 1/\sqrt{3}\) and \(w_i=1\).
\[\delta\mathbf{U}^\top \left[ \int_V \mathbf{B}^\top \mathbf{C} \mathbf{B} dV \right] \mathbf{U} = \delta\mathbf{U}^\top \mathbf{k}^{(m)} \mathbf{U} = \delta\mathbf{U}^\top \left[ \displaystyle\sum_{p}^{np} \mathbf{B}^\top \mathbf{C} \mathbf{B} J w_p \right]\mathbf{U}.\]
In the next article in this series I will show how to assemble the stiffness matrices and solve for the solution!
