Applications in Solid Mechanics

Both are two-dimensional idealisations of a 3D elasticity problem, which reduce the unknowns from six stress components to three.

Comparison

Constitutive matrices

Plane stress

$$\begin{Bmatrix} \sigma_x \\ \sigma_y \\ \tau_{xy} \end{Bmatrix} = \frac{E}{1 - \nu^2}\begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1 - \nu}{2} \end{bmatrix}\begin{Bmatrix} \varepsilon_x \\ \varepsilon_y \\ \gamma_{xy} \end{Bmatrix}$$

Plane strain

$$\begin{Bmatrix} \sigma_x \\ \sigma_y \\ \tau_{xy} \end{Bmatrix} = \frac{E}{(1 + \nu)(1 - 2\nu)}\begin{bmatrix} 1 - \nu & \nu & 0 \\ \nu & 1 - \nu & 0 \\ 0 & 0 & \frac{1 - 2\nu}{2} \end{bmatrix}\begin{Bmatrix} \varepsilon_x \\ \varepsilon_y \\ \gamma_{xy} \end{Bmatrix}$$

The plane strain matrix is obtained from the plane stress matrix by replacing $E \to \dfrac{E}{1 - \nu^2}$ and $\nu \to \dfrac{\nu}{1 - \nu}$. Plane strain is stiffer for the same $E$ and $\nu$, because the z-direction is constrained.

(a) Plane stress

$$\frac{E}{1 - \nu^2} = \frac{200\times10^3}{1 - 0.09} = 219\,780\ \text{MPa}$$ $$[D] = 219\,780\begin{bmatrix} 1 & 0.3 & 0 \\ 0.3 & 1 & 0 \\ 0 & 0 & 0.35 \end{bmatrix} = \begin{bmatrix} 219\,780 & 65\,934 & 0 \\ 65\,934 & 219\,780 & 0 \\ 0 & 0 & 76\,923 \end{bmatrix}\ \text{MPa}$$

Note that $\dfrac{1 - \nu}{2} = 0.35$ and $219\,780(0.35) = 76\,923 = G = \dfrac{E}{2(1 + \nu)}$.

(b) Plane strain

$$\frac{E}{(1 + \nu)(1 - 2\nu)} = \frac{200\times10^3}{1.3(0.4)} = 384\,615\ \text{MPa}$$ $$[D] = 384\,615\begin{bmatrix} 0.7 & 0.3 & 0 \\ 0.3 & 0.7 & 0 \\ 0 & 0 & 0.2 \end{bmatrix} = \begin{bmatrix} 269\,231 & 115\,385 & 0 \\ 115\,385 & 269\,231 & 0 \\ 0 & 0 & 76\,923 \end{bmatrix}\ \text{MPa}$$

Stresses for $\varepsilon_x = 0.001$

The shear modulus term is the same in both cases. Plane strain gives higher stresses because $\varepsilon_z = 0$ is enforced.

Answer: plane stress $\sigma_x = 219.8$ MPa, $\sigma_y = 65.9$ MPa, $\sigma_z = 0$; plane strain $\sigma_x = 269.2$ MPa, $\sigma_y = 115.4$ MPa, $\sigma_z = 115.4$ MPa.

(a) Area and $[B]$

$$2A = x_1(y_2 - y_3) + x_2(y_3 - y_1) + x_3(y_1 - y_2) = 0 + 40(30) + 0 = 1200 \Rightarrow A = 600\ \text{mm}^2$$ $$[B] = \frac{1}{1200}\begin{bmatrix} -30 & 0 & 30 & 0 & 0 & 0 \\ 0 & -40 & 0 & 0 & 0 & 40 \\ -40 & -30 & 0 & 30 & 40 & 0 \end{bmatrix}$$

(b) Strains

$$\begin{aligned} \varepsilon_x &= \frac{30(0.02)}{1200} = 5.0\times10^{-4} \\ \varepsilon_y &= \frac{40(-0.01) - 40(0)}{1200} = -3.33\times10^{-4} \\ \gamma_{xy} &= \frac{30(0.005) + 40(0.004)}{1200} = \frac{0.31}{1200} = 2.58\times10^{-4} \end{aligned}$$

(c) Stresses (plane stress)

$\dfrac{E}{1 - \nu^2} = \dfrac{210\,000}{0.9375} = 224\,000$ MPa and $G = \dfrac{E}{2(1 + \nu)} = 84\,000$ MPa.

