Linear Problem of Elasticity Theory

In this chapter, we will use matrix notation. Tensor quantities such as displacement, strain, and stress will be represented as column vectors.

Displacement vector

$$\u = \vr{u_1}{u_2}{u_3} \rightarrow \u = \vr{u_1}{u_2}{u_3}^T$$

The strain state is represented by a matrix of components. The strain tensor is symmetric; therefore, using matrix notation, it can be written as a vector with six components.

$$\ee = \mat {\varepsilon_{11}}{\varepsilon_{12}}{\varepsilon_{13}}{\varepsilon_{21}}{\varepsilon_{22}}{\varepsilon_{23}}{\varepsilon_{31}}{\varepsilon_{32}}{\varepsilon_{33}} \rightarrow \ee = \begin{bmatrix} \varepsilon_{11} & \varepsilon_{22} & \varepsilon_{33} & \varepsilon_{12} & \varepsilon_{13} & \varepsilon_{23} \end{bmatrix}^T$$

Quite often, the engineering strain is used, so the strain is given by

$$\ee = \begin{bmatrix} \varepsilon_{11} & \varepsilon_{22} & \varepsilon_{33} & \gamma_{12} & \gamma_{13} & \gamma_{23} \end{bmatrix}^T$$

In the above formula, we use the so-called engineering definitions of shear strains (angles of shear deformation), which are related to the corresponding components of the strain tensor through the following relationships:

$$\gamma_{12} = 2\varepsilon_{12}, \quad \gamma_{13} = 2\varepsilon_{13}, \quad \gamma_{23} = 2\varepsilon_{23},$$

Similarly, the stress tensor can be written as:

$$\ss = \mat {\sigma_{11}}{\sigma_{12}}{\sigma_{13}}{\sigma_{21}}{\sigma_{22}}{\sigma_{23}}{\sigma_{31}}{\sigma_{32}}{\varepsilon_{33}} \rightarrow \ss = \begin{bmatrix} \sigma_{11} & \sigma_{22} & \sigma_{33} & \sigma_{12} & \sigma_{13} & \sigma_{23} \end{bmatrix}^T$$

General equations

The problem of linear elasticity theory is described by the following system of partial differential equations:

– Navier’s equilibrium equations, where $b$ denotes the body force vector.

$$\pp{\sigma_{ij}}{x_j} + b_i = 0, \quad i=1,2,3$$

• geometric equations, i.e. relationships between displacements and strains

$$\varepsilon_{ij} = \cfrac12 \left(\pp{u_i}{x_j} + \pp{u_j}{x_i} \right), \quad i=1,2,3, \quad j=1,2,3$$

• constitutive (physical) equations, which relate strain and stress, e.g. the isotropic Hooke’s law

$$\sigma_{ij} = \cfrac{E}{1+\nu}\varepsilon_{ij} + \cfrac{\nu E}{(1+\nu)(1-2\nu)}\varepsilon_{kk}\delta_{ij}$$

To the above equations, two types of boundary conditions are added: stress (traction) and displacement boundary conditions.

$$\sigma_{ij} n_j = t_i^{*} \in S_{\sigma}, \quad u_i = u_i^{*} \in S_u$$

The problem of linear elasticity can also be formulated using matrix notation by introducing the differential operator matrix $\L$ and the stiffness matrix $\D$.

$$\L = \begin{bmatrix} \pp{}{x_1} & 0 & 0 \\ 0 & \pp{}{x_2}& 0 \\ 0 & 0 & \pp{}{x_3} \\ \pp{}{x_2} & \pp{}{x_1} & 0\\ \pp{}{x_3} & 0 & \pp{}{x_1} \\ 0 & \pp{}{x_3} & \pp{}{x_2} \\ \end{bmatrix}, \quad \D = \cfrac{E}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1-\nu & \nu & \nu & 0 & 0& 0 \\ \nu & 1-\nu & \nu & 0 & 0& 0 \\ \nu & \nu & 1-\nu & 0 & 0& 0 \\ 0 & 0 & 0 & \cfrac{(1-2\nu)}{2} & 0 & 0 \\ 0 & 0 & 0 & 0 & \cfrac{(1-2\nu)}{2} & 0 \\ 0 & 0 & 0 & 0 & 0 & \cfrac{(1-2\nu)}{2} \\ \end{bmatrix}$$

Our equations then take the following form:

