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.

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}\]

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)

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.14}\]

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.15}\]

The two expressions give the same number, since 8.7.14 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.15 simplifies to

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

which is the form that appears in most hand calculations, and the quantity that a finite element post-processor plots by default.

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}\). 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?

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

eig = S_num.eigenvects()
principal = sorted([sp.nsimplify(e[0]) for e in eig], reverse=True)
directions = [e[2][0].normalized() for e in sorted(eig, key=lambda e: -e[0])]

ltx(r"\sigma_1 &=", float(principal[0]), r"~\text{MPa}, \qquad \bm n_1 =", directions[0].evalf(3),
    r"\\ \sigma_2 &=", float(principal[1]), r"~\text{MPa}, \qquad \bm n_2 =", directions[1].evalf(3),
    aligned=True)

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

The maximum shear stress in the plane is half the difference of the principal stresses, a result that follows from 8.4.2 applied to the principal frame. A plane stress state has a third principal stress of zero out of the plane, and the largest shear stress over all planes is half the difference between the largest and the smallest of the three. The von Mises stress follows from 8.7.16.

tau_max_plane = (principal[0] - principal[1])/2
principal_3d = sorted([principal[0], principal[1], sp.Integer(0)], reverse=True)
tau_max_abs = (principal_3d[0] - principal_3d[-1])/2

sx, sy, txy = 120, -40, 50
sigma_vm = sp.sqrt(sx**2 - sx*sy + sy**2 + 3*txy**2)

ltx(r"\tau_{\max,\text{plane}} &=", float(tau_max_plane), r"~\text{MPa}",
    r"\\ \tau_{\max,\text{abs}} &=", float(tau_max_abs), r"~\text{MPa}",
    r"\\ \sigma_{\text{vM}} &=", float(sigma_vm), 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 two shear values coincide here because \(\sigma_2\) is negative, so the in-plane pair already spans the largest and smallest of the three principal stresses. Had both in-plane principal stresses been positive, the out-of-plane zero would have set the absolute maximum and the two numbers would differ. Observe also that the largest principal stress, \(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.

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.8, and with it in hand the continuum model is closed and can be specialised to any geometry we like.