8.7  Stress and strain states

The one-dimensional case as a starting point

For an infinitesimal element of length \(\Delta x\), the displacement at \(x\) is \(u(x)\) and at the other end \(u(x+\Delta x)\), so the strain is the difference divided by the length,

\[ \varepsilon := \frac{\text{elongation}}{\text{original length}} = \frac{u(x+\Delta x) - u(x)}{\Delta x} \longrightarrow \frac{du}{dx} \]

and the stress follows from Hooke’s law \(\sigma = E\varepsilon\).

Shear is the other elementary deformation. A shear stress \(\tau\) tilts the element through the shear angle \(\gamma\), and the constitutive relation is \(\tau = G\gamma\).

The general state of strain

In two dimensions the displacement is a vector field \(\bm u = (u_x, u_y)\) and the element can stretch in two directions while simultaneously shearing. Figure 8.7.3 shows the four independent deformations: two extensions and two contributions to the change of angle.

The extension in the \(x\) direction comes from the variation of \(u_x\) along \(x\), and the extension in the \(y\) direction from the variation of \(u_y\) along \(y\),

\[ \frac{u_x(x+\Delta x, y) - u_x(x,y)}{\Delta x} \longrightarrow \frac{\partial u_x}{\partial x} =: \varepsilon_x, \qquad \frac{u_y(x, y+\Delta y) - u_y(x,y)}{\Delta y} \longrightarrow \frac{\partial u_y}{\partial y} =: \varepsilon_y \]

The shear contribution is the tilt of an originally vertical or horizontal edge. Moving along \(x\), the vertical displacement changes by \(u_y(x+\Delta x, y) - u_y(x,y)\), and trigonometry gives the tilt angle \(\gamma_1\) of the originally horizontal edge as

\[ \tan\gamma_1 = \frac{u_y(x+\Delta x, y) - u_y(x,y)}{\Delta x} \longrightarrow \frac{\partial u_y}{\partial x} \]

Assuming deformations small enough that \(\tan\gamma_1 \approx \gamma_1\), and adding the tilt \(\gamma_2 = \partial u_x/\partial y\) of the originally vertical edge, the total change of the right angle between the two edges is

\[ \gamma_{xy} = \gamma_1 + \gamma_2 = \frac{\partial u_y}{\partial x} + \frac{\partial u_x}{\partial y} \tag{8.7.1}\]

The same reasoning in three dimensions adds one more extension and two more shear angles,

\[ \varepsilon_z = \frac{\partial u_z}{\partial z}, \qquad \gamma_{xz} = \frac{\partial u_x}{\partial z} + \frac{\partial u_z}{\partial x}, \qquad \gamma_{yz} = \frac{\partial u_y}{\partial z} + \frac{\partial u_z}{\partial y} \]

giving six independent strains in total. Collecting them into a matrix, with the tensor shear strains \(\tfrac{1}{2}\gamma\) off the diagonal, gives the strain tensor

\[ \bm\varepsilon = \begin{bmatrix} \varepsilon_x & \tfrac{1}{2}\gamma_{xy} & \tfrac{1}{2}\gamma_{xz}\\[4pt] \tfrac{1}{2}\gamma_{xy} & \varepsilon_y & \tfrac{1}{2}\gamma_{yz}\\[4pt] \tfrac{1}{2}\gamma_{xz} & \tfrac{1}{2}\gamma_{yz} & \varepsilon_z \end{bmatrix} \tag{8.7.2}\]

All six entries follow from the single index-notation definition

\[ \boxed{\varepsilon_{ij} = \frac{1}{2}\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right)} \tag{8.7.3}\]

or in direct notation \(\bm\varepsilon = \tfrac{1}{2}\left(\nabla\bm u + \nabla\bm u^\mathsf{T}\right)\), the symmetric part of the displacement gradient. Written out, that gradient is the matrix of every derivative of every displacement component,

\[ \nabla\bm u = \begin{bmatrix} \dfrac{\partial u_x}{\partial x} & \dfrac{\partial u_x}{\partial y} & \dfrac{\partial u_x}{\partial z}\\[8pt] \dfrac{\partial u_y}{\partial x} & \dfrac{\partial u_y}{\partial y} & \dfrac{\partial u_y}{\partial z}\\[8pt] \dfrac{\partial u_z}{\partial x} & \dfrac{\partial u_z}{\partial y} & \dfrac{\partial u_z}{\partial z} \end{bmatrix}, \qquad (\nabla\bm u)_{ij} = \frac{\partial u_i}{\partial x_j} \tag{8.7.4}\]

with the component differentiated running down the rows and the direction of differentiation across the columns. Nine derivatives go in, and adding the transpose averages each entry with its mirror across the diagonal, which leaves the six independent strains of 8.7.2. Checking a diagonal entry, \(\varepsilon_{11} = \tfrac{1}{2}(\partial u_1/\partial x_1 + \partial u_1/\partial x_1) = \partial u_x/\partial x\), recovers \(\varepsilon_x\). Checking an off-diagonal entry, \(\varepsilon_{12} = \tfrac{1}{2}(\partial u_1/\partial x_2 + \partial u_2/\partial x_1) = \tfrac{1}{2}\gamma_{xy}\), gives half the engineering shear strain.