– Navier’s equilibrium equations (matrix form)

$$\L^T\ss + \b = \boldsymbol{0}$$ $$\begin{bmatrix} \cfrac{\partial}{\partial x_1} & 0 & 0 & \cfrac{\partial}{\partial x_2} & \cfrac{\partial}{\partial x_3} & 0 \\ 0 & \cfrac{\partial}{\partial x_2} & 0 & \cfrac{\partial}{\partial x_1} & 0 & \cfrac{\partial}{\partial x_3} \\ 0 & 0 & \cfrac{\partial}{\partial x_3} & 0 & \cfrac{\partial}{\partial x_1} & \cfrac{\partial}{\partial x_2} \end{bmatrix} \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{12} \\ \sigma_{13} \\ \sigma_{23} \end{bmatrix} + \begin{bmatrix} b_x \\ b_y \\ b_z \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}$$

• geometric equations (matrix form)

$$\quad \ee = \L\u$$

\[
\begin{bmatrix}
\varepsilon_{11} \\
\varepsilon_{22} \\
\varepsilon_{33} \\
\varepsilon_{12} \\
\varepsilon_{13} \\
\varepsilon_{23}
\end{bmatrix}
=
\begin{bmatrix}
\cfrac{\partial}{\partial x_1} & 0 & 0 \\
0 & \cfrac{\partial}{\partial x_2} & 0 \\
0 & 0 & \cfrac{\partial}{\partial x_3} \\
\cfrac{\partial}{\partial x_2} & \cfrac{\partial}{\partial x_1} & 0 \\
\cfrac{\partial}{\partial x_3} & 0 & \cfrac{\partial}{\partial x_1} \\
0 & \cfrac{\partial}{\partial x_3} & \cfrac{\partial}{\partial x_2}
\end{bmatrix}
\begin{bmatrix}
u_1 \\
u_2 \\
u_3
\end{bmatrix}
\]

• constitutive equations

$$\ss = \D \ee$$ $$\begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{12} \\ \sigma_{13} \\ \sigma_{23} \end{bmatrix} = \frac{E}{(1+\nu)(1-2\nu)} \begin{bmatrix} 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} \end{bmatrix} \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ \varepsilon_{12} \\ \varepsilon_{13} \\ \varepsilon_{23} \end{bmatrix}$$

Basics of FEM derived from the principle of minimum total potential energy

The total energy of a body is given by the formula:

$$E = \Pi + K,$$

where $\Pi$ is the potential energy and $K$ is the kinetic energy. For static problems, the kinetic energy of the body is neglected. The potential energy in this case can be expressed as

$$\Pi = U-W$$

where $U$ is the potential energy of internal forces (strain energy), and $W$ is the potential energy of external forces, equal to the negative of the work done by these forces on the displacements.

$$\Pi = \cfrac12 \int_{\Omega}\boldsymbol{\sigma}^T\boldsymbol{\varepsilon} dV - \int_{\Omega} \u^T \b dV - \int_{\partial \Omega} \u^T \t dS$$

The principle of minimum total potential energy states that:
Among all geometrically admissible displacement fields that a body can undergo, the actual displacement field is the one for which the total potential energy attains its minimum value

$$\Pi \rightarrow \Pi_{min}$$

Let us assume that the finite element solution is expressed as

$$\u = \N \boldsymbol{d}$$

where matrix $\N$ contains shape functions, and vector $\d$ represents the unknown generalized displacements.

The geometric equation, after applying the displacement field approximation, becomes

$$\ee = \L \u = \L\N\boldsymbol{d} = \B \boldsymbol{d}, \quad \B = \L\N.$$

Substituting into the energy functional, we obtain

$$\Pi= \cfrac12 \int_{\Omega}(\D\B\d)^T \B\d dV - \int_{\Omega} \d^T \N^T \b dV - \int_{\partial \Omega} \d^T \N^T \t dS$$ $$\Pi= \cfrac12 \int_{\Omega}\d^T \B^T \D \B \d dV - \int_{\Omega} \d^T \N^T \b dV - \int_{\partial \Omega} \d^T \N^T \t dS$$

From the stationarity of this expression follows the equilibrium of the element:

$$\pp{}{\boldsymbol{d}} \Pi = 0 \rightarrow \left( \int_{\Omega} \B^T \D \B dV\right)\d -\int_{\Omega}\N^T \b dV - \int_{\partial \Omega} \N^T \t dS$$