$$\begin{aligned} \sigma_x &= 224\,000\,(5.0\times10^{-4} + 0.25(-3.33\times10^{-4})) = 93.3\ \text{MPa} \\ \sigma_y &= 224\,000\,(-3.33\times10^{-4} + 0.25(5.0\times10^{-4})) = -46.7\ \text{MPa} \\ \tau_{xy} &= 84\,000\,(2.58\times10^{-4}) = 21.7\ \text{MPa} \end{aligned}$$

(d) Principal stresses

$$\sigma_{1,2} = \frac{\sigma_x + \sigma_y}{2} \pm \sqrt{\left(\frac{\sigma_x - \sigma_y}{2}\right)^2 + \tau_{xy}^2} = 23.3 \pm \sqrt{70.0^2 + 21.7^2} = 23.3 \pm 73.3$$ $$\sigma_1 = 96.6\ \text{MPa},\quad \sigma_2 = -50.0\ \text{MPa}$$ $$\tan 2\theta_p = \frac{2\tau_{xy}}{\sigma_x - \sigma_y} = \frac{43.4}{140} = 0.310 \Rightarrow \theta_p = 8.6^\circ$$

Answer: $\varepsilon = (5.0\times10^{-4},\ -3.33\times10^{-4},\ 2.58\times10^{-4})$; $\sigma_x = 93.3$ MPa, $\sigma_y = -46.7$ MPa, $\tau_{xy} = 21.7$ MPa; $\sigma_1 = 96.6$ MPa, $\sigma_2 = -50.0$ MPa at $8.6^\circ$ from the x-axis.

Displacement and strain

The element has nodes 1, 2, 3 with $\{d\} = \{u_1, v_1, u_2, v_2, u_3, v_3\}^T$. The displacement field is $\{u\} = [N]\{d\}$ with linear shape functions $N_i = \dfrac{1}{2A}(a_i + b_ix + c_iy)$.

Strains:

$$\{\varepsilon\} = \begin{Bmatrix} \partial u/\partial x \\ \partial v/\partial y \\ \partial u/\partial y + \partial v/\partial x \end{Bmatrix} = [B]\{d\},\quad [B] = \frac{1}{2A}\begin{bmatrix} b_1 & 0 & b_2 & 0 & b_3 & 0 \\ 0 & c_1 & 0 & c_2 & 0 & c_3 \\ c_1 & b_1 & c_2 & b_2 & c_3 & b_3 \end{bmatrix}$$

Stresses: $\{\sigma\} = [D]\{\varepsilon\} = [D][B]\{d\}$.

Potential energy

$$\Pi = \underbrace{\frac{1}{2}\int_V \{\varepsilon\}^T[D]\{\varepsilon\}\,dV}_{U} - \underbrace{\left(\int_V \{u\}^T\{X\}\,dV + \int_S \{u\}^T\{T\}\,dS + \{d\}^T\{P\}\right)}_{W}$$

Substituting $\{\varepsilon\} = [B]\{d\}$:

$$\Pi = \frac{1}{2}\{d\}^T\left[\int_V [B]^T[D][B]\,dV\right]\{d\} - \{d\}^T\{f\}$$

Minimising, $\dfrac{\partial\Pi}{\partial\{d\}} = 0$, gives $[k]\{d\} = \{f\}$ with

$$[k] = \int_V [B]^T[D][B]\,dV$$

Evaluation

For CST, $[B]$ and $[D]$ are constant, and the integral over the volume is just the element volume $tA$ ($t$ is the thickness):

$$\boxed{[k] = t\,A\,[B]^T[D][B]}$$

This is a symmetric $6\times6$ matrix, singular until supports are applied (three rigid-body modes: two translations and one rotation).

Load vectors

• Body force $\{X\} = \{X_x, X_y\}^T$ (per unit volume): $\{f_b\} = \displaystyle\int_V [N]^T\{X\}\,dV = \dfrac{tA}{3}\{X_x,\ X_y,\ X_x,\ X_y,\ X_x,\ X_y\}^T$, i.e. one-third of the total body force at each node.
• Uniform traction $\{T\} = \{T_x, T_y\}^T$ on an edge of length $\ell$ between nodes $i$ and $j$: $\{f_s\} = \displaystyle\int_S [N]^T\{T\}\,t\,dS$ gives $\dfrac{t\,\ell}{2}T_x$ and $\dfrac{t\,\ell}{2}T_y$ at each of the two edge nodes and zero at the third node.

Axisymmetric problems