gamma = sp.symbols('gamma', positive=True)

grad_u = sp.Matrix([[0, gamma],                  # simple shear
                    [0, 0]])
mm = sp.Matrix([1, 1])/sp.sqrt(2)                # unit vector along the diagonal

eps       = (grad_u + grad_u.T)/2                # symmetric part, with the factor 1/2
eps_naive = sp.Matrix([[0, gamma], [gamma, 0]])  # full gamma off the diagonal

stretched = (sp.eye(2) + grad_u)*mm
geometric = sp.series(stretched.norm() - 1, gamma, 0, 2).removeO()

# A cell renders only its last expression, so all three rows go in one call.
ltx(r"\text{geometry:}\quad |(\bm I + \nabla\bm u)\bm m| - 1 &=", geometric, r"+ O(\gamma^2)",
    r"\\ \bm\varepsilon =", eps,
    r"\;\Rightarrow\;\; \bm m^\mathsf{T}\bm\varepsilon\,\bm m &=", sp.simplify((mm.T*eps*mm)[0]),
    r"\\ \text{naive} =", eps_naive,
    r"\;\Rightarrow\;\; \bm m^\mathsf{T}(\cdot)\,\bm m &=", sp.simplify((mm.T*eps_naive*mm)[0]),
    aligned=True)

The stress tensor

The stresses are collected the same way, with normal stresses on the diagonal and shear stresses off it,

\[ \bm\sigma = \begin{bmatrix} \sigma_x & \tau_{xy} & \tau_{xz}\\ \tau_{xy} & \sigma_y & \tau_{yz}\\ \tau_{xz} & \tau_{yz} & \sigma_z \end{bmatrix} \tag{8.7.6}\]

where \(\sigma_x\) is the normal stress on the face whose normal points along \(x\), and \(\tau_{xy}\) is the stress on that same face acting in the \(y\) direction. Reading the matrix by rows therefore gives, for each face of the element, the full traction vector acting on it.

Equilibrium of an infinitesimal element

The stresses cannot be arbitrary functions of position. Equilibrium of a small block of material relates their derivatives to the body force. Consider the element in Figure 8.7.5, of size \(\Delta x\) by \(\Delta y\), loaded by a body force per unit volume \(\bm f\).

Summing forces in the \(x\) direction, each stress is multiplied by the area of the face it acts on, and the body force by the volume of the element,

\[ \sum F_x = 0: \quad \Delta\sigma_x\, \Delta y + \Delta\tau_{yx}\, \Delta x + f_x\, \Delta x \Delta y = 0 \]

which written out in full reads

\[ \left(\sigma_x(x+\Delta x, y) - \sigma_x(x,y)\right)\Delta y + \left(\tau_{yx}(x, y+\Delta y) - \tau_{yx}(x,y)\right)\Delta x + f_x\,\Delta x \Delta y = 0 \]

Dividing by \(\Delta x \Delta y\) and letting the element shrink turns the differences into partial derivatives. Doing the same in the \(y\) direction gives the pair

\[ \begin{array}{l} -\dfrac{\partial \sigma_x}{\partial x} - \dfrac{\partial \tau_{yx}}{\partial y} = f_x\\[8pt] -\dfrac{\partial \tau_{xy}}{\partial x} - \dfrac{\partial \sigma_y}{\partial y} = f_y \end{array} \tag{8.7.7}\]

Each line is the divergence of one row of the stress tensor. Recalling the gradient operator

\[ \nabla = \begin{bmatrix} \dfrac{\partial}{\partial x} \\[8pt] \dfrac{\partial}{\partial y}\end{bmatrix} \]

the two equations collapse into the single statement

\[ \boxed{-\nabla \cdot \bm\sigma = \bm f} \tag{8.7.8}\]

which is the equilibrium equation of continuum mechanics. It is the direct descendant of \(\sum F_i = 0\) from statics, written for a body in which the internal force varies continuously from point to point. Every structural problem in this book is 8.7.8 together with a constitutive law and a set of boundary conditions, and the various beam and rod and truss theories are what 8.7.8 becomes once a particular geometry is assumed.

Boundary traction

Equilibrium governs the interior. At the surface, the stress in the material must balance whatever is applied from outside. Let \(\Omega\) denote the body and \(\partial\Omega\) its boundary, and let an external force \(\bm t\) act on part of that boundary. This force is a traction, measured in force per unit length in two dimensions and force per unit area in three.

Cut a small triangular wedge from the boundary, as in Figure 8.7.7, with the traction acting on the hypotenuse of length \(\Delta s\) and the two straight legs lying along the coordinate directions. The leg with outward normal \(-\bm e_x\) has length \(\Delta y\) and the leg with outward normal \(-\bm e_y\) has length \(\Delta x\), and the outward unit normal of the hypotenuse has components

\[ n_x = \frac{\Delta y}{\Delta s}, \qquad n_y = \frac{\Delta x}{\Delta s} \]

Force equilibrium of the wedge in the two coordinate directions gives

\[ \begin{array}{l} t_x\, \Delta s = \sigma_x\, \Delta y + \tau_{xy}\, \Delta x\\[4pt] t_y\, \Delta s = \tau_{xy}\, \Delta y + \sigma_y\, \Delta x \end{array} \]