The above equation can be rewritten as

$$\k \d = \f$$

where the stiffness matrix $\k$ is defined as

$$\boldsymbol{k} = \int_{\Omega} \B^T \D \B dV$$

and the load vector is given by

$$\f = \int_{\Omega}\N^T \b dV + \int_{\partial \Omega} \N^T \t dS.$$

This is the general FEM stiffness matrix, valid for any continuum problem (1D, 2D, 3D solids).

Bar (truss) element

For a 1D axial bar element we assume:

• differential operator

$$\L = \left[\dd{}{x} \right]$$

• nodal displacement vector

$$\d = [u_1, u_2]^T$$

• shape functions (linear interpolation)

$$\N = \left[\cfrac{x_2-x}{L^e}, \quad \cfrac{x-x_1}{L^e} \right]$$

• strain–displacement relation

$$\varepsilon = \cfrac{du}{dx} \quad \B^e = \L\N^e = \left[\cfrac{-1}{L^e}, \cfrac{1}{L^e}\right]$$

• material law (Hooke’s law)

$$\D = \left[E\right]$$

Since the integration is performed over the volume, we must include the cross-sectional area of the bar before the integral. The volume of the element in the case of a bar can be written as

$$\int_{\Omega^e} dV = A\int_{x_1}^{x_2} dx$$

Bar Element Stiffness Matrix Derivation:

$$\mathbf{k} = \int_{\Omega^e} \mathbf{B}^T \mathbf{D} \mathbf{B} \, dV$$ $$\mathbf{k} = AE \int_{x_1}^{x_2} \mathbf{B}^T \mathbf{B} \, dx$$ $$\mathbf{k} = AE \int_{x_1}^{x_2} \frac{1}{(L^e)^2} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} dx$$ $$\k = \cfrac{EA}{L^e}\matt{1}{-1}{-1}{1}$$

Euler–Bernoulli Beam Element

For a 1D Euler–Bernoulli beam element we assume:

• differential operator

$$\mathbf{L} =\begin{bmatrix}\frac{d^2}{dx^2}\end{bmatrix}$$

• nodal displacement vector

$$\mathbf{d} =\begin{bmatrix}u_1 & u_2 & u_3 & u_4\end{bmatrix}^T$$

• shape functions (cubic Hermite interpolation)

$$\mathbf{N} =\begin{bmatrix}N_1 & N_2 & N_3 & N_4\end{bmatrix}$$

with:

$$\u = N_1 u_1 + N_2 u_2 + N_3 u_3 + N_4 u_4$$

• strain–displacement relation (curvature)

$$\varepsilon = \kappa = \frac{d^2 \u}{dx^2}\quad \Rightarrow \quad\mathbf{B} = \mathbf{L}\mathbf{N}$$

• material law (Hooke’s law for bending)

$$\mathbf{D} = [EI]$$

Since the integration is performed over the volume, we include the cross-sectional area contribution through the second moment of inertia \( I \):

$$\int_{\Omega^e} dV = A \int_{x_1}^{x_2} dx \quad \Rightarrow \quad \text{bending stiffness uses } EI$$

Beam Element Stiffness Matrix Derivation

$$\mathbf{k} =\int_{\Omega^e} \mathbf{B}^T \mathbf{D} \mathbf{B} \, dV$$ $$\mathbf{k} =\int_{x_1}^{x_2} \mathbf{B}^T EI \mathbf{B} \, dx$$

\[
\mathbf{k} =
EI \int_{x_1}^{x_2} \mathbf{B}^T \mathbf{B} \, dx
\]

After substituting Hermite shape functions and computing derivatives:

\[
\mathbf{B} =
\begin{bmatrix}
\frac{d^2 N_1}{dx^2} &
\frac{d^2 N_2}{dx^2} &
\frac{d^2 N_3}{dx^2} &
\frac{d^2 N_4}{dx^2}
\end{bmatrix}
\]

\[
\mathbf{k} =
\frac{EI}{(L^e)^3}
\begin{bmatrix}
12 & 6L^e & -12 & 6L^e \\
6L^e & 4(L^e)^2 & -6L^e & 2(L^e)^2 \\
-12 & -6L^e & 12 & -6L^e \\
6L^e & 2(L^e)^2 & -6L^e & 4(L^e)^2
\end{bmatrix}
\]