A body of revolution (geometry and material symmetric about an axis) with loads and supports also symmetric about the same axis can be analysed on a single plane section, as a 2D problem in $r$-$z$. All quantities are independent of the circumferential angle $\theta$. Examples: pressure vessels, pipes, flywheels, rotors, pistons and shrink fits.

        z (axis)
        |    +------+
        |    | mesh |  <- triangular ring
        |    +------+     elements
        +------------- r

Each element is really a ring of material of volume $2\pi\bar r A$.

Unknowns and strains

Displacements: radial $u$ and axial $w$. Because $u$ displaces a ring, the circumference changes, producing a hoop strain.

$$\{\varepsilon\} = \begin{Bmatrix} \varepsilon_r \\ \varepsilon_\theta \\ \varepsilon_z \\ \gamma_{rz} \end{Bmatrix} = \begin{Bmatrix} \partial u/\partial r \\ u/r \\ \partial w/\partial z \\ \partial u/\partial z + \partial w/\partial r \end{Bmatrix}$$

For the triangular element with $u = \sum N_iu_i$:

$$[B] = \frac{1}{2A}\begin{bmatrix} \beta_1 & 0 & \beta_2 & 0 & \beta_3 & 0 \\ \frac{\alpha_1}{r} + \beta_1 + \frac{\gamma_1z}{r} & 0 & \frac{\alpha_2}{r} + \beta_2 + \frac{\gamma_2z}{r} & 0 & \frac{\alpha_3}{r} + \beta_3 + \frac{\gamma_3z}{r} & 0 \\ 0 & \gamma_1 & 0 & \gamma_2 & 0 & \gamma_3 \\ \gamma_1 & \beta_1 & \gamma_2 & \beta_2 & \gamma_3 & \beta_3 \end{bmatrix}$$

with $\alpha_i$, $\beta_i = z_j - z_k$ and $\gamma_i = r_k - r_j$ defined as for the CST using $(r, z)$.

Stress-strain matrix

$$[D] = \frac{E}{(1 + \nu)(1 - 2\nu)}\begin{bmatrix} 1 - \nu & \nu & \nu & 0 \\ \nu & 1 - \nu & \nu & 0 \\ \nu & \nu & 1 - \nu & 0 \\ 0 & 0 & 0 & \frac{1 - 2\nu}{2} \end{bmatrix}$$

Stiffness matrix

$$[k] = 2\pi\int_A [B]^T[D][B]\,r\,dr\,dz$$

The hoop-strain row of $[B]$ contains $1/r$ and $z/r$ terms, so $[B]$ varies with position, unlike the CST in the plane problem. The integral also contains the factor $r$. Therefore $[k]$ is not constant. In practice it is approximated by evaluating $[B]$ at the centroid $(\bar r, \bar z)$:

$$[k] \approx 2\pi\,\bar r\,A\,[\bar B]^T[D][\bar B]$$

This is accurate when elements are small compared with their distance from the axis. Loads are applied as total forces on the full ring, $2\pi r$ times the load per unit circumference.

Area and centroid

$$2A = r_1(z_2 - z_3) + r_2(z_3 - z_1) + r_3(z_1 - z_2) = 20(-30) + 40(30) + 0 = 600 \Rightarrow A = 300\ \text{mm}^2$$ $$\bar r = \frac{20 + 40 + 20}{3} = 26.67\ \text{mm},\qquad \bar z = \frac{0 + 0 + 30}{3} = 10\ \text{mm}$$

Coefficients

$\alpha_1 = r_2z_3 - r_3z_2 = 1200$, $\alpha_2 = r_3z_1 - r_1z_3 = -600$, $\alpha_3 = r_1z_2 - r_2z_1 = 0$.

Hoop-strain coefficient at the centroid, $f_i = \dfrac{\alpha_i}{\bar r} + \beta_i + \dfrac{\gamma_i\bar z}{\bar r}$:

$$f_1 = \frac{1200}{26.67} - 30 - \frac{200}{26.67} = 45 - 30 - 7.5 = 7.5,\quad f_2 = -22.5 + 30 + 0 = 7.5,\quad f_3 = 0 + 0 + 7.5 = 7.5$$

(Each is $7.5$ because $N_i = 1/3$ at the centroid, so $f_i = 2A/(3\bar r) = 7.5$.)

$[B]$ at the centroid

$$[\bar B] = \frac{1}{600}\begin{bmatrix} -30 & 0 & 30 & 0 & 0 & 0 \\ 7.5 & 0 & 7.5 & 0 & 7.5 & 0 \\ 0 & -20 & 0 & 0 & 0 & 20 \\ -20 & -30 & 0 & 30 & 20 & 0 \end{bmatrix}$$

Strains