and dividing through by \(\Delta s\) replaces the lengths by the components of the normal,

\[ \left.\begin{array}{l} \sigma_x n_x + \tau_{xy} n_y = t_x\\ \tau_{xy} n_x + \sigma_y n_y = t_y \end{array}\right\} \quad\Longleftrightarrow\quad \begin{bmatrix} \sigma_x & \tau_{xy}\\ \tau_{xy} & \sigma_y\end{bmatrix} \begin{bmatrix} n_x \\ n_y \end{bmatrix} = \begin{bmatrix} t_x \\ t_y \end{bmatrix} \]

which is Cauchy’s formula

\[ \boxed{\bm\sigma \cdot \bm n = \bm t} \tag{8.7.9}\]

The tensor is now doing exactly the job we asked of it at the start of the chapter. Feed it any direction \(\bm n\) and it returns the traction acting across a cut with that normal. The stress at a point is not a number and not even a list of three numbers, it is the linear map that turns cut directions into tractions.

A part of the boundary where the traction is prescribed, as a distributed load, a pressure \(\bm t = -p\,\bm n\) or a force, carries a Neumann boundary condition: through 8.7.9 it fixes the stress there, as \(EA\,du/dx(L) = F\) fixed the normal force at the end of the rod’s boundary value problem.

Example 1: Recovering the inclined cut

In Chapter 8.4 we computed the normal and shear stress on a cut inclined at an angle \(\varphi\) to the axis of a bar in uniaxial tension, and obtained 8.4.2 by writing equilibrium for that particular cut. Cauchy’s formula should reproduce the result without any new free body diagram, and it does.

A bar pulled along \(x\) with \(\sigma_0 = P/A\) has the stress tensor \(\bm\sigma = \operatorname{diag}(\sigma_0, 0)\): nothing acts in \(y\), and there is no shear in this frame. A cut at angle \(\varphi\) to the bar axis has unit normal \(\bm n = (\sin\varphi, -\cos\varphi)\) and unit tangent \(\bm m = (\cos\varphi, \sin\varphi)\). The traction on the cut is \(\bm t = \bm\sigma \bm n\), and its normal and tangential components are the projections of \(\bm t\) onto \(\bm n\) and \(\bm m\).

sigma_0, varphi = sp.symbols('sigma_0 varphi', real=True)

sigma_uniaxial = sp.Matrix([[sigma_0, 0],
                            [0,       0]])

nn = sp.Matrix([sp.sin(varphi), -sp.cos(varphi)])   # normal of the cut
mm = sp.Matrix([sp.cos(varphi),  sp.sin(varphi)])   # tangent of the cut

tt = sigma_uniaxial * nn
sigma_phi = sp.simplify((tt.T * nn)[0])
tau_phi   = sp.simplify((tt.T * mm)[0])

ltx(r"\bm t &=", tt,
    r"\\ \sigma_\varphi &=", sigma_phi,
    r"\\ \tau_\varphi &=", sp.simplify(sp.expand_trig(tau_phi)), aligned=True)

The results are \(\sigma_\varphi = \sigma_0\sin^2\varphi\) and \(\tau_\varphi = \sigma_0\sin\varphi\cos\varphi\), identical to 8.4.2 with \(\sigma_0 = P/A\). The earlier derivation needed a free body diagram for each new cut angle. The tensor needs none: one object encodes the response to every possible cut, and Cauchy’s formula extracts the answer for whichever cut we ask about.

Principal stress and strain

Consider a plate in pure tension along \(x\), as in Figure 8.7.8. Aligned with the load, the stress tensor is diagonal, with a normal stress along the load, some transverse component, and no shear at all.

Rotate the frame of reference without touching the plate or the load, and shear stresses appear, exactly as they did in Example 1. Nothing physical changed, only our description of it. Reversing the question gives something useful: given a stress tensor with shear in it, is there a rotation that removes the shear entirely? The directions in which it vanishes are the principal directions, and the normal stresses in those directions are the principal stresses.

The shear on a cut vanishes when the traction is parallel to the normal, that is when

\[ \bm t = \lambda \bm n \]

for some scalar \(\lambda\). Substituting Cauchy’s formula 8.7.9 gives

\[ \bm\sigma \cdot \bm n = \lambda \bm n \quad\Longleftrightarrow\quad (\bm\sigma - \lambda \bm I)\cdot \bm n = \bm 0 \tag{8.7.10}\]

which is an eigenvalue problem. The principal stresses are the eigenvalues of the stress tensor and the principal directions are its eigenvectors. Because the tensor is symmetric, the eigenvalues are real and the eigenvectors are mutually orthogonal, so the principal directions always form a proper coordinate frame.

Non-trivial solutions require the matrix to be singular, \(\det(\bm\sigma - \lambda\bm I) = 0\), which in two dimensions expands to

\[ \begin{vmatrix} \sigma_x - \lambda & \tau_{xy} \\ \tau_{xy} & \sigma_y - \lambda \end{vmatrix} = (\sigma_x - \lambda)(\sigma_y - \lambda) - \tau_{xy}^2 = 0 \]

or, collecting terms,

