# Deriving the element stiffness matrix for 2D linear elasticity

I'm following the derivation from Finite Element Method using Matlab 2nd Edition, pg 311-315, which derives of the local stiffness matrix for planar isotropic linear elasticity as follows:

Force Balance Equations
$\frac{\partial\sigma_x}{\partial x}+\frac{\tau_{xy}}{\partial y} + f_x=0$
$\frac{\partial\tau_{xy}}{\partial x}+\frac{\sigma_y}{\partial y} + f_y=0$

Using galerkin method, we multiply the first and second equation by test functions $w_1$ and $w_2$, respectively and integrate over the domain $\Omega$. Integrating by parts, I see that we obtain the weak formulation:

$$\int_\Omega \begin{bmatrix} \frac{\partial w_1}{\partial x} & 0 & \frac{\partial w_1}{\partial y}\\ 0 & \frac{\partial w_2}{\partial y} & \frac{\partial w_2}{\partial x} \end{bmatrix} \begin{bmatrix} \sigma_x \\ \sigma_y \\ \tau_{xy} \end{bmatrix} = \int_\Omega \begin{bmatrix} w_1 f_x \\ w_2 f_y \end{bmatrix} + \int_{\partial\Gamma} \begin{bmatrix} w_1 \Phi_x \\ w_2 \Phi_y \end{bmatrix}$$ where $\Gamma$ is the portion of the boundary with the neumann (traction) boundary condition in the x and y directions $\Phi_x$ and $\Phi_y$.

Let's just consider the integral on the left hand side of this equation. Using the linear isotropic stress strain relationship in two dimensions we can rewrite this equation as

$$\int_\Omega M D \epsilon$$

where

$M=\begin{bmatrix} \frac{\partial w_1}{\partial x} & 0 & \frac{\partial w_1}{\partial y}\\ 0 & \frac{\partial w_2}{\partial y} & \frac{\partial w_2}{\partial x}\end{bmatrix}$, $D=\frac{E}{1-\nu^2} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2} \end{bmatrix}$, and $\epsilon = \begin{bmatrix} \frac{\partial u}{\partial x} \\ \frac{\partial v}{\partial y} \\ \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \end{bmatrix}$.

Suppose the domain is tesselated into triangular elements. Consider a single element $e$ and the three basis functions (unit hat functions) on this element as $H_i(x)$ for i=1,2,3. We can charaterize the displacement functions on this element as $u(x,y)=\sum_{i=1}^3 u_iH_i$ and $v(x,y)=\sum_{i=1}^3 v_iH_i$, then we can rewrite

$$\epsilon=Bd$$