$$\begin{aligned} \varepsilon_r &= \frac{-30(0.01) + 30(0.02)}{600} = 5.0\times10^{-4} \\ \varepsilon_\theta &= \frac{7.5(0.01 + 0.02 + 0)}{600} = 3.75\times10^{-4} \\ \varepsilon_z &= \frac{-20(0) + 0(0.01) + 20(0.015)}{600} = 5.0\times10^{-4} \\ \gamma_{rz} &= \frac{-20(0.01) + (-30)(0) + 0 + 30(0.01) + 0 + 0}{600} = 1.67\times10^{-4} \end{aligned}$$

Stresses

$\dfrac{E}{(1 + \nu)(1 - 2\nu)} = 384\,615$ MPa.

$$\begin{aligned} \sigma_r &= 384\,615\,[0.7(5.0) + 0.3(3.75 + 5.0)]\times10^{-4} = 235.6\ \text{MPa} \\ \sigma_\theta &= 384\,615\,[0.7(3.75) + 0.3(5.0 + 5.0)]\times10^{-4} = 216.3\ \text{MPa} \\ \sigma_z &= 384\,615\,[0.7(5.0) + 0.3(5.0 + 3.75)]\times10^{-4} = 235.6\ \text{MPa} \\ \tau_{rz} &= 384\,615(0.2)(1.67\times10^{-4}) = 12.8\ \text{MPa} \end{aligned}$$

Answer: $A = 300\ \text{mm}^2$, centroid $(26.67,\ 10)$ mm; $\varepsilon_r = 5.0\times10^{-4}$, $\varepsilon_\theta = 3.75\times10^{-4}$, $\varepsilon_z = 5.0\times10^{-4}$, $\gamma_{rz} = 1.67\times10^{-4}$; $\sigma_r = 235.6$, $\sigma_\theta = 216.3$, $\sigma_z = 235.6$, $\tau_{rz} = 12.8$ MPa.

Thermal load in FEM

If an element is free to expand, there is no stress. If the expansion is prevented, the initial strain $\varepsilon_0 = \alpha\Delta T$ produces an equivalent nodal force vector

$$\{f_T\} = \int_V [B]^T[D]\{\varepsilon_0\}\,dV$$

For a bar element this gives $\{f_T\} = EA\alpha\Delta T\begin{Bmatrix} -1 \\ +1 \end{Bmatrix}$. Stress is found from the mechanical strain only: $\sigma = E(\varepsilon - \alpha\Delta T)$.

Element data

 wall|==[A1=500]==o==[A2=250]==|wall
      L1=300      2    L2=200
      1                       3

$$k_1 = \frac{500(200\times10^3)}{300} = 3.333\times10^5\ \text{N/mm},\qquad k_2 = \frac{250(200\times10^3)}{200} = 2.5\times10^5\ \text{N/mm}$$ $$E A_1\alpha\Delta T = 200\times10^3(500)(12\times10^{-6})(60) = 72\,000\ \text{N},\qquad EA_2\alpha\Delta T = 36\,000\ \text{N}$$

Equations

$u_1 = u_3 = 0$; only $u_2$ is unknown. The thermal force at node 2 is $+72\,000$ from element 1 and $-36\,000$ from element 2:

$$(k_1 + k_2)\,u_2 = 72\,000 - 36\,000 = 36\,000$$ $$u_2 = \frac{36\,000}{5.833\times10^5} = 0.0617\ \text{mm (to the right)}$$

Element forces and stresses

$$\begin{aligned} N_1 &= k_1(u_2 - 0) - EA_1\alpha\Delta T = 20\,571 - 72\,000 = -51\,429\ \text{N} \\ N_2 &= k_2(0 - u_2) - EA_2\alpha\Delta T = -15\,429 - 36\,000 = -51\,429\ \text{N} \end{aligned}$$

The two forces are equal, as equilibrium of the junction requires.

$$\sigma_1 = \frac{-51\,429}{500} = -102.9\ \text{MPa},\qquad \sigma_2 = \frac{-51\,429}{250} = -205.7\ \text{MPa}$$

Check

Free expansion would be $\alpha\Delta T(L_1 + L_2) = 0.36$ mm. It is cancelled by the compression of the two bars: $P\left(\dfrac{300}{500E} + \dfrac{200}{250E}\right) = \dfrac{P(1.4)}{200\times10^3} = 0.36$ gives $P = 51\,429$ N, which agrees.

Answer: $u_2 = 0.0617$ mm; $\sigma_1 = -102.9$ MPa and $\sigma_2 = -205.7$ MPa (compressive); wall reaction $= 51.43$ kN.