\[ \lambda^2 - \underbrace{(\sigma_x + \sigma_y)}_{\operatorname{tr}\bm\sigma}\lambda + \underbrace{\sigma_x\sigma_y - \tau_{xy}^2}_{\det\bm\sigma} = 0 \tag{8.7.11}\]

whose roots are the two principal stresses. The coefficients are the trace and the determinant, both of which are unchanged by rotation, which is the algebraic statement of the physical fact that the principal stresses do not depend on the frame we happened to start in. Solving the quadratic gives the closed form

\[ \sigma_{1,2} = \frac{\sigma_x + \sigma_y}{2} \pm \sqrt{\left(\frac{\sigma_x - \sigma_y}{2}\right)^2 + \tau_{xy}^2} \tag{8.7.12}\]

and the identical construction applied to the strain tensor gives the principal strains, with \(\tfrac{1}{2}\gamma_{xy}\) in place of \(\tau_{xy}\),

\[ \varepsilon_{1,2} = \frac{\varepsilon_x + \varepsilon_y}{2} \pm \sqrt{\left(\frac{\varepsilon_x - \varepsilon_y}{2}\right)^2 + \left(\frac{\gamma_{xy}}{2}\right)^2} \tag{8.7.13}\]

The shear stress is largest on the planes at \(45^\circ\) to the principal directions, where it reaches the maximum shear stress

\[ \tau_{\max} = \frac{\sigma_{\max} - \sigma_{\min}}{2} \tag{8.7.14}\]

half the difference between the largest and the smallest of the three principal stresses. Under uniaxial tension \(\sigma_0\) the principal stresses are \(\sigma_0\), \(0\) and \(0\), so \(\tau_{\max} = \sigma_0/2\) acts on the planes at \(45^\circ\) to the load, as 8.4.2 found by cutting. Mohr’s circle draws the two in-plane principal stresses as the points where a circle crosses the \(\sigma\) axis and the largest in-plane shear stress as its radius, and shows how they move as the frame turns.

Let us confirm 8.7.12 by solving the eigenvalue problem symbolically.

sigma_x, sigma_y, tau_xy = sp.symbols('sigma_x sigma_y tau_xy', real=True)

S = sp.Matrix([[sigma_x, tau_xy],
               [tau_xy,  sigma_y]])

eigs = sp.simplify(sp.Matrix(sorted(S.eigenvals().keys(), key=sp.default_sort_key)))
ltx(r"\lambda_{1,2} =", eigs)

sigma_z, tau_xz, tau_yz, lam = sp.symbols('sigma_z tau_xz tau_yz lambda', real=True)

Sigma3 = sp.Matrix([[sigma_x, tau_xy, tau_xz],
                    [tau_xy,  sigma_y, tau_yz],
                    [tau_xz,  tau_yz,  sigma_z]])

I_1 = sigma_x + sigma_y + sigma_z                                          # @eq-invariant-1
I_2 = sigma_x*sigma_y + sigma_y*sigma_z + sigma_z*sigma_x - tau_xy**2 - tau_yz**2 - tau_xz**2
I_3 = (sigma_x*sigma_y*sigma_z + 2*tau_xy*tau_yz*tau_xz
       - sigma_x*tau_yz**2 - sigma_y*tau_xz**2 - sigma_z*tau_xy**2)
cubic = lam**3 - I_1*lam**2 + I_2*lam - I_3                                # @eq-characteristic-3d

# det(sigma - lambda I) is minus the cubic, so their sum must vanish identically
plane_stress = {sigma_z: 0, tau_xz: 0, tau_yz: 0}
ltx(r"\det(\bm\sigma - \lambda\bm I) + \lambda^3 - I_1\lambda^2 + I_2\lambda - I_3 &=",
    sp.expand((Sigma3 - lam*sp.eye(3)).det() + cubic),
    r"\\ \text{plane stress:}\quad \lambda^3 - I_1\lambda^2 + I_2\lambda - I_3 &=",
    sp.factor(cubic.subs(plane_stress)),
    aligned=True)

Effective stress and the von Mises criterion

A design check needs a single number to compare against the yield strength, and a stress state has six components. Reducing them to one number requires a hypothesis about what causes yielding. The hypothesis that works best for ductile metals is that yielding is driven by distortion rather than by volume change, which is the von Mises criterion. In terms of the three principal stresses it reads

\[ \sigma_{\text{vM}} = \sqrt{\tfrac{1}{2}\left[(\sigma_1 - \sigma_2)^2 + (\sigma_2 - \sigma_3)^2 + (\sigma_3 - \sigma_1)^2\right]} \tag{8.7.21}\]

or, in terms of the components in any frame,

\[ \sigma_{\text{vM}} = \sqrt{\tfrac{1}{2}\left[(\sigma_x - \sigma_y)^2 + (\sigma_y - \sigma_z)^2 + (\sigma_z - \sigma_x)^2 + 6(\tau_{xy}^2 + \tau_{yz}^2 + \tau_{zx}^2)\right]} \tag{8.7.22}\]

The two expressions give the same number, since 8.7.21 is written in the principal frame where the shear terms vanish. Note that adding the same pressure to all three normal stresses leaves \(\sigma_{\text{vM}}\) unchanged, because only differences of principal stresses appear. Hydrostatic pressure does not cause yielding, which is why material deep in the ocean does not yield under its own ambient pressure.