where $B=\begin{bmatrix} \frac{\partial H_1}{\partial x} & 0 & \frac{\partial H_2}{\partial x} & 0 & \frac{\partial H_3}{\partial x} & 0 \\ 0 & \frac{\partial H_1}{\partial y} & 0 & \frac{\partial H_2}{\partial y} & 0 & \frac{\partial H_3}{\partial y} \\ \frac{\partial H_1}{\partial y} & \frac{\partial H_1}{\partial x} & \frac{\partial H_2}{\partial y} & \frac{\partial H_2}{\partial x} & \frac{\partial H_3}{\partial y} & \frac{\partial H_3}{\partial x}\end{bmatrix}$ and $d=\begin{bmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ u_3 \\ v_3 \end{bmatrix}$.

Thus, we can write the integral over each element as

$$\int_e M D \epsilon= \int_e MDBd$$.

The author claims that the matrix $M$ becomes $B^T$ when we only consider the test functions equivalent to basis functions with support on element $e$. That is, when $w_1,w_2 = H_i$ for $i=1,2,3$ we obtain

$$\int_e B^TDBd$$

It's not immediately obvious to me why the matrix $M$ becomes $B^T$ over the element $e$. How can I arrive at this conclusion just by letting the test functions be the basis functions? Any help with this would be greatly appreciated! :)

My answer will address plane strain (instead of plane stress) using a different notation. I use Lamé parameters instead of Young's modulus and Poisson ratio for notational convenience. Consider

$$\DeclareMathOperator{\trace}{tr} -\nabla\cdot \big[2 \mu D u + \lambda \trace(Du) I \big]$$

where $D u = \frac 1 2 \Big[ \nabla u + (\nabla u)^T \Big]$ is the strain tensor (frequently denoted by $\varepsilon$, but my notation will be cleaner here; note that the cross-term is scaled different from your $\epsilon$) written as a linear differential operator $D$ ("symmetric gradient") applied to the displacement vector $u$. Now multiply by test functions $v$, integrate by parts, and look at the interior term

$$\int_\Omega \nabla v : \Big[2 \mu D u + \lambda \trace(Du) I\Big] = \int_\Omega 2\mu \nabla v:D u + \lambda \trace(Du) \nabla v:I \\ = \int_\Omega 2\mu Dv : Du + \lambda \trace(Dv) \trace(Du)$$

where we have used the property that a tensor contraction of nonsymmetric tensor with a symmetric tensor $N : S$ is equivalent to the symmetrized form $\frac 1 2 (N + N^T) : S$ to replace instances of $\nabla v$ with $D v$.

This latter form is obviously symmetric and what I recommend computing with. To get back to your notation, choose your ordering convention for the three entries in the symmetric tensors $Du$ and $Dv$, discretize the test and trial functions $v_h$ and $u_h$ as well as the differential operator $D_h$, and expand the integral as

$$2 \mu v_h^T D_h^T \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 2\end{pmatrix} D_h u_h + \lambda v_h^T D_h^T \begin{pmatrix} 1 & 1 & 0 \\ 1 & 1 & 0 \\ 0 & 0 & 0 \end{pmatrix} D_h u_h$$

Differentiating yields the Jacobian

$$D_h^T \begin{pmatrix} 2\mu + \lambda & \lambda & 0 \\ \lambda & 2\mu + \lambda & 0 \\ 0 & 0 & 4 \mu \end{pmatrix} D_h$$

which is the stiffness matrix contribution associated with a quadrature point (single-point quadrature integrates this exactly for $P_1$ triangular elements).

• So then, your martix $D_h$ is equivalent to my B matrix if i choose the alternating ordering $u_1, v_1, u_2, v_2, u_3, v_3$ – Paul Nov 4 '12 at 18:41
• Almost, but the last row (which produces the off-diagonal entry of the strain tensor) is scaled differently. Your notation sums the two entries which produces (in my opinion) confusing factors of 2. (I also hate that the second Lamé parameter $\mu$ is defined with the factor of 2, but sadly, that convention is too ingrained to change.) I never defined the ordering of basis functions (columns of $D_h$) as that choice doesn't matter for this discussion (anything is okay, you just have to be consistent). – Jed Brown Nov 4 '12 at 19:44
• In continuum mechanics there are two notations for stress/strain: tensor and engineering. Jed is using tensor notation: $\varepsilon_{ij} = \frac12 (u_{i,j}+u_{j,i})$. Paul is using engineering notation (which is commonly adopted in FEM codes): $(\varepsilon_x, \varepsilon_y, \gamma_{xy}) = (\varepsilon_{11}, \varepsilon_{22}, 2\varepsilon_{12})$. The factor 2 accounts for the fact that virtual work (contraction) has to be the same in both notations: $\sigma_{ij}\varepsilon_{ij} = \sigma_x\varepsilon_x + \sigma_y\varepsilon_y + \tau_{xy}\gamma_{xy}$. (Einstein summation for tensor notation). – Stefano M Nov 7 '12 at 12:10

The derivation presented by the OP in the question is, in my opinion, if not wrong, at least highly confusing. In fact after introducing test functions and a weak formulation, it fails to sum the two residuum equations obtained.

If we define the force-balance residuum components as \begin{align*} r_x &= \frac{\partial\sigma_x}{\partial x}+\frac{\tau_{xy}}{\partial y} + f_x, \\ r_y &= \frac{\partial\tau_{xy}}{\partial x}+\frac{\sigma_y}{\partial y} + f_y, \end{align*} we note that the weak formulation is written in the question above as two distinct scalar equations: \begin{equation} \int_\Omega \begin{bmatrix} w_1\,r_x \\ w_2 \,ry\end{bmatrix}\,d\Omega = \begin{bmatrix} 0 \\ 0 \end{bmatrix} \end{equation} or more simply \begin{align*} \int_\Omega w_1\,r_x \: d\Omega &= 0, & \int_\Omega w_2\,r_y \: d\Omega &= 0. \end{align*} Since we are interested in obtaining a formulation in which an energy norm is defined, the weak formulation should be written instead as a single scalar equation: \begin{equation} \int_\Omega w_1\,r_x \: d\Omega + \int_\Omega w_2\,r_y \: d\Omega = 0 \end{equation} or in vector notation \begin{equation} \int_\Omega \begin{bmatrix} w_1 \\ w_2 \end{bmatrix}^T \begin{bmatrix} r_x \\ r_y \end{bmatrix} \: d\Omega = 0. \end{equation} It is easy to show now that after integration by parts we have \begin{equation} -\int_\Omega \begin{bmatrix} \frac{\partial w_1}{\partial x} \\ \frac{\partial w_2}{\partial y} \\ \frac{\partial w_1}{\partial y} + \frac{\partial w_2}{\partial x} \end{bmatrix}^T \; D \; \begin{bmatrix} \frac{\partial u}{\partial x} \\ \frac{\partial v}{\partial y} \\ \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \end{bmatrix} \: d\Omega + \dots = 0 \end{equation} From this symmetric expression ($D=D^T$ and positive definite) it is clear why the expression of the element stiffness matrix results in $\int_e B^TDB \,d\Omega$ if both $(w_1, w_2)$ and $(u, v)$ are discretized by the same basis (shape) functions $H_i$.

Please note also that \begin{equation} a(w,u) = \int_\Omega \begin{bmatrix} \frac{\partial w_1}{\partial x} \\ \frac{\partial w_2}{\partial y} \\ \frac{\partial w_1}{\partial y} + \frac{\partial w_2}{\partial x} \end{bmatrix}^T \; D \; \begin{bmatrix} \frac{\partial u}{\partial x} \\ \frac{\partial v}{\partial y} \\ \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \end{bmatrix} \: d\Omega \end{equation} is the inner product in the strain energy norm, so that $\frac12 a(u,u)$ is the strain energy associated with displacement field $u$.

• Actually, this is the notation used in the book I cited above. – Paul Nov 11 '12 at 20:28
• @Paul I understand you obtained the notation from the cited book: unfortunately I have no way to check the reference, so I prefer to make statements only on what is written on this thread. No blame is intended on you, for sure, since it seems that you were actually confused from what seems a poor explanation of FEM in linear elasticity. Note that the trick is simply summing the two rows of matrix $M$, and everything becomes symmetric: no need actually to introduce tensor notation, as in the (however correct and fine) answer by Jed. – Stefano M Nov 11 '12 at 23:45