Deriving the element stiffness matrix for 2D linear elasticity

Question

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! :)

Answer 1

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).

Answer 2

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$.