For plane stress, where \(\sigma_z = \tau_{yz} = \tau_{zx} = 0\), 8.7.22 simplifies to

\[ \boxed{\sigma_{\text{vM}} = \sqrt{\sigma_x^2 - \sigma_x\sigma_y + \sigma_y^2 + 3\tau_{xy}^2}} \tag{8.7.23}\]

which is the form that appears in most hand calculations, and the quantity that a finite element post-processor plots by default. Chapter 8.10 derives the criterion from the distortion energy and compares it with the Tresca and maximum normal stress criteria.

Example 2: Evaluating a plane stress state

A point in a loaded plate is in the plane stress state \(\sigma_x = 120~\text{MPa}\), \(\sigma_y = -40~\text{MPa}\), \(\tau_{xy} = 50~\text{MPa}\), Figure 8.7.11. Determine the principal stresses, the principal directions, the maximum shear stress and the von Mises stress. The material yields at \(250~\text{MPa}\): is the point safe?

The normal \(\bm e_z\) of the plate carries no stress, so it is a principal direction with principal stress zero. The in-plane principal stresses follow from 8.7.12,

\[ \sigma_{\text{I},\text{II}} = \frac{120 - 40}{2} \pm \sqrt{\left(\frac{120 + 40}{2}\right)^2 + 50^2}~\text{MPa} = 40 \pm 94.3~\text{MPa} \]

and the direction \(\bm n_\text{I} = (\cos\theta_\text{I}, \sin\theta_\text{I})\) from the first row of \((\bm\sigma - \sigma_\text{I}\bm I)\bm n_\text{I} = \bm 0\), 8.7.10,

\[ (120 - \sigma_\text{I})\cos\theta_\text{I} + 50\sin\theta_\text{I} = 0 \quad\Longrightarrow\quad \tan\theta_\text{I} = \frac{\sigma_\text{I} - 120}{50} \]

with \(\bm n_\text{II}\) perpendicular to it. SymPy solves the eigenvalue problem of the \(2\times 2\) tensor.

S_num = sp.Matrix([[120, 50],
                   [50, -40]])

# eigenvects gives (value, multiplicity, [vector]); the larger in-plane value is sigma_I
(sigma_I, _, (n_I,)), (sigma_II, _, (n_II,)) = sorted(S_num.eigenvects(), key=lambda e: -e[0])
n_I, n_II = n_I.normalized(), n_II.normalized()
n_I = n_I if n_I[0] > 0 else -n_I        # an eigenvector's sign is free; point it along +x
theta_I = sp.deg(sp.atan2(n_I[1], n_I[0]))

ltx(r"\sigma_\text{I} &=", float(sigma_I), r"~\text{MPa}, \qquad \bm n_\text{I} =", n_I.evalf(3),
    r"\\ \sigma_\text{II} &=", float(sigma_II), r"~\text{MPa}, \qquad \bm n_\text{II} =", n_II.evalf(3),
    r"\\ \theta_\text{I} &=", float(theta_I), r"^\circ",
    aligned=True)

Figure 8.7.12 draws the result onto the element of Figure 8.7.11. Turned by \(\theta_\text{I}\), the element carries only the two principal stresses: tension along \(\bm n_\text{I}\) and compression along \(\bm n_\text{II}\), and no shear on any of its faces.

from matplotlib.patches import FancyArrowPatch, Polygon
from mechanicskit import sketch as sk

def arrow(ax, p0, p1, color=sk.LOAD, lw=1.8):
    ax.add_patch(FancyArrowPatch(p0, p1, arrowstyle='-|>', mutation_scale=12, color=color,
                                 lw=lw, shrinkA=0, shrinkB=0, zorder=6))

def element(ax, theta, normal, shear, h=1.0, k=1/120, g=0.12):
    """A square of half side h turned by theta, with a normal stress on each pair of faces
    (normal[0] on the faces along e1, normal[1] on those along e2) and the shear on all four."""
    c, s = np.cos(theta), np.sin(theta)
    e1, e2 = np.array([c, s]), np.array([-s, c])
    ax.add_patch(Polygon([h*(-e1 - e2), h*(e1 - e2), h*(e1 + e2), h*(-e1 + e2)], closed=True,
                         facecolor=sk.STEEL, edgecolor=sk.EDGE, lw=sk.LW, zorder=2))
    for n, sig in ((e1, normal[0]), (-e1, normal[0]), (e2, normal[1]), (-e2, normal[1])):
        L = 0.25 + k*abs(sig)
        face = (h + g)*n
        # tension leaves the face, compression points into it
        arrow(ax, *((face, face + L*n) if sig > 0 else (face + L*n, face)))
        sk.label(ax, face + (L + 0.25)*n, f'${abs(sig):.0f}$', color=sk.LOAD, fontsize=12)
    if shear != 0:
        Ls = 0.25 + k*abs(shear)
        for n, t in ((e1, e2), (-e1, -e2), (e2, e1), (-e2, -e1)):
            mid = (h + 0.1)*n
            arrow(ax, mid - Ls/2*t, mid + Ls/2*t)
            sk.label(ax, mid + (Ls/2 + 0.05)*t + 0.2*n, f'${abs(shear):.0f}$',
                     color=sk.LOAD, fontsize=12)
    return e1, e2

