\(\newcommand{\beps}{\boldsymbol\varepsilon}\) \(\newcommand{\bsig}{\boldsymbol\sigma}\) \(\newcommand{\ud}{\mathrm{d}}\) \(\newcommand{\us}{\mathrm{s}}\) \(\newcommand{\ba}{\mathbf{a}}\) \(\newcommand{\bb}{\mathbf{b}}\) \(\newcommand{\bc}{\mathbf{c}}\) \(\newcommand{\bt}{\mathbf{t}}\) \(\newcommand{\bu}{\mathbf{u}}\) \(\newcommand{\bx}{\mathbf{x}}\) \(\newcommand{\bw}{\mathbf{w}}\) \(\newcommand{\bN}{\mathbf{N}}\) \(\newcommand{\bB}{\mathbf{B}}\) \(\newcommand{\bD}{\mathbf{D}}\) \(\newcommand{\bK}{\mathbf{K}}\) \(\newcommand{\pder}[2]{\frac{\partial #1}{\partial #2}}\) \(\newcommand{\iD}{\boldsymbol{\mathcal{D}}}\) \(\newcommand{\mbf}[1]{\mathbf{#1}}\) \(\newcommand{\mrm}[1]{\mathrm{#1}}\) \(\newcommand{\bs}[1]{\boldsymbol{#1}}\) \(\newcommand{\T}{^\mathrm{T}}\) \(\newcommand{\inv}{^{-1}}\) \(\newcommand{\myVec}[1]{\left\{ \begin{matrix} #1 \end{matrix} \right\}}\) \(\newcommand{\myMat}[1]{\left[ \begin{matrix} #1 \end{matrix} \right]}\) \(\newcommand{cA}[1]{\textcolor[RGB]{1,113,136}{#1}}\) \(\newcommand{cB}[1]{\textcolor[RGB]{195,49,47}{#1}}\) \(\newcommand{cC}[1]{\textcolor[RGB]{0,102,162}{#1}}\) \(\newcommand{cD}[1]{\textcolor[RGB]{0,183,211}{#1}}\) \(\newcommand{cE}[1]{\textcolor[RGB]{0,163,144}{#1}}\) \(\newcommand{cF}[1]{\textcolor[RGB]{97,164,180}{#1}}\) \(\newcommand{cG}[1]{\textcolor[RGB]{130,215,198}{#1}}\) \(\newcommand{cH}[1]{\textcolor[RGB]{153,210,140}{#1}}\) \(\newcommand{cI}[1]{\textcolor[RGB]{235,114,70}{#1}}\) \(\newcommand{cJ}[1]{\textcolor[RGB]{241,190,62}{#1}}\) \(\newcommand{cK}[1]{\textcolor[RGB]{231,41,138}{#1}}\)

3.1. Strong form for linear elasticity

In this chapter we treat mechanical equilibrium problems in a more general way. In the previous chapter we have seen how the 1D Poisson equation can be directly used to treat structural bar elements in extension. Here we give the equilibrium problem a more thorough treatment, starting from continuum mechanics and arriving at a strong form PDE we can treat with FEM.

Kinematic relations

For equilibrium problems, we are concerned with modeling how solids deform as loads are applied. Consider the following domain \(\Omega\) of arbitrary shape:

Fig. 3.1 Deformation of an arbitrary two-dimensional domain \(\Omega\).

where the solid is subjected to body loads \(\mbf{b}\), prescribed tractions \(\mbf{t}\) at surface \(\Gamma_N\) and prescribed displacements at surface \(\Gamma_D\). Assume the body moves as the loads are applied, and we represent this movement by the displacement vector \(\bu\) shown above. In the more general 3D case, we can decompose \(\bu\) into three components along our coordinate system axes, and we are therefore interested in solving for a vector field:

(3.1)\[\begin{split} \bu(\bx) \equiv \myVec{u_x(x,y,z)\\u_y(x,y,z)\\u_z(x,y,z)} \end{split}\]

Part of this displacement does not generate deformations, and are therefore called rigid-body motions. We are actually more interested in how the body actually deforms. For this we need to adopt a measure of strain. In this chapter we will be confining ourselves to small displacements and small strains. Under these assumptions, we can define strains as a tensor:

(3.2)\[ \varepsilon_{ij} = \displaystyle\frac{1}{2}\left(u_{i,j}+u_{j,i}\right) \]

where \(i,j\in\{1,2,3\}\). In matrix/vector notation this can also be written as:

(3.3)\[ \beps = \nabla^\mrm{s}\bu = \displaystyle\frac{1}{2}\left(\nabla\bu+\left(\nabla\bu\right)\T\right) \]

where \(\nabla^\mrm{s}\) is the symmetric gradient operator. We can also look at the strain tensor component by component:

(3.4)\[\begin{split} \beps = \myMat{ \varepsilon_{xx} & \varepsilon_{xy} & \varepsilon_{xz}\\ \varepsilon_{yx} & \varepsilon_{yy} & \varepsilon_{yz}\\ \varepsilon_{zx} & \varepsilon_{zy} & \varepsilon_{zz}\\ } = \myMat{ u_{x,x} & \frac{1}{2}(u_{x,y}+u_{y,x}) & \frac{1}{2}(u_{x,z}+u_{z,x})\\ \frac{1}{2}(u_{y,x}+u_{x,y}) & u_{y,y} & \frac{1}{2}(u_{y,z}+u_{z,y})\\ \frac{1}{2}(u_{z,x}+u_{x,z}) & \frac{1}{2}(u_{z,y}+u_{y,z}) & u_{z,z}\\ } \end{split}\]

from where we can see that \(\beps\) is symmetric, as expected. We can also switch to Voigt notation and represent strains as a vector instead:

(3.5)\[\begin{split} \beps = \myVec{\varepsilon_{xx}\\\varepsilon_{yy}\\\varepsilon_{zz}\\2\varepsilon_{xy}\\2\varepsilon_{yz}\\2\varepsilon_{zx}} = \myVec{\varepsilon_{xx}\\\varepsilon_{yy}\\\varepsilon_{zz}\\\gamma_{xy}\\\gamma_{yz}\\\gamma_{zx}} = \myVec{u_{x,x}\\u_{y,y}\\u_{z,z}\\u_{x,y}+u_{y,x}\\u_{y,z}+u_{z,y}\\u_{z,x}+u_{x,z}} \end{split}\]

Constitutive relations

Stress is a measure of internal forces inside a material, and is intrinsically linked to strains. Stresses are represented by a tensor \(\bsig\) and can be related to surface tractions through:

(3.6)\[ \mbf{t} = \bsig\cdot\mbf{n} \]

which are in turn related related to loads the body \(\Omega\) is subjected to. In the expression above, \(\mbf{n}\) is the unit normal vector to the surface on which the tractions act. The stress tensor has components:

(3.7)\[\begin{split} \bsig = \myMat{ \sigma_{xx} & \sigma_{xy} & \sigma_{xz}\\ \sigma_{yx} & \sigma_{yy} & \sigma_{yz}\\ \sigma_{zx} & \sigma_{zy} & \sigma_{zz}\\ } \end{split}\]

where the first index of each component associates it to a specific surface in an infinitesimal cubic material volume, and the second index indicates the direction of the stress.

For linear-elastic materials, there is a constant mapping between stresses and strains that depends on the Young’s modulus \(E\) and the Poisson’s ratio \(\nu\). Working with full strain and stress tensors, this mapping would be a fourth-order tensor:

(3.8)\[ \bsig = \bs{\mathcal{C}}:\beps \]

which is known as Hooke’s Law. The tensor \(\bs{\mathcal{C}}\) can be written in index notation as:

(3.9)\[ \mathcal{C}_{ijkl} = \mu(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}) + \lambda\delta_{ij}\delta_{kl} \]

where \(\lambda\) and \(\mu\) are the so-called Lamé constants:

(3.10)\[ \lambda=\displaystyle\frac{\nu E}{(1+\nu)(1-2\nu)} \quad\quad \mu = \displaystyle\frac{E}{2(1+\nu)} \]

Working in Voigt notation and exploiting symmetries in \(\bs{\mathcal{C}}\), we can reach a matrix-vector version of Hooke’s Law in three dimensions:

(3.11)\[\begin{split} \myVec{\sigma_{xx}\\\sigma_{yy}\\\sigma_{zz}\\\sigma_{xy}\\\sigma_{yz}\\\sigma_{zx}} = \frac{E}{(1+\nu)(1-2\nu)} \myMat{1-\nu & \nu & \nu & 0 & 0 & 0\\ \nu & 1-\nu & \nu & 0 & 0 & 0\\ \nu & \nu & 1-\nu & 0 & 0 & 0\\ 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 & 0\\ 0 & 0 & 0 & 0 & \frac{1-2\nu}{2} & 0\\ 0 & 0 & 0 & 0 & 0 & \frac{1-2\nu}{2}\\ } \myVec{\varepsilon_{xx}\\\varepsilon_{yy}\\\varepsilon_{zz}\\\gamma_{xy}\\\gamma_{yz}\\\gamma_{zx}} \end{split}\]

which we can write more compactly as:

(3.12)\[ \bsig = \mbf{D}\beps \]

When working in two dimensions (i.e. in the \(x-y\) plane), we need an extra assumption of how the material behaves out of plane. There are two main possibilities:

Plane strain assumption

Assume that all out-of-plane strain components are zero:

\[ \varepsilon_{zz} = \gamma_{yz} = \gamma_{zx} = 0 \]

Hooke’s Law reduces to:

\[\begin{split} \myVec{\sigma_{xx}\\\sigma_{yy}\\\sigma_{xy}} = \frac{E}{(1+\nu)(1-2\nu)} \myMat{1-\nu & \nu & 0\\ \nu & 1-\nu & 0\\ 0 & 0 & \frac{1-2\nu}{2}} \myVec{\varepsilon_{xx}\\\varepsilon_{yy}\\\gamma_{xy}} \end{split}\]

Plane stress assumption

Assume that the all out-of-plane stress components are zero:

\[ \sigma_{zz} = \sigma_{yz} = \sigma_{zx} = 0 \]

Hooke’s Law reduces to:

\[\begin{split} \myVec{\sigma_{xx}\\\sigma_{yy}\\\sigma_{xy}} = \frac{E}{1-\nu^2} \myMat{1 & \nu & 0\\ \nu & 1 & 0\\ 0 & 0 & \frac{1-\nu}{2}} \myVec{\varepsilon_{xx}\\\varepsilon_{yy}\\\gamma_{xy}} \end{split}\]

Equilibrium equations

To get to a final strong-form PDE we need to enforce stress equilibrium over our domain \(\Omega\). We can arrive at a result in two ways. First, consider the small volume of material \(\Delta\Omega\) shown below:

Fig. 3.2 Stresses acting on a small piece of continuum \(\Omega\).

where we use only two dimensions for simplicity. Translational equilibrium in \(x\) and \(y\) directions within this domain therefore requires:

(3.13)\[ \displaystyle\frac{\partial\sigma_{xx}}{\partial x} + \frac{\partial\sigma_{yx}}{\partial y} + b_x = 0 \quad\quad \displaystyle\frac{\partial\sigma_{yy}}{\partial y} + \frac{\partial\sigma_{xy}}{\partial x} + b_y = 0 \]

These equilibrium relations can be summarized in a single expression as:

(3.14)\[ \nabla\cdot\bsig+\mbf{b}=\mbf{0} \]

For the rotational equilibrium the easiest is to compute the moments at the center of the volume. This leads to:

(3.15)\[ \sigma_{xy} = \sigma_{yx} \]

or, in other words \(\bsig=\bsig\T\), which means the stress tensor should be symmetric.

Another way to arrive at the same result is to consider the equilibrium of the whole body \(\Omega\) in an integral sense:

(3.16)\[ \displaystyle\int_\Gamma\bsig\cdot\mbf{n}\,\ud\Gamma + \int_\Omega\mbf{b}\,\ud\Omega = 0 \]

With the divergence theorem we can turn the integral in \(\Gamma\) into an integral in \(\Omega\). Then, putting both terms together we have:

(3.17)\[ \displaystyle\int_\Omega\nabla\cdot\bsig + \mbf{b}\,\ud\Omega = \mbf{0} \]

and since this must hold also not only for the whole body \(\Omega\), but also for every subdomain inside the body, it follows from this expression that the earlier obtained Equation (3.14) must hold at every point:

(3.18)\[ \nabla\cdot\bsig+\mbf{b}=\mbf{0} \quad\text{in}\,\,\Omega \]

Enforcing rotational equilibrium in an integral sense is a bit more involved, but results again in the requirement that \(\bsig\) is symmetric.