th = float(sp.rad(theta_I))
fig, (ax, bx) = plt.subplots(1, 2, figsize=(10, 4.8))
for a in (ax, bx):
    a.set_aspect('equal')
    a.axis('off')
    a.set_xlim(-2.7, 2.7)
    a.set_ylim(-2.5, 2.5)

# the given state, as in the problem statement
element(ax, 0.0, (120, -40), 50)
# the principal state: sigma_I along n_I, sigma_II along n_II, no shear
e1, e2 = element(bx, th, (float(sigma_I), float(sigma_II)), 0)
for a in (ax, bx):                                   # the frame at the lower left
    o = np.array([-2.35, -2.1])
    for d, name in (((1, 0), '$x$'), ((0, 1), '$y$')):
        arrow(a, o, o + 0.6*np.array(d), color=sk.EDGE, lw=1.2)
        sk.label(a, o + 0.78*np.array(d), name, fontsize=12)
# the principal directions through the centre, and the angle from the x axis
for d, name, side in ((e1, r'$\mathbf{n}_\mathrm{I}$', e2), (e2, r'$\mathbf{n}_\mathrm{II}$', -e1)):
    bx.plot(*np.c_[-2.4*d, 2.4*d], color=sk.BLUE, lw=0.9, ls=(0, (6, 3)), zorder=1)
    sk.label(bx, 2.1*d + 0.38*side, name, color=sk.BLUE, fontsize=12)
# the angle from the x direction to n_I, drawn outside the element
r0, r1, ra = 1.3, 2.2, 1.95
bx.plot([r0, r1], [0, 0], color=sk.GREY, lw=0.9, zorder=3)
t = np.linspace(0, th, 30)
bx.plot(ra*np.cos(t), ra*np.sin(t), color=sk.GREY, lw=0.9, zorder=3)
sk.label(bx, (ra + 0.35)*np.array([np.cos(th/2), np.sin(th/2)]) + np.array([0.05, -0.02]),
         rf'${float(theta_I):.0f}^{{\circ}}$', fontsize=12)
sk.label(ax, (-2.55, 2.35), '(a)', fontsize=12, ha='left')
sk.label(bx, (-2.55, 2.35), '(b)', fontsize=12, ha='left')
fig.tight_layout()
plt.show()

The in-plane principal stresses are \(\sigma_\text{I} = 134~\text{MPa}\) in tension and \(\sigma_\text{II} = -54~\text{MPa}\) in compression. The direction \(\bm n_\text{I}\) makes an angle of \(16^\circ\) with the \(x\) axis, which is the amount the frame must be rotated for the shear to disappear.

The pair has opposite signs, the second case of the plane stress ordering, so the zero normal to the plate falls between them and the ordered principal stresses 8.7.20 are

\[ \sigma_1 = \sigma_\text{I} = 134~\text{MPa}, \qquad \sigma_2 = 0, \qquad \sigma_3 = \sigma_\text{II} = -54~\text{MPa} \]

The maximum shear stress in the plane is half the difference of the in-plane pair, a result that follows from 8.4.2 applied to the principal frame, and the largest shear stress over all planes is 8.7.14,

\[ \begin{aligned} \tau_{\max,\text{plane}} &= \frac{\sigma_\text{I} - \sigma_\text{II}}{2}\\[4pt] \tau_{\max} &= \frac{\sigma_1 - \sigma_3}{2} \end{aligned} \]

The von Mises stress follows from 8.7.23,

\[ \sigma_{\text{vM}} = \sqrt{120^2 - 120\cdot(-40) + (-40)^2 + 3\cdot 50^2}~\text{MPa} = \sqrt{28\,300}~\text{MPa} \]

and 8.7.21 with the ordered three must give the same number. The factor of safety against yielding is \(\eta = 250~\text{MPa}/\sigma_{\text{vM}}\).

# the zero normal to the plate joins the in-plane pair; sorting puts the three in order
sigma_1, sigma_2, sigma_3 = sorted([sigma_I, sigma_II, sp.Integer(0)], reverse=True)

tau_max_plane = (sigma_I - sigma_II)/2
tau_max = (sigma_1 - sigma_3)/2                                              # @eq-max-shear-stress

sx, sy, txy = 120, -40, 50
sigma_vm = sp.sqrt(sx**2 - sx*sy + sy**2 + 3*txy**2)                         # @eq-von-mises-plane
sigma_vm_123 = sp.sqrt(((sigma_1 - sigma_2)**2 + (sigma_2 - sigma_3)**2
                        + (sigma_3 - sigma_1)**2)/2)                         # @eq-von-mises-principal

ltx(r"\sigma_1, \sigma_2, \sigma_3 &=", float(sigma_1), r",\;", float(sigma_2), r",\;",
    float(sigma_3), r"~\text{MPa}",
    r"\\ \tau_{\max,\text{plane}} &=", float(tau_max_plane), r"~\text{MPa}, \qquad \tau_{\max} =",
    float(tau_max), r"~\text{MPa}",
    r"\\ \sigma_{\text{vM}} &=", float(sigma_vm), r"~\text{MPa}, \qquad \text{from } \sigma_{1,2,3}:\;",
    float(sigma_vm_123), r"~\text{MPa}",
    r"\\ \eta &=", float(250/sigma_vm), aligned=True)

The von Mises stress is \(168~\text{MPa}\) against a yield strength of \(250~\text{MPa}\), so the point is elastic with a factor of safety of \(1.49\). the chapter on Mohr’s circle draws this state as a circle, and its interactive figure loads it as a preset. The two shear values coincide here because the in-plane pair has opposite signs: the zero normal to the plate is \(\sigma_2\), and the pair itself already spans \(\sigma_1\) and \(\sigma_3\). Had both in-plane principal stresses been positive, the zero would have been \(\sigma_3\), the absolute maximum would have been \(\sigma_\text{I}/2\), and the two numbers would differ. Observe also that the largest principal stress, \(\sigma_1 = 134~\text{MPa}\), is smaller than the effective stress. A check based on the largest normal stress alone would have reported a factor of safety of \(1.9\) and been optimistic by a quarter. Combined stress states have to be reduced through a failure criterion, not by picking the largest component.

Example 3: Principal stresses of a three-dimensional state

A finite element model of a forged steel part reports, at a point in its interior, the stress components

\[ \sigma_x = 90, \quad \sigma_y = -30, \quad \sigma_z = 40, \quad \tau_{xy} = 40, \quad \tau_{yz} = -20, \quad \tau_{xz} = 30~\text{MPa} \]

shown in Figure 8.7.13. Determine the principal stresses \(\sigma_1 \ge \sigma_2 \ge \sigma_3\) and their directions, the maximum shear stress and the von Mises stress.

All three shear stresses are non-zero, so no principal direction can be read off and the full cubic 8.7.15 has to be solved. Its coefficients follow from 8.7.16, 8.7.17 and 8.7.18,

\[ \begin{aligned} I_1 &= 90 - 30 + 40~\text{MPa} = 100~\text{MPa}\\[4pt] I_2 &= 90\cdot(-30) + (-30)\cdot 40 + 40\cdot 90 - 40^2 - (-20)^2 - 30^2~\text{MPa}^2 = -3200~\text{MPa}^2\\[4pt] I_3 &= 90\cdot(-30)\cdot 40 + 2\cdot 40\cdot(-20)\cdot 30 - 90\cdot(-20)^2 - (-30)\cdot 30^2 - 40\cdot 40^2~\text{MPa}^3\\ &= -229\,000~\text{MPa}^3 \end{aligned} \]

so that, with \(\lambda\) in MPa, the characteristic equation is

\[ \lambda^3 - 100\lambda^2 - 3200\lambda + 229\,000 = 0 \]

SymPy computes the invariants from the components, as written above, and finds the three roots numerically.

s_x, s_y, s_z = 90, -30, 40                 # MPa
t_xy, t_yz, t_xz = 40, -20, 30

I_1 = s_x + s_y + s_z                                                        # @eq-invariant-1
I_2 = s_x*s_y + s_y*s_z + s_z*s_x - t_xy**2 - t_yz**2 - t_xz**2              # @eq-invariant-2
I_3 = (s_x*s_y*s_z + 2*t_xy*t_yz*t_xz
       - s_x*t_yz**2 - s_y*t_xz**2 - s_z*t_xy**2)                            # @eq-invariant-3
cubic = lam**3 - I_1*lam**2 + I_2*lam - I_3

# the roots are irrational, so we let nroots find them and put them in the order sigma_1 >= ...
roots = sorted(sp.Poly(cubic, lam).nroots(), reverse=True)

ltx(r"I_1 &=", I_1, r"~\text{MPa}, \quad I_2 =", I_2, r"~\text{MPa}^2, \quad I_3 =", I_3,
    r"~\text{MPa}^3",
    r"\\ \sigma_1, \sigma_2, \sigma_3 &=", float(roots[0]), r",\;", float(roots[1]), r",\;",
    float(roots[2]), r"~\text{MPa}",
    aligned=True)

The roots are \(\sigma_1 = 110.2~\text{MPa}\), \(\sigma_2 = 40.8~\text{MPa}\) and \(\sigma_3 = -51.0~\text{MPa}\). One is compressive and two are tensile, and the intermediate one is far from zero, so that this state is not close to plane stress in any frame. Figure 8.7.14 shows the cubic crossing zero at the three roots.

lam_range = np.linspace(-80, 140, 441)
p = sp.lambdify(lam, cubic, 'numpy')

fig, ax = plt.subplots(figsize=(6.5, 3.8))
ax.plot(lam_range, p(lam_range)/1e5, color='C0', lw=1.8)
ax.axhline(0, color='0.4', lw=0.8)
for k, r in enumerate(roots, start=1):
    ax.plot(float(r), 0, 'o', color='#c00000', ms=6, zorder=5)
    ax.annotate(rf'$\sigma_{k}$', (float(r), 0), xytext=(float(r) + 3, 0.5), fontsize=11,
                color='#c00000')
ax.set_xlim(-80, 140)
ax.set_ylim(-6, 5)
ax.set_xlabel(r'$\lambda$ [MPa]')
ax.set_ylabel(r'$\lambda^3 - I_1\lambda^2 + I_2\lambda - I_3$ [$10^5$ MPa$^3$]')
ax.grid(True, alpha=0.3)
end_ticks(fig)
plt.show()

The principal directions are the eigenvectors of 8.7.10. We collect them as the columns of the matrix \(\bm N = [\bm n_1\;\bm n_2\;\bm n_3]\), which numpy.linalg.eigh returns, with the eigenvalues, in ascending order; reversing the columns puts them in the order of 8.7.20. In the principal frame the tensor is diagonal,

\[ \bm N^\mathsf{T}\bm\sigma\bm N = \begin{bmatrix}\sigma_1 & 0 & 0\\ 0 & \sigma_2 & 0\\ 0 & 0 & \sigma_3\end{bmatrix} \]

and the directions are orthonormal, \(\bm N^\mathsf{T}\bm N = \bm I\). The eigenvalues must also agree with the roots of the cubic and reproduce the invariants by 8.7.19.

Sigma_ex = np.array([[s_x,  t_xy, t_xz],
                     [t_xy, s_y,  t_yz],
                     [t_xz, t_yz, s_z]], dtype=float)

w, V = np.linalg.eigh(Sigma_ex)            # ascending eigenvalues, orthonormal columns
w, N = w[::-1], V[:, ::-1]                 # reversed: sigma_1 >= sigma_2 >= sigma_3
s1, s2, s3 = w

ltx(r"\bm N = [\bm n_1\;\bm n_2\;\bm n_3] &=", N,
    r"\\ \bm N^\mathsf{T}\bm\sigma\bm N &=", np.round(N.T @ Sigma_ex @ N, 9) + 0.0,
    r"~\text{MPa}, \qquad \bm N^\mathsf{T}\bm N =", np.round(N.T @ N, 12) + 0.0,
    aligned=True, precision=4)

ltx(r"\max_k|\sigma_k^\text{eigh} - \sigma_k^\text{cubic}| &=",
    float(max(abs(a - float(b)) for a, b in zip(w, roots))), r"~\text{MPa}",
    r"\\ \sigma_1 + \sigma_2 + \sigma_3 &=", s1 + s2 + s3, r"\qquad I_1 =", I_1,
    r"\\ \sigma_1\sigma_2 + \sigma_2\sigma_3 + \sigma_3\sigma_1 &=", s1*s2 + s2*s3 + s3*s1,
    r"\qquad I_2 =", I_2,
    r"\\ \sigma_1\sigma_2\sigma_3 &=", s1*s2*s3, r"\qquad I_3 =", I_3,
    aligned=True, precision=4)

The eigenvalues agree with the roots to round-off, and the three invariants are reproduced. The first column of \(\bm N\) is the direction of the largest tensile stress, \(\bm n_1 \approx (0.92, 0.22, 0.33)\) up to sign, close to \(x\) and tipped toward \(z\), where \(\sigma_x\) and \(\sigma_z\) are the two largest normal components and \(\tau_{xz}\) couples them.

The maximum shear stress is 8.7.14, and the von Mises stress is 8.7.21,

\[ \begin{aligned} \tau_{\max} &= \frac{\sigma_1 - \sigma_3}{2}\\[4pt] \sigma_{\text{vM}} &= \sqrt{\tfrac{1}{2}\left[(\sigma_1 - \sigma_2)^2 + (\sigma_2 - \sigma_3)^2 + (\sigma_3 - \sigma_1)^2\right]} \end{aligned} \]

and 8.7.22 must give the same von Mises stress from the components, without the principal stresses.

tau_max = (s1 - s3)/2                                                        # @eq-max-shear-stress
sigma_vm = np.sqrt(((s1 - s2)**2 + (s2 - s3)**2 + (s3 - s1)**2)/2)          # @eq-von-mises-principal
sigma_vm_xyz = np.sqrt(((s_x - s_y)**2 + (s_y - s_z)**2 + (s_z - s_x)**2
                        + 6*(t_xy**2 + t_yz**2 + t_xz**2))/2)                # @eq-von-mises-general

ltx(r"\tau_{\max} &=", tau_max, r"~\text{MPa}",
    r"\\ \sigma_{\text{vM}} &=", sigma_vm, r"~\text{MPa}, \qquad \text{from the components: }",
    sigma_vm_xyz, r"~\text{MPa}",
    aligned=True, precision=4)

The largest shear stress is \(80.6~\text{MPa}\), on the planes that bisect the angle between \(\bm n_1\) and \(\bm n_3\), and the von Mises stress is \(140~\text{MPa}\) by both routes. The intermediate principal stress \(\sigma_2 = 40.8~\text{MPa}\) does not enter \(\tau_{\max}\) but does enter \(\sigma_{\text{vM}}\), and on this point the two yield criteria of Chapter 8.10 part ways.

What comes next

We now have a complete kinematic description in 8.7.3, a complete equilibrium statement in 8.7.8, and a boundary condition in 8.7.9. What is missing is the link between stress and strain in more than one dimension. Hooke’s law \(\sigma = E\varepsilon\) relates one component to one component; the generalisation to tensors is the subject of Chapter 8.9, and with it in hand the continuum model is closed and can be specialised to any geometry we like.