332: Mechanical Behavior of Materials
Department of Materials Science and EngineeringNorthwestern University
1 Catalog Description
Plastic deformation and fracture of metals, ceramics, and polymeric materials; structure/property relations. Role of imperfections, state of stress, temperatures, strain rate. Lectures, laboratory. Prerequisites: 316 1; 316 2 (may be taken concurrently).
2 Course Outcomes
At the conclusion of the course students will be able to:
 Apply basic concepts of linear elasticity, including multiaxial stressstrain relationships through elastic constants for single and polycrystals.
 Quantify the different strengthening mechanisms in crystalline materials, based on interactions between dislocations and obstacles, such as: point defect (solid solution strengthening), dislocations (work hardening), grain boundaries (boundary strengthening) and particles (precipitation and dispersion strengthening).
 Apply fracture mechanics concepts to determine quantitatively when existing cracks in a material will grow.
 Describe how composite toughening mechanisms operate in ceramic matrix and polymer matrix composties.
 Derive simple relationships for the composite stiffness and strength based on those of the constituent phases.
 Exhibit a quantitative understanding of high temperature deformation in metals and ceramics, based on various creep mechanisms relate to diffusional and dislocation flow (Coble, NabarroHerring and Dislocation creep/climb mechanisms).
 Exhibit a basic understanding of factors affecting fatigue in engineering materials, as related to crack nucleation and propagation, as well as their connection to macroscopic fatigue phenomena.
 Describe the interplay between surface phenomena (environmental attack) and stresses leading to material embrittlement.
 Use the finite element method to calculate the stress and strain states for simple test cases, including a cantilever beam and a material with a circulat hole that is placed in tension.
 Use complex moduli to solve mechanics problems involving an oscillatory stress.
 Prepare and characterize specimens for measurement of mechanical properties.
 Write results from a laboratory project in the form of a journal article, and present their work orally as would be required in a technical forum.
 Select materials based on design requirements.
3 Introduction
In Figure
3.1, we show the materials available to mankind from the beginning of time until the current day.[
3] The data are plotted as socalled 'Ashby plots', where two material properties are chosen for the x and y axes. In Figure
3.1 we use density to separate different materials along the x axis, and the fracture strength,
${\sigma}_{f}$, to separate materials on the y axis. When the story begins, all we had to work with are the materials we could find around us or dig out of the ground. Over time, we Figure d out how to produce more materials for all kinds of purposes. The biggest developments come after World War II, a situation that coincides with the emergence of materials science as a discipline. We no longer rely on empiricism for materials development, but can actually begin to design new materials with the properties we want, or at least design materials that re better than anything we had before. The principles underlying the development of materials with the mechanical response forms the basis for this course.
4 Stress and Strain
The mechanical properties of a material are defined in terms of the strain response of material after a certain stress is applied. In order to properly understand mechanical properties, we have to have a good understanding of stress and strain, so that's where we begin.
Some Notes on Notation:
There are different ways to represent scalar quantities, vectors and matrices. Here's how we do it in this text:
 Scalar quantities are straight up symbols, like ${P}_{1}$, ${\sigma}_{12}$, etc.
 Vectors are indicated with an arrow over the symbol, like $\overrightarrow{P}$.
 Unit vectors are indicated with a caret above the symbol, like $\stackrel{\u02c6}{n}\mathrm{.}$
 Matrices are enclosed in square brackets, like $\left[\sigma \right]$
4.1 Tensor Representation of Stress
 i: surface normal (i= x, y, z)
 j: direction of force (j=x, y, z)
$$\begin{array}{cc}\left[\sigma \right]=\left[\begin{array}{cc}{\sigma}_{xx}& {\sigma}_{xy}\\ {\sigma}_{yx}& {\sigma}_{yy}\end{array}\right]& (4.1)\end{array}$$
 2 normal stresses, ${\sigma}_{xx}$, ${\sigma}_{yy}$. These are referred to as 'normal' stresses because the force acts perpendicular to the plane that it is referred to.
 A single shear stress, ${\sigma}_{xy}$.
In three dimensions we add a $z$ axis to the existing $x$ and $y$ axes, so the stress state is defined by a symmetric 3x3 tensor. The full stress tensor can be used to define the stresses acting on any given plane. To simplify the notation a bit we label the three orthogonal directions by numbers (1, 2 and 3) instead of letters ($x$, $y$ and $z$). The stress tensor gives the components of the force (${P}_{1},\phantom{\rule{6px}{0ex}}{P}_{2}$ and ${P}_{3}$) acting on a given plane. The plane is specified by the orientation of the unit vector, $\stackrel{\u02c6}{n}$ that is perpendicular to the plane. This vector has components ${n}_{1}$, ${n}_{2}$ and ${n}_{3}$ in the 1, 2 and 3 directions, respectively. It's a 'unit' vector because the length of the vector is 1, i.e. ${\left({n}_{1}^{2}+{n}_{2}^{2}+{n}_{3}^{2}\right)}^{1/2}=1.$ The relationship between $\overrightarrow{P}$, $\sigma $ and $\overrightarrow{n}$ is as follows:
or in more compact matrix notation:
$$\begin{array}{cc}\overrightarrow{P}=A\left[\sigma \right]\stackrel{\u02c6}{n}& (4.3)\end{array}$$
In graphical form the relationship is as shown in Figure 4.2. Like the 2dimensional stress tensor mentioned above, the 3dimensional stress tensor must also be symmetric in order for static equilibrium to be achieved. There are therefore 6 independent components of the threedimensional stress tensor:
 3 normal stresses, ${\sigma}_{11}$, ${\sigma}_{22}$ and ${\sigma}_{33}$, describing stresses applied perpendicular to the 1, 2 and 3 faces of the cubic volume element.
 3 shear stresses: ${\sigma}_{12}$, ${\sigma}_{23}$ and ${\sigma}_{13}$.
The three dimensional stress tensor is a 3x3 matrix with 9 elements (though only 6 are independent), corresponding to the three stress components acting on each of the three orthogonal faces of cube in the Cartesian coordinate system used to define the stress components. The 1 face has
${n}_{1}$=1,
${n}_{2}=0$ and
${n}_{3}=0.$ By setting
$\overrightarrow{n}=\left(1,0,0\right)$ in Eq.
4.2, we get the following for the stresses acting on the 1 face of the volume element:
$$\begin{array}{cc}\begin{array}{c}{P}_{1}/A={\sigma}_{11}\\ {P}_{2}/A={\sigma}_{12}\\ {P}_{3}/A={\sigma}_{13}\end{array}& (4.4)\end{array}$$Equivalent expressions can be obtained for the stresses acting on the 2 and 3 faces, by setting $\overrightarrow{n}=\left(0,1,0\right)$ and $\overrightarrow{n}=\left(0,0,1\right)$, respectively.
4.2 Tensor Transformation Law
The stress experienced by a material does not depend on the coordinate system used to define the stress state. The stress tensor will look very different if we chose a different set of coordinate axes to describe it, however, and it is important to understand how changing the coordinate system changes the stress tensor. We begin in this section by describing the procedure for obtaining the stress tensor that emerges from a given change in the coordinate system. We then describe the method for obtaining a specific set of coordinate axis which gives a diagonalized tensor where only normal stresses are present (Section
4.3).
In general, we consider the case where our 3 axes (which we refer to simply as axes 1, 2 and 3) are moved about the origin to define a new set of coordinate axes that we refer to as ${1}^{\prime}$, ${2}^{\prime}$ and ${3}^{\prime}$. As an example, consider the simple counterclockwise rotation around the 3 axis by an angle $\phi ,$ shown schematically in Figure 4.3. In general, the relative orientation of the transformed (rotated) and untransformed coordinate axes are given by a set of 9 angles between the 3 untransformed axes and the three transformed axes. In our notation we specify these angles as
${\theta}_{ij}$, where
$i$ specifies the transformed axes (
${1}^{\prime}$,
${2}^{\prime}$ or
${3}^{\prime}$) and
$j$ specifies the untransformed axis (1, 2 or 3). In our simple example, the angle between the 1 and
${1}^{\prime}$ axes is
$\phi $, so
${\theta}_{11}=\phi $. The angle between the 2 and
${2}^{\prime}$ axes is also
$\phi $, so
${\theta}_{22}=\phi $. The 3 axis remains unchanged in our rotation example, so
${\theta}_{33}=0.$ The
$3/{3}^{\prime}$ axis remains perpendicular to the 1,
${1}^{\prime}$ , 2,
${2}^{\prime}$ axes, so we have
${\theta}_{31}={\theta}_{32}={\theta}_{13}={\theta}_{23}=9{0}^{\u25cb}\mathrm{.}$ Finally, we see that the angle between the
${1}^{\prime}$ and the 2 axis is
$90\phi $ (
${\theta}_{12}=90\phi $) and the angle between the
${2}^{\prime}$ and 1 axis is
$90+\phi $ (
${\theta}_{21}=90+\phi $). The full
$\left[\theta \right]$ matrix in this case is as follows:
$$\begin{array}{cc}\left[\theta \right]=\left[\begin{array}{ccc}\phi & 90\phi & 90\\ 90+\phi & \phi & 90\\ 90& 90& 0\end{array}\right]& (4.5)\end{array}$$
Note that the $\left[\theta \right]$ matrix is NOT symmetric (${\theta}_{ij}\ne {\theta}_{ji}$), so you always need to make sure the first index, $i,$ (denoting the row in the $\left[\theta \right]$ matrix) corresponds to the transformed axes, and the second index, $j$ (denoting the column in the $\left[\theta \right]$ matrix) corresponds to the original, untransformed axes.
4.2.2 Expressions for the Stress Components
Once we specify all the different components of $\left[\theta \right]$, we can use the following general expression to obtain the stresses in the new (primed) coordinate system as a function of the stresses in the original coordinate system:
For each component of the stress tensor, we have to sum 9 individual terms (all combinations of $k$ and $l$ from 1 to 3). For example, ${\sigma}_{12}^{\prime}$ is given as follows:
The calculation is breathtakingly tedious if we do it all by hand, so it makes sense to automate this and do the calculation via computer, in our case with MATLAB or Python. In this example we'll start with a simple stress state corresponding to uniaxial extension in the 1 direction, with the following untransformed stress tensor:
%rotate45.m MATLAB version
sig=zeros(3); % create stress tensor and set to zero
sig(1,1)=5e6; % this is the only nonzero component
sigp=zeros(3); % initalize rotated streses to zero
phi=45;
theta=[phi,90phi,90;90+phi,phi,90;90,90,0];
for i=1:3
for j=1:3
for k=1:3
for l=1:3
sigp(i,j)=sigp(i,j)+cosd(theta(i,k))*cosd(theta(j,l))*sig(k,l);
end % cosd takes degrees instead of radians as the unit
end
end
end
sigp %display the transformed tensor components
The Python version of the code looks like this:
#!/usr/bin/env python3
# * coding: utf8 *
import numpy as np
sig=np.zeros((3, 3)) #% create stress tensor and set to zero
sig[0, 0] = 5e6; # this is the only nonzero component
sigp=np.zeros((3, 3)) # initalize rotated streses to zero
phi = 45
theta = [[phi,90phi,90], [90+phi,phi,90], [90,90,0]]
theta = np.deg2rad(theta) # trig functions need angles in radians
for i in [0, 1, 2]:
for j in [0, 1, 2]:
for k in [0, 1, 2]:
for l in [0, 1, 2]:
sigp[i,j]=sigp[i,j]+np.cos(theta[i,k])*np.cos(theta[j,l])*sig[k,l]
print(sigp) # display the transformed tensor components
Note the following:
 The normal stresses in the 1 and 2 directions are equal to one another.
 The transformed shear stress in the 12 plane is half the original tensile stress.
 The sum of the normal stresses (the sum of the diagonal components of the stress tensor) is unchanged by the coordinate transformation
4.2.3 An Easier Way: Transformation by Direct Matrix Multiplication
$$\begin{array}{cc}\left[Q\right]=cos\left[\theta \right]& (4.10)\end{array}$$
For the simple case of rotation about the z axis, the angles are given by Eq.
4.5, so that
$\left[Q\right]$ is given as:
$$\begin{array}{cc}\left[Q\right]=\left[\begin{array}{ccc}cos\phi & cos\left(90\phi \right)& cos90\\ cos\left(90+\phi \right)& cos\phi & cos90\\ cos90& cos90& cos0\end{array}\right]=\left[\begin{array}{ccc}cos\phi & sin\phi & 0\\ sin\phi & cos\phi & 0\\ 0& 0& 1\end{array}\right]& (4.11)\end{array}$$
The transformed stress is now obtained by the following simple matrix multiplication:
where the${\left[Q\right]}^{T}$ is the transpose of$\left[Q\right]$:
$$\begin{array}{cc}{Q}^{T}\left(i,j\right)=Q\left(j,i\right)& (4.13)\end{array}$$
For the rotation by $\phi $ around the z axis, ${\left[Q\right]}^{T}$ is given by the following:
$$\begin{array}{cc}{\left[Q\right]}^{T}=\left[\begin{array}{ccc}cos\phi & sin\phi & 0\\ sin\phi & cos\phi & 0\\ 0& 0& 1\end{array}\right]& (4.14)\end{array}$$
Equation 4.12 is much easier to deal with than Eq.
4.7, at least in MATLAB, where matrices can be multiplied as regular numbers. The MATLAB code to take a uniaxial stress state and rotate the coordinate system by 45
${}^{\u25cb}$ about the 3 axis looks like this if we base it on Eq.
4.12:
%rotate45_matrix.m
sig=zeros(3); % create stress tensor and set to zero
sig(1,1)=5e6; % this is the only nonzero component
phi=45;
theta=[phi,90phi,90;90+phi,phi,90;90,90,0];
Q=cosd(theta);
QT=transpose(Q);
sigp=Q*sig*QT
Running this script gives the output shown in Figure
4.4,
i.e. we obtain exactly the same result we obtained by using Eq.
4.7.
4.3 Principal Stresses
Note that we still need 6 independent parameters to specify a stress state: the 3 principal stresses, in addition to three parameters that specify the orientation of the principal axes. The stress tensor depends on our definition of the axes, but it is always possible to chose the axes so that all of the shear components of the stress tensor vanish, so that the stress tensor looks like the following:
The following MATLAB script results in the output shown in Figure
4.6.
%rotatesym.m
clear all
syms phi; % declare as symbolic variable
sig=sym(zeros(3)); % create sybolic stress tensor and set to zero
syms sig1p sig2p sig3p % these are the normal stresses
sig(1,1)=sig1p; % put normal stresses into the tensor
sig(2,2)=sig2p;
sig(3,3)=sig3p;
theta=[phi,pi/2+phi,pi/2;pi/2phi,phi,pi/2;pi/2,pi/2,0];
Q=cos(theta);
QT=transpose(Q);
sigp=Q*sig*QT;
sigp=simplify(sigp);
pretty(sigp);
The components of the stress tensor in the rotated coordinate system are:
$$\begin{array}{cc}\left[{\sigma}^{\prime}\right]=\left[\begin{array}{ccc}{\sigma}_{1}^{p}{\sigma}_{1}^{p}{sin}^{2}\phi +{\sigma}_{2}^{p}{sin}^{2}\phi & sin\left(2\phi \right)\left({\sigma}_{1}^{p}{\sigma}_{2}^{p}\right)/2& 0\\ sin\left(2\phi \right)\left({\sigma}_{1}^{p}{\sigma}_{2}^{p}\right)/2& {\sigma}_{2}^{p}+{\sigma}_{1}^{p}{sin}^{2}\phi {\sigma}_{2}^{p}{sin}^{2}\phi & 0\\ 0& 0& 0\end{array}\right]& (4.16)\end{array}$$
4.3.1 Mohr's Circle Construction
The Mohr circle is a graphical construction that can be used to describe a two dimensional stress state. A two dimensional stress state is specified by two principal stresses, ${\sigma}_{1}^{p}$ and ${\sigma}_{2}^{p}$, and by the orientation of the principal axes. The Mohr circle is drawn with a diameter of ${\sigma}_{1}^{p}{\sigma}_{2}^{p}$, centered at $\left({\sigma}_{1}^{p}+{\sigma}_{2}^{p}\right)$. We can use MATLAB's symbolic math capability to derive the Mohr circle construction as follows:
% Mohr's circle derivation  beginning from previous calculation of sigp
center=(sig1p+sig2p)/2; % the center of Mohr's circle
radius=(sig1psig2p)/2; % radius of Mohr's circle
% now we shift the origin to the center of Mohr's circle, and
% normalize by the radius of Mohrs circle
sigp(1,1)=(sigp(1,1)center)/radius;
sigp(2,2)=(sigp(2,2)center)/radius;
% shear stresses are plotted on the y axis, so we don't need to shift the
% origin, we just need to normlized by the radius
sigp(1,2)=sigp(1,2)/radius;
sigp(2,1)=sigp(1,2);
simplify(sigp)
In general, there are three principal stresses,
${\sigma}_{1}^{p}$,
${\sigma}_{2}^{p}$ and
${\sigma}_{3}^{p}$, and we can draw the Mohr circle construction with any combination of these 3 stresses. We end up with 3 different circles, as shown in Figure
4.9. Note that the convention is that
${\sigma}_{1}^{p}$ is the largest principal stress and that
${\sigma}_{3}^{p}$ is the smallest principal stress,
i.e. ${\sigma}_{1}^{p}>{\sigma}_{2}^{p}>{\sigma}_{3}^{p}$. An important result is that the largest shear stress,
${\tau}_{max}$, is given by the difference between the largest principal stress and the smallest one:
$$\begin{array}{cc}{\tau}_{max}=\frac{1}{2}\left({\sigma}_{1}^{p}{\sigma}_{3}^{p}\right)& (4.17)\end{array}$$
This maximum shear stress is an important quantity, because it determines when a material will deform plastically (much more on this later). In order to determine this maximum shear stress, we need to first Figure out what the principal stresses are. In some cases this is easy. In a uniaxial tensile test, one of the principals stresses is the applied stress, and the other two principal stresses are equal to zero.
The individual Mohr's circles in Figure 4.9 correspond to rotations in the around the individual principal axes. Circle
${C}_{1}$ corresponds to rotation around the direction in which
${\sigma}_{1}^{p}$ is directed,
${C}_{2}$ corresponds to rotation around the direction in which
${\sigma}_{2}^{p}$ is directed, and
${C}_{3}$ corresponds to rotation around the direction in which
${\sigma}_{3}^{p}$ is directed. A consequence of this is that is always possible to use the Mohr's circle construction to determine the principal stresses if there is only one nonzero shear stress. Consider, for example, the following stress state:
In this case there is only one nonzero shear stress (${\sigma}_{13}$), so we can determine the principals stresses in the following manner:
 One of the three principal stresses is the normal stress in the direction that does not involve either of the directions in the nonzero shear stress. Since the nonzero shear stress in our example is ${\sigma}_{13},$ one of the principal stresses is ${\sigma}_{22}$=3 MPa.
 Draw a Mohr's circle construction using the two normal stresses and the nonzero shear stress, in this case ${\sigma}_{11}$, ${\sigma}_{33}$ and ${\sigma}_{13}$. Mohr's circle is centered at the the average of these two normal stresses, in our case at $\left({\sigma}_{11}+{\sigma}_{33}\right)/2$ = 4 MPa. We refer to this stress as ${\sigma}_{c}$, as illustrated below in Figure 4.10.
 Determine the radius of the circle, ${R}_{M}$, is given as:
$${R}_{m}=\sqrt{{\left({\sigma}_{33}{\sigma}_{c}\right)}^{2}+{\sigma}_{13}^{2}}=\sqrt{{\left(54\right)}^{2}+{\sigma}_{13}^{2}}=2.24\phantom{\rule{6px}{0ex}}MPa$$
 The principal stresses are given by the intersections of the circle with the horizontal axis:
$${\sigma}^{p}={\sigma}_{c}\pm {R}_{m}=6.24\phantom{\rule{6px}{0ex}}MPa,\phantom{\rule{6px}{0ex}}1.76\phantom{\rule{6px}{0ex}}MPa$$
The third principal stress is 3 MPa, as we already determined.
 The maximum shear stress is half the difference between the largest principal stress (6.64 MPa) and the smallest one (1.76), or 2.24 MPa.
4.3.2 Critical Resolved Shear Stress for Uniaxial Tension
$$\begin{array}{cc}{\tau}_{rss}={\sigma}_{1}^{p}cos\phi cos\left(9{0}^{\u25cb}\phi \right)& (4.20)\end{array}$$
We can use the identities $cos\left(90\phi \right)=sin\phi $ and $sin\left(2\phi \right)=2sin\phi cos\phi $ to obtain the following:
$$\begin{array}{cc}{\tau}_{rss}=\frac{{\sigma}_{1}^{p}}{2}sin\left(2\phi \right)& (4.21)\end{array}$$
We can get the same thing from the Mohr's circle construction to redefine the axes by a rotation of $\phi $. The shear stress is simply the radius of the circle (${\phi}_{1}^{p}/2$ in this case) multiplied by $sin\left(2\phi \right)$. Mohr's circle also gives us the normal stresses:
$$\begin{array}{cc}\begin{array}{c}{\sigma}_{11}=\frac{{\sigma}_{1}^{p}}{2}+\frac{{\sigma}_{1}^{p}}{2}cos\left(2\phi \right)\\ {\sigma}_{22}=\frac{{\sigma}_{1}^{p}}{2}\frac{{\sigma}_{1}^{p}}{2}cos\left(2\phi \right)\end{array}& (4.22)\end{array}$$
The untransformed 2dimensional stress tensor looks like this:
$$\begin{array}{cc}\left[\sigma \right]=\left[\begin{array}{cc}{\sigma}_{1}^{p}& 0\\ 0& 0\end{array}\right]& (4.23)\end{array}$$
The transformed stress tensor (after rotation by $\phi $ to give the resolved shear stress) looks like this:
$$\begin{array}{cc}\left[{\sigma}^{\prime}\right]=\left[\begin{array}{cc}\frac{{\sigma}_{1}^{p}}{2}+\frac{{\sigma}_{1}^{p}}{2}cos\left(2\phi \right)& \frac{{\sigma}_{1}^{p}}{2}sin\left(2\phi \right)\\ \frac{{\sigma}_{1}^{p}}{2}sin\left(2\phi \right)& \frac{{\sigma}_{1}^{p}}{2}\frac{{\sigma}_{1}^{p}}{2}cos\left(2\phi \right)\end{array}\right]& (4.24)\end{array}$$
To illustrate, we'll start with the stress state specified by Eq.
4.9. which we got by starting with a simple uniaxial extension in the 1 direction, and rotating the coordinate system by 45
${}^{\u25cb}$ about the 3 axis. The MATLAB script to do this is very simple and is as follows:
Here's the output generated by this script:
3 principal axes returned as column vectors. In this case there is a single normal stress, acting in a direction midway between the original x and y axes. The original uniaxial stress state is recovered in this example, as it should be. To summarize:
 Principal Stresses: Eigenvalues of the stress tensor
 Principal Stress directions: Eigenvectors of the stress tensor
4.3.4 Stress Invariants
$$\begin{array}{cc}p=\frac{1}{3}\left({\sigma}_{11}+{\sigma}_{22}+{\sigma}_{33}\right)=\frac{1}{3}\left({\sigma}_{1}^{p}+{\sigma}_{2}^{p}+{\sigma}_{3}^{p}\right)& (4.25)\end{array}$$
$$\begin{array}{cc}{I}_{1}={\sigma}_{11}+{\sigma}_{22}+{\sigma}_{33}& (4.26)\end{array}$$
The second and third stress invariants,${I}_{2}$ and${I}_{3}$, are also independent of the way the axes are defined:
$$\begin{array}{cc}{I}_{2}={\sigma}_{11}{\sigma}_{22}+{\sigma}_{22}{\sigma}_{33}+{\sigma}_{33}{\sigma}_{11}{\sigma}_{12}^{2}{\sigma}_{13}^{2}{\sigma}_{23}^{2}& (4.27)\end{array}$$
$$\begin{array}{cc}{I}_{3}={\sigma}_{11}{\sigma}_{22}{\sigma}_{33}{\sigma}_{11}{\sigma}_{23}^{2}{\sigma}_{22}{\sigma}_{13}^{2}{\sigma}_{33}{\sigma}_{12}^{2}+2{\sigma}_{12}{\sigma}_{13}{\sigma}_{23}& (4.28)\end{array}$$
The fact that these invariants exist is useful, but we're not going to do a lot more with them in this class. The most important invariant is the hydrostatic pressure, $p$.
There are 3 related definitions of the strain:
 Engineering strain
 Tensor strain
 Generalized strain (large deformations, also referred to as 'finite strain')
 ${P}_{1}$ moves by an amount$\left(u,\phantom{\rule{6px}{0ex}}v,\phantom{\rule{6px}{0ex}}w\right)$
 ${P}_{2}$ moves by$\left(u+du,v+dv,w+dw\right)$
Strain describes how much farther point to moved in three different directions, as a function of how far ${P}_{2}$ was from ${P}_{1}$ initially. For small strains we can ignore higher order terms in a Taylor expansion for $du$, $dv$ and $dw$ and maintain only the first, partial derivative terms as follows:
$$\begin{array}{cc}\begin{array}{c}du=\frac{\partial u}{\partial x}dx+\frac{\partial u}{\partial y}dy+\frac{\partial u}{\partial z}dz\\ dv=\frac{\partial v}{\partial x}dx+\frac{\partial v}{\partial y}dy+\frac{\partial v}{\partial z}dz\\ dw=\frac{\partial w}{\partial x}dx+\frac{\partial w}{\partial y}dy+\frac{\partial w}{\partial z}dz\end{array}& (4.29)\end{array}$$
The three normal components of the strain correspond to the change in the displacement in a given direction corresponds to a change in initial separation between the points of interest in the same direction:
The engineering shear strains are defined as follows:
Note: shear strains are often represented by the lower case Greek gamma to distinguish them from normal strains:
$$\begin{array}{cc}{\gamma}_{xy}\equiv {e}_{xy};\phantom{\rule{6px}{0ex}}{\gamma}_{yz}\equiv {e}_{yz};\phantom{\rule{6px}{0ex}}{\gamma}_{xz}={e}_{xz}& (4.32)\end{array}$$
4.4.2 Tensor Shear Strains
Engineering strains relate two vectors to one another, ($dx,\phantom{\rule{6px}{0ex}}dy,\phantom{\rule{6px}{0ex}}dz$) and ($du,\phantom{\rule{6px}{0ex}}dv,\phantom{\rule{6px}{0ex}}dw$), just as a tensor does, but the the transformation law between different coordinate systems is not obeyed for the engineering strains. For this reason the engineering strains are NOT tensor strains. Fortunately, all we need to do to change engineering strains to tensor strains is to divide the shear components by 2. In our notation we use $e$ to indicate engineering strain and $\epsilon $ to indicate tensor strains. The tensor normal strains are exactly the same as the engineering normal strains:
$$\begin{array}{cc}\begin{array}{c}{\epsilon}_{xx}={e}_{xx}\\ {\epsilon}_{yy}={e}_{yy}\\ {\epsilon}_{zz}={e}_{zz}\end{array}& (4.33)\end{array}$$
Engineering shear strains (${e}_{yz},\phantom{\rule{6px}{0ex}}{e}_{xz},\phantom{\rule{6px}{0ex}}{e}_{zy}$) are divided by two to give tensor shear strains:
Note that the tensor strains must be used in coordinate transformations (axis rotation, calculation of principal strains, ${e}_{1}^{p}$, ${e}_{2}^{p}$, ${e}_{3}^{p}$).
4.4.3 Generalized Strain
$$\begin{array}{cc}\begin{array}{c}{\lambda}_{1}=1+{e}_{1}^{P}\\ {\lambda}_{2}=1+{e}_{2}^{P}\\ {\lambda}_{3}=1+{e}_{3}^{P}\end{array}& (4.35)\end{array}$$
$$\begin{array}{cc}{e}_{1}^{t}=ln\left({\lambda}_{1}\right)& (4.36)\end{array}$$
This expression for the true strain can be obtained by recognizing that the incremental strain is always given by the fractional increase in length $d\ell /\ell $, where $\ell $ is the current length of the material as it is being deformed. If the initial length is ${\ell}_{0}$ and the final, deformed length is ${\ell}_{f}$, then the total true strain, ${e}_{1}^{true}$ is obtained by integrating the incremental strains accumulated throughout the entire deformation history:
The extension ratios provide a useful description of the strain for both small and large values of the strain. A material with isotropic mechanical properties has the same coordinate axes for the principal stresses and the principal strains.
Now that we've formally defined stress and strain we can give some specific examples where these definitions are used, and begin to define some elastic constants. We'll begin with the two most fundamental deformation states: simple shear and hydrostatic compression. These are complementary strain states  for an isotropic material simple shear changes the shape but not the volume, and hydrostatic compression changes the volume but not the shape. We'll eventually show that for an isotropic material there are only two independent elastic constants, so if we know how an isotropic material behaves in response to these two stress states, we have a complete understanding of the elastic properties of the material.
4.5.1 Simple Shear
The stress tensor looks like this:
$$\begin{array}{cc}\sigma =\left[\begin{array}{ccc}0& {\sigma}_{xy}& 0\\ {\sigma}_{xy}& 0& 0\\ 0& 0& 0\end{array}\right]& (4.38)\end{array}$$
$$\begin{array}{cc}{e}_{xy}=\frac{u}{d}& (4.39)\end{array}$$
We need to divide the engineering shear strains by 2 to get the tensor strains, so we get the following:
$$\begin{array}{cc}G\phantom{\rule{6px}{0ex}}\left(shear\phantom{\rule{6px}{0ex}}modulus\right)=\frac{{\sigma}_{xy}}{{e}_{xy}}& (4.41)\end{array}$$
Note that the volume of the material is not changed, but it's shape has. In very general terms we can view the shear modulus of a material as a measure of its resistance to a change in shape under conditions where the volume remains constant.
4.5.2 Simple Shear and the Mohr's Circle Construction for Strains
I
4.5.3 Hydrostatic Compression
$$\begin{array}{cc}{K}_{b}=V\frac{dp}{dV}& (4.43)\end{array}$$
The hydrostatic stress state corresponds to the stress state where there are no shear stresses, and each of the normal stresses are equal. Compressive stresses are defined as negative, whereas a compressive pressure is positive, so the stress state for hydrostatic compression looks like this:
$$\begin{array}{cc}\sigma =\left[\begin{array}{ccc}p& 0& 0\\ 0& p& 0\\ 0& 0& p\end{array}\right]& (4.44)\end{array}$$
4.5.4 Uniaxial Extension
Uniaxial extension corresponds to the application of a normal stress along one direction, which we define here as the 3 direction so that the stress tensor looks like this:
$$\begin{array}{cc}\sigma =\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& {\sigma}_{33}\end{array}\right]& (4.45)\end{array}$$
We can measure two separate strains from this experiment: the longitudinal strain in the same direction that we apply the stress, and the transverse strain, ${e}_{22}$, measured in the direction perpendicular to the applied stress (we'll assume that the sample is isotropic in the 12 plane, so ${e}_{11}={e}_{22}$. The strains are given by the fractional changes in the length and width of the sample:
$$\begin{array}{cc}{e}_{33}=\frac{\Delta \ell}{{\ell}_{0}}& (4.46)\end{array}$$
$$\begin{array}{cc}{e}_{22}=\frac{\Delta w}{w}& (4.47)\end{array}$$
$$\begin{array}{cc}\nu ={e}_{22}/{e}_{33}& (4.49)\end{array}$$
4.5.5 Longitudinal Compression
A final deformation state that we will consider is longitudinal compression. In this state all of the compression is in one direction, which we will specify as the 3 direction. The strains in the other two direction are constrained to be zero, so the strain state is as follows:
$$e=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& {e}_{33}\end{array}\right]$$
Note that the strain state is similar to that of uniaxial extension or compression (Figure 4.18), but in the current case we have a single nonzero strain instead of a single nonzero stress. Finite values of
${\sigma}_{11}$ and
${\sigma}_{22}$ must exist in order for sample in order for this strain state to be maintained, but we're not going to worry about those for now. Instead, we'll use the following relationship for the longitudinal elastic modulus,
${E}_{\ell}$ which is the ratio of
${\sigma}_{33}$ to
${e}_{33}$ for this strain state. Note that this deformation state changes both the shape and volume of the material, so
${E}_{\ell}$ involves both
$G$ and
${K}_{b}$:
$$\begin{array}{cc}{E}_{\ell}=\frac{{\sigma}_{33}}{{e}_{33}}={K}_{b}+\frac{4}{3}G& (4.50)\end{array}$$
4.6 Representative Moduli
A few typical values for
$G$ and
$K$ are listed in Table
elsewhere. Liquids do not have a shear modulus, but they do have a bulk modulus.
4.7 Case Study: Speed of Sound
$$\begin{array}{cc}{V}_{sound}=\sqrt{\frac{M}{\rho}}& (4.51)\end{array}$$
Here $\rho $ is the density of the material. The modulus that we need to use depends on the type of sound wave that is propagating. The two most common are a shear wave and a longitudinal compressional wave:
 Longitudinal compressional wave: $M={E}_{\ell}$
 Shear wave: $M=G$
In a liquid or gas (like air), $G=0$ and shear waves cannot propagate. In this case there is a single sound velocity obtained by setting $M={K}_{b}$. For an ideal gas:
$$\begin{array}{cc}P=\frac{n}{V}RT& (4.52)\end{array}$$
If the compression is applied slowly enough so that the temperature of the gas can equilibrate, we have:
5 Matrix representation of Stress and Strain
As usual, we begin by replacing the directions (x, y, and z) with numbers:
$x\to 1$,
$y\to 2$,
$z\to 3$. Once we do this we have 6 stress components, and six strain components. We then number these components from 16, so that 1, 2 and 3 are the normal components and 4, 5 and 6 are the shear components. We do this for both stress and strain as shown in Table
2.
A series of elastic constants relate the stresses to the strains. We can do calculations in either of the following two ways:
 Start with a column vector consisting of the 6 elements of an applied stress, and use the compliance matrix to calculate the strains.
 Start with a column vector consisting of the 6 elements of an applied strain, and use the stiffness matrix to calculate the stresses.
In each case we use a$6\times 6$ matrix to relate two 6element column vectors to one another. The procedure in each case is outlined below.
5.1 Compliance matrix
$$\begin{array}{cc}\left[\begin{array}{c}{e}_{1}\\ {e}_{2}\\ {e}_{3}\\ {e}_{4}\\ {e}_{5}\\ {e}_{6}\end{array}\right]=\left[\begin{array}{cccccc}{s}_{11}& {s}_{12}& {s}_{13}& {s}_{14}& {s}_{15}& {s}_{16}\\ {s}_{21}& {s}_{22}& {s}_{23}& {s}_{24}& {s}_{25}& {s}_{26}\\ {s}_{31}& {s}_{32}& {s}_{33}& {s}_{34}& {s}_{35}& {s}_{36}\\ {s}_{41}& {s}_{42}& {s}_{43}& {s}_{44}& {s}_{45}& {s}_{46}\\ {s}_{51}& {s}_{52}& {s}_{53}& {s}_{54}& {s}_{55}& {s}_{56}\\ {s}_{61}& {s}_{62}& {s}_{63}& {s}_{64}& {s}_{65}& {s}_{66}\end{array}\right]\left[\begin{array}{c}{\sigma}_{1}\\ {\sigma}_{2}\\ {\sigma}_{3}\\ {\sigma}_{4}\\ {\sigma}_{5}\\ {\sigma}_{6}\end{array}\right]& (5.1)\end{array}$$
The matrix must be symmetric, with ${s}_{ij}={s}_{ji}$, so there are a maximum of 21 independent compliance coefficients:
$$\begin{array}{cc}\left[\begin{array}{c}{e}_{1}\\ {e}_{2}\\ {e}_{3}\\ {e}_{4}\\ {e}_{5}\\ {e}_{6}\end{array}\right]=\left[\begin{array}{cccccc}{s}_{11}& {s}_{12}& {s}_{13}& {s}_{14}& {s}_{15}& {s}_{16}\\ {s}_{12}& {s}_{22}& {s}_{23}& {s}_{24}& {s}_{25}& {s}_{26}\\ {s}_{13}& {s}_{23}& {s}_{33}& {s}_{34}& {s}_{35}& {s}_{36}\\ {s}_{14}& {s}_{24}& {s}_{34}& {s}_{44}& {s}_{45}& {s}_{46}\\ {s}_{15}& {s}_{25}& {s}_{35}& {s}_{45}& {s}_{55}& {s}_{56}\\ {s}_{16}& {s}_{26}& {s}_{36}& {s}_{46}& {s}_{56}& {s}_{66}\end{array}\right]\left[\begin{array}{c}{\sigma}_{1}\\ {\sigma}_{2}\\ {\sigma}_{3}\\ {\sigma}_{4}\\ {\sigma}_{5}\\ {\sigma}_{6}\end{array}\right]& (5.2)\end{array}$$
Note that the compliance coefficients have the units of an inverse stress (Pa${}^{1}$).
5.2 Stiffness Matrix
The stiffness matrix (c) is the inverse of compliance matrix (note the somewhat confusing notation in that the compliance matrix is $s$ and the stiffness is $c$, backwards from what you might expect). The stiffness coefficients have units of stress.$$\begin{array}{cc}\left[\begin{array}{c}{\sigma}_{1}\\ {\sigma}_{2}\\ {\sigma}_{3}\\ {\sigma}_{4}\\ {\sigma}_{5}\\ {\sigma}_{6}\end{array}\right]=\left[\begin{array}{cccccc}{c}_{11}& {c}_{12}& {c}_{13}& {c}_{14}& {c}_{15}& {c}_{16}\\ {c}_{12}& {c}_{22}& {c}_{23}& {c}_{24}& {c}_{25}& {c}_{26}\\ {c}_{13}& {c}_{23}& {c}_{33}& {c}_{34}& {c}_{35}& {c}_{36}\\ {c}_{14}& {c}_{24}& {c}_{34}& {c}_{44}& {c}_{45}& {c}_{46}\\ {c}_{15}& {c}_{25}& {c}_{35}& {c}_{45}& {c}_{55}& {c}_{56}\\ {c}_{16}& {c}_{26}& {c}_{36}& {c}_{46}& {c}_{56}& {c}_{66}\end{array}\right]\left[\begin{array}{c}{e}_{1}\\ {e}_{2}\\ {e}_{3}\\ {e}_{4}\\ {e}_{5}\\ {e}_{6}\end{array}\right]& (5.3)\end{array}$$
5.3 Symmetry requirements on the compliance (or stiffness) matrix.
5.3.1 Orthorhombic symmetry
 ${E}_{1}=1/{s}_{11}$, ${E}_{2}=1/{s}_{22}$ and ${E}_{3}=1/{s}_{33}$, Young's moduli for extension in the 1, 2 and 3 directions, respectively.
 ${G}_{1}=1/{s}_{44}$, ${G}_{2}=1/{s}_{55}$ and ${G}_{3}=1/{s}_{66}$, Shear moduli for shear in the planes perpendicular to the 1, 2 and 3 directions, respectively.
 ${s}_{12}$, ${s}_{13}$ and ${s}_{23}$, which relate stresses in one direction to strains in the perpendicular direction.
5.3.2 Fiber Symmetry
For a material with fiber symmetry, one of the axes is unique (in this case the 3 axis) and the material is isotropic in the orthogonal plane. Since the 1 and 2 axes are identical, there are now 5 independent elastic constants ${s}_{33}$, ${s}_{13}$, ${s}_{44}$, ${s}_{11}$, ${s}_{12}$:
Examples of materials with fiber symmetry include the following:
 Many liquid crystalline polymers (e.g. Kevlar).
 Materials after cold drawing (plastic deformation to high strains, carried out below the glass transition temperature or melting temperature of the material.)
To better understand the significance of the 5 elastic constants for fiber symmetry, it is useful to consider the types of experiment we would need to conduct to measure each of them for a cylindrical fiber. The necessary experiments are described below.
5.3.2.1 Fiber extension along 3 direction: measurement of ${s}_{33}$ and ${s}_{13}$
We then obtain
${s}_{33}$ and
${s}_{13}$ from Eq.
5.5, recalling that
${\sigma}_{3}$ is the only nonzero stress component in this situation:
$$\begin{array}{cc}\begin{array}{c}{s}_{33}={e}_{33}/{\sigma}_{3}\\ {s}_{13}={e}_{1}/{\sigma}_{3}\end{array}& (5.7)\end{array}$$
5.3.2.2 Fiber Torsion: Measurement of ${s}_{44}$
We define a cylindrical system with a z axis along the fiber axis. The other axes in this coordinate system are the distance $r$ from this axis of symmetry, and the angle $\theta $ around the z axis. The shear strain in the$\theta z$ plane depends only on$r$, and is given by :
$$\begin{array}{cc}{e}_{\theta z}=r\frac{d\theta}{dz}=r\frac{{\theta}_{0}}{\ell}& (5.8)\end{array}$$
The corresponding shear stress is obtained by multiplying by the shear modulus, $G$ characterizing deformation in the 12 and 23 planes:
$$\begin{array}{cc}{\sigma}_{\theta z}=Gr{\theta}_{0}/\ell & (5.9)\end{array}$$
We integrate the shear stress to give the torque, $T$:
So once we know the torsional stiffness of the fiber (the relationship between the applied $T$ and ${\theta}_{0}$$)$ we know the shear modulus, $G$. This shear modulus is simply the inverse of${s}_{44}$:
$$\begin{array}{cc}G=\frac{1}{{s}_{44}}& (5.11)\end{array}$$
5.3.3 Fiber compression in 12 plane: determination of ${s}_{11}$ and ${s}_{12}$
where $P$ is the force applied to the fiber, ${d}_{0}$ is its undeformed diameter and $\ell $ is its length.
$$\begin{array}{cc}{e}_{3}={s}_{13}\left({\sigma}_{1}^{p}+{\sigma}_{2}^{p}\right)+{s}_{33}{\sigma}_{3}& (5.13)\end{array}$$
Setting ${e}_{3}$ to zero in this equation gives the following for ${\sigma}_{3}$:
$$\begin{array}{cc}{\sigma}_{3}=\frac{{s}_{13}}{{s}_{33}}\left({\sigma}_{1}^{p}+{\sigma}_{2}^{p}\right)& (5.14)\end{array}$$
A consequence of this stress is that the frictionless expression for $b$ gets modified to the following:
$$\begin{array}{cc}{b}^{2}=\frac{2F{d}_{0}}{\ell}\left({s}_{11}\frac{{s}_{13}^{2}}{{s}_{33}}\right)& (5.15)\end{array}$$
The remaining constant, ${s}_{12}$, is determined from a measurement of $d/{d}_{0}$, the ratio of the fiber width at the midplane to the original width of the fiber. This relationship is complicated, and involves several of the different elastic constants.
$$\begin{array}{cc}\frac{d}{{d}_{0}}=f\left(P,\phantom{\rule{6px}{0ex}}{s}_{11},\phantom{\rule{6px}{0ex}}{s}_{13},\phantom{\rule{6px}{0ex}}{s}_{33},\phantom{\rule{6px}{0ex}}{s}_{12}\right)& (5.16)\end{array}$$
5.3.4 Cubic Symmetry
For a material with cubic symmetry the 1,2 and 3 axes are all identical to one another, and we end up with 3 independent elastic constants:
5.3.5 Isotropic systems
For an isotropic material all axes are equivalent, and the properties are invariant to any rotation of the coordinate axes. In this case there are two independent elastic constants, and the compliance matrix looks like this:
The requirement that the material properties be invariant with respect to any rotation of the coordinate axes results in the requirement that ${s}_{44}=2\left({s}_{11}{s}_{12}\right)$, so there are two independent elastic constants. The shear modulus, $G$, Young's modulus $E$ and Poisson's ratio, $\nu $ are given as follows:
$$\begin{array}{cc}G=1/2\left({s}_{11}{s}_{12}\right);\phantom{\rule{6px}{0ex}}\phantom{\rule{6px}{0ex}}E=1/{s}_{11};\phantom{\rule{6px}{0ex}}\phantom{\rule{6px}{0ex}}\nu ={s}_{12}/{s}_{11}& (5.19)\end{array}$$
5.3.5.1 Bulk Modulus for an Isotropic Material
The bulk modulus, ${K}_{b}$, of a material describes it's resistance to a change in volume (or density) when we apply a hydrostatic pressure, $p$. It is defined in the following way:
$$\begin{array}{cc}{K}_{b}=V\frac{dP}{dV}& (5.20)\end{array}$$
The stress state in this case is as follows:
From the compliance matrix (Eq. 5.18) we get
${e}_{1}={e}_{2}={e}_{3}=p\left({s}_{11}+2{s}_{12}\right)$. The change in volume,
$\Delta V$ can be written in terms of the three principal extension ratios,
${\lambda}_{1}$,
${\lambda}_{2}$ and
${\lambda}_{3}$:
$$\begin{array}{cc}\frac{\Delta V}{{V}_{0}}=\frac{V}{{V}_{0}}1={\lambda}_{1}{\lambda}_{2}{\lambda}_{3}1=\left(1+{e}_{1}\right)\left(1+{e}_{2}\right)\left(1+{e}_{3}\right)1\approx {e}_{1}+{e}_{2}+{e}_{3}& (5.22)\end{array}$$
Now we use the fact that for small $x$, ${\left(1+x\right)}^{3}\approx 1+3x$, so we have:
$$\begin{array}{cc}\frac{\Delta V}{{V}_{0}}={e}_{1}+{e}_{2}+{e}_{3}=3p\left({s}_{11}+2{s}_{12}\right)& (5.23)\end{array}$$
Recognizing that the derivative $dP/dV$ in the definition of ${K}_{b}$ can be written as the limit of $p/\Delta V$ for very small $p$ allows us to obtain the expression we want for ${K}_{b}$:$$\begin{array}{cc}{K}_{b}={lim}_{p\to 0}\frac{p}{\Delta V/{V}_{0}}=\frac{1}{3\left({s}_{11}+2{s}_{12}\right)}& (5.24)\end{array}$$
5.3.5.2 Relationship between the Isotropic Elastic Constants:
Because there are only two independent elastic constants for an isotropic system $E$ and $\nu $ can be expressed in terms of some combination of ${K}_{b}$ and $G$. For $E$ the relevant relationship is as follows.
Note that if ${K}_{b}\gg G$, $E\approx 3G$ and $\nu \approx 0.5$.
6 Other SymmetryRelated Constitutive Relationships
Property

Symbol

Field (row)

Response (column)

tensor properties of rank 0 (scalars)

specific heat

$C$ (1x1)

$\Delta T$

$T\Delta S$

tensor properties of rank 1 (vectors)

pyroelectricity

${p}_{i}^{\prime}$ (3x1)

$\Delta T$

${D}_{i}$ (electric displacement)

tensor properties of rank 2

Thermal expansion

${\alpha}_{i}$ (6x1)

$\Delta T$

${e}_{i}$

Dielectric permittivity

${\kappa}_{ij}$ (3x3)

${E}_{j}$ (elec. field)

${D}_{i}$ (electric displacement)

Electrical conductivity

${\sigma}_{ij}$ (3x3)

${E}_{j}$

${j}_{i}$ (current density)

Thermoelectricity

${\Sigma}_{ij}$ (3x3)

$\frac{\partial T}{\partial {x}_{j}}$ (Temp. gradient)

${E}_{i}$ (electric field)

tensor properties of rank 3

Piezoelectricity

${d}_{ij}$ (3x6)

${\sigma}_{j}$ (stress)

${D}_{i}$ (electric displacement)

Converse piezeoelectricty

${d}_{ij}^{\prime}$ (6x3)

${E}_{j}$ (elec. field)

${e}_{i}$ (strain)

tensor properties of rank 4

Elasticity

${s}_{ij}$ (6x6)

${\sigma}_{j}$

${e}_{i}$

Table 3: Some symmetryrelated properties.
7 Contact Mechanics
7.1 Sign conventions
Sign conventions have a tendency to lead to confusion. This issue is particularly problematics in contact mechanics because compressive loads are considered to be positive, but a compressive stress is negative. Here's a summary of the sign conventions relevant to our treatment of contact mechanics:
 $P$ (force): a positive force is compressive
 $\delta $ (displacement): a positive displacement is compressive
 $\sigma $ (stress): a positive stress is tensile
 $e$ (strain): a positive strain is tensile
In order not to get too hung up in issues related to the sign, we define${P}_{t}$ and${\delta}_{t}$ as the tensile loads and displacements:
$$\begin{array}{cc}\begin{array}{c}{P}_{t}=P\\ {\delta}_{t}=\delta \end{array}& (7.1)\end{array}$$
7.2 Flat Punch Indentation
Consider a flatended cylindrical punch with a radius of
$a$ in contact with another material of thickness,
$h$, as shown schematically in Figure
7.2. The material being indented by a punch rests on a rigid substrate. We are interested in the compressive force,
$P$, that accompanies a compressive displacement,
$\delta $, applied to the indenter.
7.2.1 Flat Punch: Approximate Result for an elastic half space.
For an elastic half space ($h\to \infty $), the strain field under the indenter is nonuniform. The largest strains are confined to a region with characteristic dimensions defined by the punch radius, $a$. We can get a very approximate expression for the relationship between the compressive load, $P$, and the compressive displacement, $\delta $ from the following approximate concepts:
 The average strain in the highly deformed region of the sample must increase linearly with $\delta $. Because strain is dimensionless we need to divide $\delta $ by some length scale in the problem to get a strain. For an elastic half space with $h=\infty $ the only length scale in the problem is the punch radius, $a$. So the strain fields must depend on $\delta /a$. We'll take this one step further and assume that $\delta /a$ is an average in a region of volume $\approx {a}^{3}$ under the punch:
$$\begin{array}{cc}{e}_{avg}=\delta /a& (7.2)\end{array}$$
 The average contact stress, ${\sigma}_{avg}$ under the punch can be quantitatively defined by dividing the compressive load by the stress:$$\begin{array}{cc}{\sigma}_{avg}=P/\pi {a}^{2}& (7.3)\end{array}$$
 An approximate relationship between $P$ and $\delta $ is obtained by assuming that the stress and strain are related through the elastic modulus, i.e. ${\sigma}_{avg}=E{e}_{avg}$. Using the equations above for the compliance, ${C}_{0}$:
$$\begin{array}{cc}{C}_{0}\equiv \frac{\delta}{P}\approx \frac{1}{\pi Ea}& (7.4)\end{array}$$
7.2.2 Flat punch  Detailed Result
In a more general situation both of the contacting materials (the indenter and the substrate) may deform to some extent, so the compliance depends on the properties of both materials. If the materials have Young's moduli of ${E}_{1}$ and ${E}_{2}$ and Poisson's ratios of ${\nu}_{1}$ and ${\nu}_{2}$, then the expression for ${C}_{0}$ is:
$$\begin{array}{cccc}\frac{1}{{E}_{r}}& =& \frac{1{\nu}_{1}^{2}}{{E}_{1}}+\frac{1{\nu}_{2}^{2}}{{E}_{2}}& (7.6)\end{array}$$
Note that for a stiff indenter, $\left({E}_{2}\gg {E}_{1}\right)$ we have ${E}_{r}=\frac{{E}_{1}}{1{\nu}_{1}^{2}}$. This is the plane strain modulus that appears in a variety of situations, which we derive below.
7.2.3 Plane strain modulus.
The plane strain modulus, ${E}_{r}$, describes the response of a material when it cannot contract in one of the directions that is perpendicular to an applied tensile stress. It's easy to derive this by using the compliance matrix for an amorphous material, which must look like this;
In writing the compliance matrix this way, we have used the fact that Young's modulus is $1/{s}_{11}$ and the Poisson's ratio is ${s}_{12}/{s}_{11,}$ so we have ${s}_{11}=1/E$ and ${s}_{12}=\nu /E$. Suppose we apply a stress in the 1 direction, and require that the strain in the 2 direction is 0. This requires that a nonzero stress develop in the 2 direction. If we assume that ${\sigma}_{3}=0$, we have:
If ${e}_{2}$ is constrained to be zero, we have:
$$\begin{array}{cc}{\sigma}_{2}=\nu {\sigma}_{1}& (7.9)\end{array}$$
$$\begin{array}{cc}{e}_{1}=\frac{1}{E}\left({\sigma}_{1}{\nu}^{2}{\sigma}_{1}\right)& (7.10)\end{array}$$
The plane strain modulus relates ${\sigma}_{1}$ to ${e}_{1}$, which for the case assumed above (${e}_{2}=0$) gives:
$$\begin{array}{cc}{E}_{r}=\frac{{\sigma}_{1}}{{e}_{1}}=\frac{E}{1{\nu}^{2}}& (7.11)\end{array}$$
7.3 Flat Punch Detachment and the Energy Release Rate
 $a={a}_{0}$: this is the initial contact condition, where the substrate is in contact with the full surface of the indenter.
 $a={a}_{0}/2$: the contact radius has been reduced to half its initial value.
7.3.1 Energy Release Rate for a Linearly Elastic Material
Specifying the stress field is the same as specifying the stored elastic energy. Fracture occurs when available energy is sufficient to drive a crack forward, or equivalently in our punch problem, to reduce the contact area between the punch and the substrate. To begin we define the following variables:
 $W$= work done on system by external stresses
 ${U}_{el}$= elastically stored energy
 $W{U}_{el}$= energy available to drive crack forward.
where ${A}_{c}$ is the crack area. Fracture occurs when the applied energy release rate exceeds a critical value characteristic of the material, defined as the critical energy release rate, ${G}_{C}$. The fracture condition is therefore:
$$\begin{array}{cc}G={G}_{c}& (7.14)\end{array}$$
The lowest possible value of ${G}_{c}$ is 2$\gamma $ where $\gamma $ is surface energy of the material. That's because the minimal energy to break a material into two pieces is the thermodynamic energy associated with the two surfaces. Some typical values for the surface energy of different materials are listed below (note that 1 mJ/m^{2} = 1 erg/cm^{2} = 1 dyne/cm).
 Polymers: 2050 mJ/m^{2} Van der Waals bonding between molecules
 Water: 72 mJ/m^{2} Hydrogen bonding between molecules
 Metals:$\approx $1000 mJ/m^{2} Metallic bonding
We can derive a simple expression for the energy release rate if we assume that the material has a linear elastic response. Consider, for example, an experiment where we apply a tensile force, ${P}_{t}$, to a sample, resulting in a tensile displacement, ${\delta}_{t}$, as illustrated in Figure 7.4a. If the material has a linear elastic response, the behavior is as illustrated in Figure
7.4. Suppose that the crack area remains constant as the material is loaded to a tensile force
${P}_{t}$. The sample compliance,
$C$, is given by the slope of the displacementforce curve:
$$\begin{array}{cc}C={.\frac{d{\delta}_{t}}{d{P}_{t}}}_{{A}_{c}}& (7.15)\end{array}$$
$$\begin{array}{cc}\delta W=\frac{1}{2}{P}_{t}^{2}\delta C& (7.16)\end{array}$$
Because the sample at the beginning an end of this process is unstrained, we have ${U}_{el}$ =0. We can now take the limit where $\delta {A}_{c}$ becomes very small to replace $\delta W/\delta C$ with the derivative, $dW/dC$ to obtain:
7.3.2 Stable and Unstable Contact
If $dG/da<0,$ a decrease in $A$ gives rise to an increase in $G$, and the contact is unstable, so that the indenter rapidly detaches from the indenter once $a$ starts to decrease.
If $dG/dA>0,$ a decrease in $a$ corresponds to a decrease in $G\mathrm{.}$ In this case the contact is stable, and the load (or displacement) must be increased further to continue to decrease the contact area. Detachment in this case occurs gradually as the load continues to increase.
7.3.3 Application of the Griffith Approach to the Flat Punch Problem
$$\begin{array}{cc}\frac{dC}{dA}=\frac{dC}{d{A}_{c}}& (7.19)\end{array}$$
With
$C={C}_{0}=1/2{E}_{r}a$ (Eq.
7.5) we obtain the following expression for
$G$:
$$\begin{array}{cc}G=\frac{{E}_{r}{\delta}_{t}^{2}}{2\pi a}& (7.22)\end{array}$$
It is useful at this point to make the following general observations:
7.3.4 Detachment: Size Scaling
An interesting aspect of Eq. 7.24 is that the pulloff force scales with
${a}^{3/2}$, whereas the punch cross sectional area scales more strongly with
$a$ (
$A=\pi {a}^{2}$). This behavior has some interesting consequences, which we can obtain by dividing
${P}_{c}$ by the punch cross sectional area to obtain a critical pulloff stress,
${\sigma}_{c}$:
$$\begin{array}{cc}{\sigma}_{c}=\frac{{P}_{c}}{\pi {a}^{2}}={\left(\frac{8{E}_{r}{G}_{c}}{\pi a}\right)}^{1/2}& (7.25)\end{array}$$
Note that the detachment stress increases with decreasing punch size.
Let's put in some typical numbers to see what sort of average stresses we end up with:
 ${E}_{r}\approx 1{0}^{9}$ Pa (typical of glassy polymer)
 ${G}_{c}\approx 0.1\phantom{\rule{6px}{0ex}}J/{m}^{2}$ (twice the surface energy of a typical organic material)
 $a\approx 100$ nm (smallest reasonably possible value)
 ${\sigma}_{c}\approx 50$ MPa
In order for stresses to be obtained, the pillars must be separated so that the stress fields in substrate don't overlap. This decreases the maximum detachment stress from the previous calculation by about a factor of 10, so that the largest stress we could reasonably expect is$\approx $5 MPa. That's still a pretty enormous stress, corresponding to 500 N (49 Kg) over a 1 cm${}^{2}$ area. This is still difficult to achieve, however, because it requires that the pillar array be extremely well aligned with the surface of interest, a requirement that is very difficult to meet in practice. Nevertheless, improvements in the pulloff forces can be realized by structuring the adhesive layer, and this effect is largely responsible for the adhesive behavior of geckos and other creatures with highly structured surfaces.
7.3.5 Thickness Effects
When the thickness,
$h$, of the compliant layer between a rigid cylindrical punch punch and a rigid, flat substrate decreases, the mechanics change in a way that makes it more difficult to pull the indenter out of contact with the compliant layer. For the geometry shown in Figure
7.8 we can write the compliance of the material in the following way:
$$\begin{array}{cc}C=\frac{1}{2{E}_{r}a}{f}_{C}& (7.26)\end{array}$$
$$\begin{array}{cc}{f}_{C}={\left[1+1.33\left(a/h\right)+1.33{\left(a/h\right)}^{3}\right]}^{1}& (7.27)\end{array}$$
The behavior of
${f}_{C}$ as a function of
$a/h$ is plotted in Figure
7.9. A series of geometric correction factors can be derived from this expression for
${f}_{C}\mathrm{.}$ The first of these is a correction factor for the compliance of the energy release rate expression with the tensile load,
${P}_{t}$, as the independent variable. In this case we use Eq.
7.18 to get write the expression for the energy release rate in the following form:
Here ${f}_{Gp}$ accounts for deviations in the compliance derivative due to the confinement effects, in this case determined by the ratio between the actual value of $dC/da$ and the value of this quantity for $a/h=0$:
$$\begin{array}{cc}{f}_{Gp}=\frac{dC}{da}/{.\frac{dC}{da}}_{a/h=0}& (7.29)\end{array}$$
Finally, we can use the the fact that ${P}_{t}={\delta}_{t}/C$ to get a similar expression for $G$ in terms of the tensile displacement, ${\delta}_{t}$:
$$\begin{array}{cc}G=\frac{{E}_{r}{\delta}_{t}^{2}}{2\pi a}{f}_{G\delta}& (7.30)\end{array}$$
$${f}_{G\delta}=\frac{{f}_{Gp}}{{f}_{c}^{2}}$$
The confinement functions
${f}_{Gp}$ and
${f}_{G\delta}$ are both equal to one for
$a/h=0$ and are plotted as a function of
$a/h$ for
$\nu =0.5$ in Figure
7.9. A practical consequence of the decrease in
${f}_{Gp}$ with decreased
$h$ is that a larger tensile force is required in order to remove the cylinder from its contact with the compliant layer. With a small value of
${f}_{Gp}$, a larger tensile load needs to be applied in order for
$G$ to exceed the critical energy release rate,
${G}_{c}$.
7.4 Contact of Paraboloids
7.4.1 NonAdhesive Case
Suppose that the indenter is not flat, but has a parabolic profile that can be described by the following expression:
$$\begin{array}{cc}z=R\left(1\pm \sqrt{1{\left(r/R\right)}^{2}}\right)& (7.33)\end{array}$$
For small $x,$ $\sqrt{1x}\approx 1x/2$, so for $r\ll R$ we have:
We can use Eq. 7.35 to substitute for
$a$ and obtain a relationship between
${P}_{h}$ and
${\delta}_{h}$:
$$\begin{array}{cc}{u}_{z}=\frac{\delta}{\pi}\left\{\left(2{\left(r/a\right)}^{2}\right)arcsin\left(a/r\right)+\left(r/a\right){\left(1{\left(a/r\right)}^{2}\right)}^{1/2}\right\}& (7.38)\end{array}$$
7.4.2 Effects of Adhesion on Contact
The easiest way to understand the effect of adhesion on the contact between a parabolic is to consider a hypothetical situation where we turn off the adhesion and bring the indenter into contact with the surface, resulting in the deformation illustrated in Figure
7.10. Now we we turn on the adhesion, and begin retracting the indenter from the surface, maintaining a fixed projected contact radius
$a$. The situation for the case where we have retracted the tip to the point where the tip apex is level with the undeformed surface (
${\delta}_{t}=0$) is illustrated in Figure
7.11. The applied compressive load required to reach a given contact radius is less than the value of
${P}_{h}$ given by Eq.
7.36 (
$P<{P}_{h}$). Similarly, the compressive displacement required to reach a given contact radius is less than the value given by Eq.
7.35 (
$\delta <{\delta}_{h}$) . These deviations from
${\delta}_{h}$ and
${P}_{h}$ are related by the system compliance, which for this geometry is
${C}_{0}$ as given by Eq.
7.5:
$$\begin{array}{cc}\delta =\frac{{a}^{2}}{3R}+\frac{P}{2{E}_{r}a}& (7.40)\end{array}$$
This expression is the one that needs to be used in order to obtain the reduced modulus in situations where adhesive forces between the indenter and the substrate modify the contact radius. It use requires that the contact radius be measured independently. This is easy to do when the contact area is big enough to visualize directly, but is a very difficulty problem for very small contacts (as in atomic force microscopy) where the contact is too small to visualize optically.
$$\begin{array}{cc}G=\frac{{\left({P}_{t}+{P}_{h}\right)}^{2}}{2}\frac{dC}{dA}=\frac{{\left({P}_{t}+{P}_{h}\right)}^{2}}{8\pi {E}_{r}{a}^{3}}& (7.41)\end{array}$$
This equation can be rearranged to give
${a}^{3}$ as a function of the compressive load,
$P$ (
$P={P}_{t}$), to give an expression that was derived in 1971 by Johnson, Kendall and Roberts [
7] and commonly referred to as the
JKR equation
:
$$\begin{array}{cc}{a}^{3}=\frac{3R}{4{E}_{r}}\left(P+3\pi GR+{\left(6\pi GRP+{\left(3\pi GR\right)}^{2}\right)}^{1/2}\right)& (7.42)\end{array}$$
7.5 Indentation with Berkovich Trips
8 Fracture
 T_{1}: Brittle behavior. This is generally observed at sufficiently low temperatures.
 T_{2}: Ductile behavior (yield before fracture)
 T_{3}: cold drawing (stable neck)
 T_{4}: uniform deformation
Here we are concerned with brittle behavior$\left({T}_{1}\right)$, or in some cases situations where there is a small degree of ductility in the sample (${T}_{2}$).. There are two equivalent approaches for describing the fracture behavior. The first of these is the energy based approach described in the previous section, where an existing crack in a material grows when the applied energy release rate is larger than some critical value. In this section we explore the second approach, where characteristic stress field in the vicinity of a crack exceeds some critical value.
8.1 Fracture Modes
8.2 Stress Concentrations
Note that we recover the behavior described above for a circular whole, where ${a}_{c}={b}_{c}$ and ${\sigma}_{max}/{\sigma}_{0}$=3. We can also write this in terms of the radius of curvature of the ellipse, ${\rho}_{c}$, at the point of maximum stress:
$$\begin{array}{cc}{\sigma}_{max}\propto {\sigma}_{0}\sqrt{{a}_{c}}& (8.4)\end{array}$$
This combination of parameters, with the applied stress multiplied by the square root of the crack length, plays a very important role in fracture mechanics, as we describe in more detail below.
8.3 Stress Intensity Factor
$$\begin{array}{cc}\sigma =\frac{K}{\sqrt{2\pi d}}f\left(\theta \right)& (8.5)\end{array}$$
where $K$ is the stress intensity factor, $d$ is the distance from the crack tip and $f\left(\theta \right)$ is some function of the angle $\theta $ that reduces to 1 for the direction directly in front of a crack ($\theta =0$). Different functional forms exist for $f\left(\theta \right)$ for the different stress components ${\sigma}_{xx}$, ${\sigma}_{yy}$, etc. The detailed stress fields depend on the loading mode (Mode I, II or II, or some combination of these), and the corresponding stress fields are specified by the appropriate value of $K$ (${K}_{I}$ for mode I, ${K}_{II}$ for mode II or ${K}_{III}$ for mode III).
Mode I loading
This compact notation is used to specify the three relevant values of $f\left(\theta \right)\mathrm{.}$ For example, for ${\sigma}_{11}$ we have the following:
$$\begin{array}{cc}{\sigma}_{11}=\left({K}_{I}/\sqrt{2\pi d}\right)cos\left(\theta /2\right)\left(1sin\frac{\theta}{2}sin\frac{3\theta}{2}\right)& (8.7)\end{array}$$
The mode I stress intensity factor for this geometry is given by the applied stress, ${\sigma}_{0}$ and the crack length ${a}_{c}$:
Mode II loading
Mode III loading
While mode III loading is often encountered in practical applications, it is generally avoided in experiments aimed at assessing the fracture behavior of materials, and is not considered further in this text.
8.4 Fracture condition
In the stressbased theory of fracture, the material fails when the stress intensity factor reaches a critical value that depends on the material. For mode I loading, we refer to this critical stress intensity factor as
${K}_{IC}$. Setting
${\sigma}_{0}$ to the fracture stress,
${\sigma}_{f}$, and setting
${K}_{I}$ to
${K}_{IC}$ in Eq.
8.8 gives:
This expression is really only valid for a crack tip with a welldefined radius of curvature, which is often not the case. Models that aim to predict and understand the fracture toughness of materials are all based on understanding the details of the yielding processes very close to the crack tip, and the resulting crack shape. We'll return to this issue later. For now we can summarize the stressbased approach fracture mechanics as follows:
 With the exception of a very small region near the crack tip, all of the strains are elastic.
 There is a very small plastic zone in the vicinity of a crack tip, with a characteristic dimension,$\rho $ that is determined by the details of the way the material plastically deforms.
 Fracture occurs when the stress field defined by ${K}_{I}$ reaches a critical value.
8.5 General relationship between $K$ and $G$
The stress intensity factor and the energy release rate are related to one another through the following expression:
Here ${E}_{r}$ is the reduced modulus that is slightly different for plane stress and plane strain conditions:
$$\begin{array}{cc}\begin{array}{c}{E}_{r}=E\phantom{\rule{6px}{0ex}}\phantom{\rule{6px}{0ex}}\left(Plane\phantom{\rule{6px}{0ex}}stress\phantom{\rule{6px}{0ex}}conditions\right)\\ {E}_{r}=\frac{E}{1{\nu}^{2}}\phantom{\rule{6px}{0ex}}\left(Plane\phantom{\rule{6px}{0ex}}strain\phantom{\rule{6px}{0ex}}conditions\right)\end{array}& (8.15)\end{array}$$
Plane stress conditions generally apply for very thin samples, whereas plane strain conditions apply for thick samples, and also for the axisymmetric punch problems that we have discussed earlier in this text.
$$\begin{array}{cc}{U}_{f}=\frac{1}{2}{\sigma}_{f}{e}_{f}=\frac{1}{2}\frac{{\sigma}_{f}^{2}}{{E}_{r}}& (8.16)\end{array}$$
The stress intensity factor, ${K}_{I}$ at the fracture point is proportional to the stress, ${\sigma}_{f}$, and the strain energy release rate, $G$, at the point of failure is proportional to the total stored elastic energy, ${U}_{f}$. This means that the following proportionality must hold:
$$\begin{array}{cc}G\propto \frac{{K}_{I}^{2}}{{E}_{r}}& (8.17)\end{array}$$
This is consistent with Eq.
8.14, but we need to do a more detailed analysis to get the prefactor exactly right.
8.6 Some Specific Geometries
8.6.1 Double cantilever beam geometry
The double cantilever beam geometry
illustrated in Figure 8.8 is a common test used to measure crack propagation in materials. It is commonly used to measure the adhesion between two materials that have been glued together. It consists of to beams, each with width,
$w$, and thickness,
$t$. The crack length,
${a}_{c}$ in this geometry is the distance between the parts of the beam where the force is applied and the beginning of the region of the sample where the two beams are in contact with one another.
For the double cantilever beam geometry the compliance is given by the following expression:
The crack area, ${A}_{c}$ is obtained by the crack length, ${a}_{c}$ by the width of the sample:
$$\begin{array}{cc}{A}_{c}=w{a}_{c}& (8.19)\end{array}$$
So we have:
$$\begin{array}{cc}G=\frac{{P}_{t}^{2}}{2w}\frac{dC}{d{a}_{c}}=\frac{12{a}_{c}^{2}{P}_{t}^{2}}{E{w}^{2}{t}^{3}}& (8.21)\end{array}$$
At fixed load, $G$ increases as the crack length increases  unstable geometry!
can use We Eq.
8.18 to substitute
${\delta}_{t}$ for
${P}_{t}$ and write the compliance in the following way.
$$\begin{array}{cc}G=\frac{3{\delta}_{t}^{2}{t}^{3}E}{16{a}_{c}^{4}}& (8.22)\end{array}$$
At a fixed displacement, crack will grow until $G={G}_{c}$ and then stop. This is a better way to do the experiment.
8.6.2 Flat Punch Geometry: Thick Compliant Layer
For the flat punch case, the following analytic expression exists for the shape of the normal tress distribution directly under the punch (the plane with
$z={\delta}_{t}$ in Figure
7.3):
Here the average stress,${\sigma}_{avg}$ is defined as follows:
$$\begin{array}{cc}{\sigma}_{avg}\equiv \frac{P}{\pi {a}^{2}}& (8.24)\end{array}$$
Note that ${\sigma}_{zz}$ diverges at the edge of the punch ($r=a$). We know that this must be the case because of the stress concentration that exists at the edge of the punch. To get an expression for stress concentration,${K}_{I}$ at this edge, we first define$d$ as the distance from the punch edge:
$$\begin{array}{cc}d\equiv ar& (8.25)\end{array}$$
$$\begin{array}{cc}\frac{{\sigma}_{zz}}{{\sigma}_{avg}}=\frac{1}{2}{\left[2d/a{\left(d/a\right)}^{2}\right]}^{1/2}& (8.26)\end{array}$$
The stress intensity factor describes the stress field near the contact edge, where $a/d$ is small. We can ignore the term involving the square of $a/d$ to obtain the following expression${\sigma}_{zz}$ that is valid near the contact edge:$$\begin{array}{cc}\frac{{\sigma}_{zz}}{{\sigma}_{avg}}\approx \frac{1}{{2}^{3/2}}{\left(\frac{d}{a}\right)}^{1/2}& (8.27)\end{array}$$
Now we can use the following equation to obtain the following expression for$G$ (assuming${K}_{II}={K}_{III}=0$):
This is the same result that we got above (see Eq. 7.20), so everything checks out okay. Note that the extra factor of two in the relationship between
$G$ and
${K}_{I}$ in Eq.
8.29 comes from the fact the punch is rigid, so it has no stored elastic energy. Because elastic energy is stored only on one side of the interface, the value of
$G$ for crack propagation at the interface with the rigid indenter (Figure
7.3) is half the value of
$G$ for a crack propagating through an elastic material (Figure
8.5, for example).
8.6.3 Flat Punch Geometry: Thin Compliant Layer
Decreasing the thickness of the compliant layer also changes the distribution of normal stresses in contact with the layer. These normal stresses are plotted for different values of
$a/h$ in Figure
8.10.
 Max. stress in center for thin, incompressible layers ($\nu =0.5)$.
 Decrease in edge stress singularity (decrease ${K}_{I}$) for thin layers.
Since for very thin layers failure does not initiate from the edge because the driving force for this 'edge crack' vanishes as $h$ becomes very thin. Instead, failure initiates from small defects within the central part of the contact zone where ${\sigma}_{zz}$ is the highest. The mode I stress intensity factor for a circular, internal crack of radius ${a}_{c}$ is given by the following expression:
$$\begin{array}{cc}G=\frac{{K}_{I}^{2}}{2{E}_{r}}=\frac{2{a}_{c}{\sigma}_{zz}^{2}}{\pi {E}_{r}}& (8.31)\end{array}$$
8.7 Fracture Toughness of Materials
In the Griffith (energybased) model of fracture, material fracture occurs when the applied energy release rate,
$G$, exceeds a threshold value,
${G}_{C}$, which is characteristic of the material. This value is called the critical strain energy release rate, and is a measure of the fracture toughness of the material, just as the critical stress intensity factor is a measure of the fracture toughness. The critical values of
$K$ and
$G$ are related to one another through Eq.
8.14. For mode I fracture,
${K}_{II}={K}_{III}=0$, and this equation reduced to the following:
$$\begin{array}{cc}{G}_{IC}=\frac{{K}_{IC}^{2}}{{E}_{r}}& (8.32)\end{array}$$
Note that we have added the subscript 'I' to $G$ to remind ourselves that this number corresponds to mode I fracture condition.
Values of ${G}_{C}$ are a bit easier to understand conceptually than values of ${K}_{IC}$, since${G}_{c}$ is simply the energy required to break a sample. We can obtain some estimates of ${G}_{C}$ by making some assumptions about where energy goes. Typical values of ${G}_{c}$ are as follows:
 ${G}_{C}=2\gamma \phantom{\rule{6px}{0ex}}\left(\approx 0.1\phantom{\rule{6px}{0ex}}J/{m}^{2}\right)$ if only work during fracture is to break Van der Waals bonds
 ${G}_{C}\approx 12\phantom{\rule{6px}{0ex}}J/{m}^{2}$ if only work during fracture is to break covalent or metallic bonds across interface
 ${G}_{C}\gg 1\phantom{\rule{6px}{0ex}}J/{m}^{2}$ if fracture is accompanied by significant plastic deformation of the sample. For the whole fracture mechanics formulation we are using to be valid, the zone of plastic deformation where this energy dissipation is occurring should be small compared to the overall sample size.
Actual values of
${G}_{C}$ are much larger than these values (see Table
4) because a significant amount of plastic deformation occurs near the crack tip.
We can use numbers from Table 4 to say something about the size of the plastically deformed zone in front of a crack tip. We know that in the elastic region directly in front of a propagating crack, the stress scales as
${K}_{I}/\sqrt{2\pi d}$, where
$d$ is the distance in front of the crack tip. If the maximum stress is equal to the yield stress,
${\sigma}_{y}$, then this the material must be yielded for values of
$r$ that give a stress exceeding
${\sigma}_{y}$. A propagating crack has
${K}_{I}={K}_{IC}$, so if the size of the plastic zone is
${h}_{p}$, we have:
$$\begin{array}{cc}{\sigma}_{y}\approx \frac{{K}_{IC}}{\sqrt{2\pi {h}_{p}}}& (8.33)\end{array}$$
Rearranging gives:
$$\begin{array}{cc}{h}_{p}\approx \frac{{K}_{I}^{2}}{2\pi {\sigma}_{y}^{2}}\approx \frac{{G}_{IC}}{{\sigma}_{y}}\frac{E}{2\pi {\sigma}_{y}}& (8.34)\end{array}$$
This formula is approximate because it neglects the fact that yielding of the material actually changes the stress distribution. The details of what is going on in the plastic zone depend on the materials system of interest. Below we give a case study for what happens for some common amorphous, glassy polymers (noncrystalline polymers deformed at temperatures below their glass transition temperature).
8.8 Case Studies in Fracture
8.8.1 Case Study: Transformation Toughening of Zirconia
First of all, why do we care about a material like zirconia (
$Zr{O}_{2})$? It makes a pretty artificial diamond if we can get it in its cubic crystal form, but it is also an important material in a certain class of fuel cells. Fuel cells are classified by the type of electrolyte that enables ions to be transported between the anode and cathode of the fuel cell. In a solid oxide fuel cell, the electrolyte is typically zirconia, heated to a high temperature typically
$\sim 1000{\phantom{\rule{6px}{0ex}}}^{\u25cb}$C, in order to provide sufficient mobility of the oxygen ions, which move through the electrolyte via a vacancy diffusion mechanism. The electrolyte is part of composite structure, in contact with the anodes and cathodes, and with multiple cells stacked in series to give a structure of the sort shown in Figure
8.11b. Thermal stresses arising as a result of the different thermal expansion coefficients for the different elements can be substantial, and need to be managed appropriately.
The phase diagram for the Yttria/Zirconia system is shown in Figure 8.13. A typical composition has 8% of the Zr atoms replaced with Y. At room temperature this composition consists of the cubic phase in equilibrium with a substantial volume fraction of monoclinic phase particles. In general, however, the tetragonal phase particles formed at high temperature remain in the material as metastable particles. When a crack begins to propagate through the material, the tensile stresses in the vicinity of the crack transform these metastable tetragonal particles to the highervolume, monoclinic particles. The material expansion of associated with this transformation reduces the stress field in front of the crack as illustrated in Figure
8.14, resulting in a substantial toughening of the material.
8.8.2 Tempered Glass
https://www.youtube.com/watch?v=xef4gokRBs
8.8.3 Fiber Composites
8.8.4 Nondestructive Testing of a Fiber Laminate Composite
9 Weibull Analysis of Failure
 b) A crack is nucleated at a defect site at the interface between the film and the substrate at a region where the hydrostatic pressure is maximized.
 c) This crack propagates under the indenter as the material is the indenter is pulled away from the substrate.
 d) Eventually the entire film in has detached from the substrate over the region where it is contact with the indenter. The remainder of the film begins to peel away from the substrate surface.
 e) The film breaks at the edge of the area of contact with the indenter.
$$\begin{array}{cc}ln\left[ln\left(1/{P}_{s}\right)\right]=Mln{p}_{max}Mln{p}_{cav}& (9.2)\end{array}$$
This means that we can obtain the Weibull modulus as the slope of a plot of $ln\left[ln\left(1/{P}_{s}\right)\right]$ vs. ${p}_{max}$. The procedure for obtaining ${P}_{s}$ as a function of ${p}_{max}$ is as follows:
 Start with a data set that includes the measured values of the critical stress at which the sample failed. In our example this would be a list of values of ${p}_{max}$ at the point where the sample failed.
 Organize this list from the lowest value of ${p}_{max}$ to the highest value.
 Use the list to obtain the survival probability, ${P}_{s}$, for each value of ${p}_{max}$. The survival probability is the fraction of samples that did NOT fail for the given value of ${p}_{max}\mathrm{.}$ For example, if our data set has 50 samples in it then the survival probability for the lowest measured value of ${p}_{max}$ is 49/50. The survival probability for the next highest value of ${p}_{max}$ is 48/50, etc. We do this for all of the samples except for the one with the highest value of ${p}_{max}$, since a value of 0 for ${P}_{s}$ would cause problems in the analysis.
In our example the stress is expressed as local maximum in the hydrostatic tension,
${p}_{max}$, and
${p}_{cav}$ is a characteristic value of
${p}_{max}$ for which a substantial fraction of samples have failed. From Eq.
9.1 we see that the survival probability is 1/e=1/2.72=0.37. A Weibull analysis can be applied in a range of situations, including tensile failure of brittle glass rods. In this case the Weibull analysis applies, with the tensile stress
$\sigma $ as the dependent variable, and with a characteristic tensile stress,
${\sigma}_{avg}$ for which the survival probability is 37%:
The point of the Weibull analysis is to obtain an expression that can be sensibly extrapolated to survival probabilities that are very close to 1. With $exp\left(x\right)\approx 1x$ for small $x,$ and with ${P}_{f}=1{P}_{f}$, we have the following expression for failure probability, assuming ${P}_{f}$ is low:
10 Toughening Mechanisms
$$\begin{array}{cc}{G}_{C}=\frac{{K}_{IC}}{{E}_{r}}& (10.1)\end{array}$$
$$\begin{array}{cc}{K}_{I}={\sigma}_{0}\sqrt{\pi a}={K}_{tip}+{K}_{s}& (10.2)\end{array}$$
Here ${K}_{tip}$ describes the stress distribution very close to the crack tip in the usual way. For example, the tensile stress directly in front of a crack tip is given as:
$$\begin{array}{cc}{\sigma}_{zz}\left(d\right)=\frac{{K}_{tip}}{\sqrt{2\pi d}}& (10.3)\end{array}$$
10.1 Hydraulic Fracturing (Fracking)
https://www.youtube.com/watch?v=VY34PQUiwOQ
11 Yield Behavior
11.1 Critical Resolved Shear Stress
$$\begin{array}{cc}{\sigma}_{y}=\frac{{\tau}_{CRSS}}{cos\phi cos\lambda}& (11.2)\end{array}$$
For polycrystalline materials, the situation is more complicated, since all crystal orientations will have a different value of the Schmid factor. In a polycrystalline material, each grain has a different value of the Schmid factor, so the situation is more complicated. However, we can get a reasonable approximation by using an appropriate average value for the Schmid factor instead of the maximum possible value that this value can have. We actually need the quantity, $M$, which is the average of the reciprocal Schmid factor for all of the grains in a material:
$$\begin{array}{cc}M=\u27e8\frac{1}{cos\phi cos\lambda}\u27e9& (11.4)\end{array}$$
Here the average is taken over all possible grain orientations of the material. The tensile yield stress is given by the following expression:
$$\begin{array}{cc}{\sigma}_{y}=M{\tau}_{CRSS}=\frac{M}{2}{\sigma}_{y}^{min}& (11.5)\end{array}$$
11.2 Yield Surfaces
The full stress state of a material is defined by the 3 principal stresses, ${\sigma}_{1}^{p}$, ${\sigma}_{2}^{p}$ and ${\sigma}_{3}^{p}$. The stress state of a material can therefore be specified on a 3dimensionsal space where the values of these 3 principal stresses are plotted on three orthogonal axes. The yield surface is the surface in this space that separates the stress states where yielding will occur from those where it will not occur. Here we describe two of the most common yield surfaces, those defined by the Tresca and Von Mises yield criteria.
11.2.1 Tresca Yield Criterion
$$\begin{array}{cc}\frac{{\left{\sigma}_{i}^{p}{\sigma}_{j}^{p}\right}_{max}}{2}>{\tau}_{crit}& (11.6)\end{array}$$
where the 'max' subscript indicates that we take the principal stress difference (${\sigma}_{1}^{p}{\sigma}_{2}^{p}$, ${\sigma}_{2}^{p}{\sigma}_{3}^{p}$ or ${\sigma}_{1}^{p}{\sigma}_{3}^{p}$) with the largest magnitude. The yield stress, ${\sigma}_{y}$, is typically measured in a uniaxial tensile experiment, where ${\sigma}_{2}^{p}={\sigma}_{3}^{p}=0$, and plastic yielding of the material occurs when ${\sigma}_{1}^{p}>{\sigma}_{y}$. In a uniaxial tensile experiment, the maximum shear stress is half the applied tensile stress, so ${\tau}_{crit}={\sigma}_{y}/2$.
11.2.2 Von Mises Yield Criterion
$$\begin{array}{cc}{\sigma}_{e}=\frac{\sqrt{2}}{2}\sqrt{{\left({\sigma}_{2}^{p}{\sigma}_{1}^{p}\right)}^{2}+{\left({\sigma}_{3}^{p}{\sigma}_{1}^{p}\right)}^{2}+{\left({\sigma}_{3}^{p}{\sigma}_{2}^{p}\right)}^{2}}& (11.7)\end{array}$$
Viewed down the hydrostatic line (${\sigma}_{1}^{p}={\sigma}_{2}^{p}={\sigma}_{3}^{p}$), this is what the criteria look like:
11.3 Localized Deformation
A material that obeys the strain hardening law of Eq. 11.17 will fail when a portion of the sample becomes thinner than the remainder of the sample. The overall behavior of the sample is a balance between the fact that the strain is larger in this region, and can therefore support a larger true stress, but the cross section is larger, so that a larger true stress is needed just to maintain a constant force along the length of the sample. To understand how to think about this we need to consider the relationship between the true stress and the engineering stress.
Consider a sample that is being deformed in uniaxial extension, as illustrated in Figure
11.4. A sample with an undeformed cross sectional area of
${A}_{0}$ and undeformed length of
${\ell}_{0}$ is stretched with a force,
$P$. The engineering tensile stress,
${\sigma}_{eng}$, is obtained by dividing the load by the undeformed cross section, and the true tensile stress,
${\sigma}_{t}$ is obtained by dividing the load by the actual cross sectional area of the deformed sample:
$$\begin{array}{cc}{\sigma}_{eng}=P/{A}_{0}& (11.8)\end{array}$$$$\begin{array}{cc}{\sigma}_{true}=P/A& (11.9)\end{array}$$
In general, the bulk modulus of a material is much larger than its yield stress, so the applied stresses associated with yield phenomena are not large enough to significantly change the volume. As a result, the sample deforms at constant volume, so we have:
$$\begin{array}{cc}A\ell ={A}_{0}{\ell}_{0}& (11.10)\end{array}$$
The relationship between the true stress and the engineering stress is therefore as follows:
11.3.1 ConsidÃ©re Construction
The second possibility is that the increased strain in the necked region leads to substantial strain hardening, so that this region of the sample is able to support the larger true stress in that region. Under the appropriate conditions the cross section of the necked region will stabilize at a value that is determined by the stress/strain relationship for the material. The sample deforms by 'drawing' new material into this necked region, as illustrated on the right side of Figure
11.5.
To understand when stable or unstable necking occur, we begin by recognizing that the onset of neck formation corresponds to a maximum tensile force that the material is able to sustain. In mathematical terms:
$$\begin{array}{cc}\frac{d{\sigma}_{eng}}{d\lambda}=0& (11.13)\end{array}$$
We can rewrite this expression in terms of the true stress by recognizing that
${\sigma}_{eng}={\sigma}_{t}/\lambda $ (see Eq.
11.11), from which we obtain:
$$\begin{array}{cc}\lambda \frac{d{\sigma}_{t}}{d\lambda}{\sigma}_{t}=0& (11.14)\end{array}$$
which we rearrange to the following:
This condition is met when a line drawn fro the origin of a plot of ${\sigma}_{t}$ vs. $\lambda $ is tangent to the curve.
11.3.2 Stable and Unstable Necking
Use of the ConsidÃ©re construction is illustrated in Figure 11.6, where we show curves of the true stress vs. extension ratio for a material that does not form a stable neck (part a) and one that does form a stable neck (part b). In part a it is only possible to draw one line originating from the origin that is tangent to the stressstrain curve (point A in Figure
11.6). A necking instability forms when the true stress reaches this value, resulting in a thinneddown region of the sample. This region continues to thin down until the sample breaks. The maximum engineering stress that the sample sees prior to failure is given by the slope of the tangent line in Figure
11.6a.
In Figure
11.6b, it is possible to draw two lines from the origin that are tangent to the curve, with tangent points labeled as A and B. The tangent at point B represents a maximum in the applied force (the stressstrain curve lies below the tangent line) and the tangent at point B represents a minimum in the applied force (the stressstrain curve lies above the tangent line). At point B the neck stabilizes. Additional material is drawn into the necked region with a characteristic draw ratio that given by the value of
$\lambda $ at point B. The engineering stress at which this drawing occurs is less than the engineering stress required to form the neck in the first place. This means that the load during the drawing process to form the stable neck is lower than the stress required to form the neck in the first place. This phenomenon is generally observed in glassy polymeric materials (
$T<{T}_{g}$) or semicrystalline polymers for which stable neck formation is observed.
11.3.3 Case Study  Power Law Strain Hardening
The following equation is often used to describe the behavior of a material after the yield point:
where:
Limiting behavior:
$n=1$ is perfectly elastic behavior, whereas
$n=0$ corresponds to perfectly plastic behavior. Actual values of
$n$ fall somewhere between these two extremes, and are listed in Table
6 for a variety of metals.
11.3.4 Necking Instability in a PowerLaw Strain Hardening Material
$$\begin{array}{cc}{\sigma}_{t}=K{\left(ln\left(\lambda \right)\right)}^{n}& (11.18)\end{array}$$
Differentiating with respect to $\lambda $ gives:
$$\frac{d{\sigma}_{t}}{d\lambda}=\frac{nK{\left(ln\left(\lambda \right)\right)}^{n1}}{\lambda}=\frac{nK{e}_{t}^{n1}}{\lambda}$$
and then using Eq.
11.15 to equate
$d{\sigma}_{\lambda}/d\lambda $ to
$\sigma /\lambda $ gives the following:
$$\frac{nK{e}_{t}^{n1}}{\lambda}=\frac{\sigma}{\lambda}=\frac{K{e}_{t}^{n}}{\lambda}$$
We can see that this equation is only true when the following condition holds:
$$\begin{array}{cc}n={\epsilon}_{t}& (11.19)\end{array}$$
So a measurement of the strain where the necking instability is observed can be used to determine the strain hardening exponent.
11.4 Deformation Twinning
A useful demonstration of deformation twinning is available from here: http://www.doitpoms.ac.uk/ldplib/acoustic/intro.php
12 Strengthening Mechanisms in Metals
Yield in metals occurs by dislocation motion. This means that if we want to increase the yield stress of a material we have two choices:
 We can make a material with a very low dislocation density, thereby eliminating dislocation as a yield mechanism.
 We can do something to the material to make it more difficult for the dislocations to move.
In almost all cases, we choose option two. This is because mechanisms exist for dislocations to be created during the deformation of a bulk material. More specifically, dislocations multiply when they are pinned in some way, as illustrated by the representation of the operation of a FrankRead dislocation source in Figure 12.1. (See the 3161 text to review FrankRead sources if you don't recall the details.) As a result, even if we manage to create a material with a very low dislocation density, a sufficient number of dislocations will always be created once deformation starts. An exception to this rule is single crystalline whiskers, which are so small that the dislocation density can be reduced essentially to zero. This is because it is energetically favorable for dislocations to move to the free surface of the material, maintaining a very low density of dislocations in the material itself. Many of the details of dislocation structure and the relationship to material deformation was covered in 3162. Our focus on this section is on the factors that control the stress required for dislocations to move through a material. Our discussion is relatively brief, and should be supplemented by the reading the corresponding Wikipedia articles on the different strengthening mechanisms[
15,
18,
19,
22], which are written at an appropriate level of detail for 332.
12.1 Review of some Dislocation Basics
The basics of dislocations are covered in 3161. Here are the most important concepts to remember and use:
 Dislocations correspond to a discontinuity of the displacement field. This continuity is quantified by the Burgers vector, $\overrightarrow{b}$.
 Dislocations move along slip planes that contain both $\overrightarrow{b}$ and the dislocation line.
 When a dislocation exits the sample, portions of the sample on either side of the slip plane are displaced by $\overrightarrow{b}$.
 The energy per length of a dislocation, ${T}_{s}$, is proportional to $G{b}^{2}$, where $b$ is the magnitude of $\overrightarrow{b}$. This energy per unit length has the dimensions of a force, and can be viewed as a 'line tension' acting along the dislocation.
 The resolved shear stress required to bend a dislocation through a radius of curvature of $r$ is ${T}_{s}/rb$.
 Dislocations interact with each other through their stress fields ('like' dislocations repel, 'opposite' dislocations attract.
Strengthening schemes based on reducing the mobility of dislocations can be divided into the following general categories, which we will discuss in more detail below:
 Work hardening
 Grain boundary strengthening
 Solid solution hardening
 Precipitation hardening
12.2 Interactions Between Dislocations
Dislocations interact through the long range strain fields induced by the dislocations. Consider two screw dislocations with Burgers vectors $\overrightarrow{{b}_{1}}$ and ${\overrightarrow{b}}_{2}$ that separated by $d$. The stress at dislocation 2 due to the presence of dislocation 1 is given by the following expression:$$\begin{array}{cc}\overrightarrow{\tau}=\frac{G\overrightarrow{b}}{2\pi d}& (12.1)\end{array}$$
We use a vector notation for the stress here to remind us that the force associated with the shear stress is directed along $\overrightarrow{b}$.
This stress induces a force on the dislocation, given by the following expression:
If ${\overrightarrow{b}}_{1}$ and ${\overrightarrow{b}}_{2}$ are pointing in the same direction, the force is positive and the interaction is repulsive. If they are pointing in the opposite direction, the force is negative (attractive). Because the force scales as $1/d$ it has a very long range.
We can also make a qualitative argument based on the energetics to see what is going on. From Eq. 12.2 we see that
$E\propto {b}^{2}$. If the dislocations have opposite signs, the dislocation come together and the net
$b$ is zero. If they have the same sign, then they form a total Burgers vector with a magnitude of
$2b$, and an energy of
$4{b}^{2}$. So the energy is twice as large as it was originally. The energy as a function of separation must look the plot shown in Figure
12.2. If the energy of each individual dislocation is
${E}_{1}$ for
$d\to \infty $, then the total dislocation energy at
$d=0$ is zero for dislocations of opposite sign and
$4{E}_{1}$ for dislocations of the same sign.
Screw dislocations are easy to think about because the stress field is axially symmetric about the dislocation line, and the stress field is always pure shear. We already know from our previous discussion of stress fields that edge dislocations are more complicated. A simple limiting case involves two edge dislocations on the same slip plane, since within the slip plane we are in a state of pure shear. In this case the discussion from the previous section is still valid, and we get the following expression for the interaction force:
$$\begin{array}{cc}{F}_{s}^{\tau}=\frac{G{\overrightarrow{b}}_{2}\cdot {\overrightarrow{b}}_{1}}{2\pi \left(1\nu \right)d}& (12.3)\end{array}$$
The expression is the same as the expression for the screw dislocations, with the extra factor of $1\nu $.
Dislocations with the same sign repel each other when they are in the same glide plane (see Figure 12.2), but they move toward each other so that they line up on top of one another when they are in different glide planes (see Figure
12.3). The overall situation is summarized in Figure
12.4, which shows the regions of attractive and repulsive interactions for the dislocations. These interactions can be understood in terms of the edge dislocation stress field. The different regions correspond to the different orientations of the shear stress. The important aspect is the direction of the shear force (right or left, in the blue boxes) acting on the top part of the each diagram showing the shear orientation. This can be viewed as the shear force that is acting on the extra half plane on the second dislocation, as a result of the stress field setup up by the first dislocation. This second dislocation moves to the left or the right in response to this force.
As a result of this impurity segregation, dislocations collect 'clouds' of impurity atoms. For a dislocation to move, it must either break away from this atmosphere or move the impurity atoms along with it. The net result in either case is an increase in the critical resolved shear stress, a process referred to as
solid solution strengthening
.
12.3 Work Hardening
(Required supplementary reading:
Work hardening originates from interactions between dislocations. Dislocations impede the motion of other dislocations, so the multiplication of dislocations that occurs as a material is deformed results in a harder material. The effect is illustrated in Figure
12.6, which shows the relationship between the shear strength of a material and the dislocation density. In the absence of dislocations, the strength approaches the maximum theoretical shear strength of
$\approx G/6$, implying a tensile strength of
$\approx G/12$. The introduction of dislocation reduces the strength of the material because the critical resolved shear stress to required to move a dislocation in its slip plan is much less than
$G/6$. As more dislocations become available to introduce a macroscopic plastic strain, the strength of the material goes down as shown in Figure
12.6. Above some critical dislocation density, dislocations interact with each other through their stress fields and the material becomes stronger as the dislocation density increases. The dislocation density at which the strength is minimized is actually quite low. In typical materials the dislocation density is high enough so that the material becomes stronger as it is deformed and the dislocation density increases. In other words, the typical range for actual density is in the increasing part of the strength vs. dislocation density curve, as illustrated in Figure
12.6.
12.4 Grain Boundary Strengthening
(Required supplementary reading:
Grain boundaries act as barriers to dislocation motion. As a result samples with small grain sizes, which have more grain boundaries in a given volume of sample, have higher yield stresses than samples with large grain sizes (as long as the grain size isn't extremely small  in the range of tens of nanometers). The HallPetch equation is commonly used to describe the relationship between the yield strength, ${\sigma}_{y}$ and the grain size, $d$:
Here${\sigma}_{i}$ and ${k}_{y}$ are factors determined by fitting to experimental data. Note that the HallPetch equation is based on experimental observations, and is not an expression derived originally from theoretical considerations.
12.5 Precipitation Hardening
(Supplementary reading:
When a dislocation encounters a second phase particle it can either cut directly through the particle and continue on its way through the matrix phase, or it can bend around the particle, as shown for example in the illustration of the FrankRead source in Figure
12.1. We know from 3161 bending mechanism is going to require the larges shear stress when the particles are small, and close together. If the particles are too close together, however, it may be easier for the dislocation to move through the particle itself. The sheared particle has a slightly larger surface area, and the extra energy required to create this surface gives a contribution to the stress required to move the dislocation through the particle, although this contribution to the stress is generally not a primary contribution to the resolved shear stress.
12.6 Solid Solution Strengthening
(Required supplementary reading:
Solid solution hardening is somewhat related to work hardening in that dislocations are interacting by stress fields originating from defects in the material. In work hardening these defects are other dislocations, whereas for solid solutions the stress fields originate from the presence of either substitutional or interstitial impurities in the material.
13 High Temperature Creep in Crystalline Materials
by 'high' temperature for a metal, we generally mean temperatures above half the absolute melting temperature, which we list below in Figure
13.1.
Different creep regimes during an engineering (constant load) creep rupture test are illustrated in Figure 13.4. At the beginning off the experiment the dislocation density increases with increasing strain, the sample strain hardens, and the strain rate decreases with time. This stage in the experiment is referred to as stage I, or primary creep. Following this there is often a linear regime (stage II in Figure
13.4, referred to as secondary creep) where the strain rate and dislocation density are constant. In this regime dislocations are being created and destroyed at the same rate. The strain rate in this steady state regime is the steady state creep rate,
${\dot{\epsilon}}_{ss}$. Steady state creep data for
$Ti{O}_{2}$ at different times and temperatures are shown in Figure
13.5. The final stage in the creep experiment is the tertiary regime, Regime III, where the disloction density is again increasing, ultimately leading to the formation of voids in the sample and to eventual fracture. In many cases the steady state creep regime corresponds to the majority of the experiment. If this is the case, and the failure is independent of the applied stress, then the creep rupture time will be inversely proportional to the steady state creep rate. We'll focus on the secondary, steadystate creep regime in the following discussion.
The data in Figure 13.3 indicate that the effects of stress and temperature are separable, with the same activation energy of 280 J/mol obtained for each of the different stress levels. The overall steady state creep rate in this case can be expressed in the following way:
For the specific case of power law creep, the stress dependence is given by a simple power law:
where $n$ is an empirically determined exponent, with typical values ranging between 1 and 7.
13.1 Dislocation Glide
The primary deformation mechanism we have talked about so far is dislocation glide, where a dislocation moves laterally in a plane containing the Burgers vector. Dislocation glide can only occur when the resolved shear stress acting on it exceeds some critical value that depends on the curvature of the dislocation. The deformation rate due to dislocation glide is obtained from following relationship between strain rate for dislocation glide, the dislocation density, ${\rho}_{d}$ (the total dislocation length per volume, which has units of 1/m${}^{2}$), the magnitude of the Burgers vector for the dislocations, $b$, and the dislocation glide velocity, ${V}_{g}$:
$$\begin{array}{cc}{.\frac{de}{dt}}_{glide}={\rho}_{d}{V}_{g}b& (13.3)\end{array}$$
$$\begin{array}{cc}{V}_{g}\propto sinh\left(\left(\sigma {\sigma}_{crit}\right)v/2{k}_{B}T\right)exp\left({Q}_{glide}/{k}_{B}T\right)& (13.4)\end{array}$$
There is a temperature dependence of dislocation glide, since there is a finite activation energy for dislocation motion, which we refer to as ${Q}_{glide}$. In most cases however, the stress dependence is much more important than the temperature dependence, and we simply say that the sample deforms plastically for $\sigma >{\sigma}_{crit}$. Some deformation mechanisms are still available below. These require atomic diffusion and have an activation energy that is larger than ${Q}_{glide}$, but become important at high temperatures (greater than about half the absolute melting temperature. These mechanisms are discussed below. The first of these (NabarroHerring creep and Coble creep) involve changes in shape of the overall object that mimic the changes shapes of individual grains, with grain shape changing as atoms diffuse between portions of the grain boundaries that are experiencing different stress states. The third mechanism involves dislocation climb.
13.2 NabarroHerring Creep: Bulk Diffusion Within a Grain
NabarroHerring creep is the mechanism by which vacancies move away from regions of the largest tensile stress in the sample. The mechanism is as illustrated in Figure
13.7. The mechanism is based on the effect that stresses have on the local concentration of vacancies in the system.
In the absence of an applied stress the equilibrium number of concentration of vacancies, ${C}_{v}^{0}$ is given by the vacancy formation energy, ${G}_{v}$:
$$\begin{array}{cc}{C}_{v}^{0}=\frac{1}{\Omega}exp\left(\frac{{G}_{v}}{{k}_{B}T}\right)& (13.5)\end{array}$$
where $\Omega $ is the atomic volume (also equal to the volume of a vacancy). Now suppose we have a grain that is experiencing a tensile stress of magnitude $\sigma $ at the top and bottom, and a compressive force of this same magnitude at the sides. We'll call the vacancy concentrations in the tensile and compressive regions of the sample ${C}_{v}^{t}$ and ${C}_{v}^{c}$, respectively, so we have:
$$\begin{array}{cc}\begin{array}{c}{C}_{v}^{t}={C}_{v}^{0}exp\left(\sigma \Omega /{k}_{B}T\right)\\ {C}_{v}^{c}={C}_{v}^{0}exp\left(\sigma \Omega /{k}_{B}T\right)\end{array}& (13.6)\end{array}$$
We can now use Fick's first law to to calculate the steady state flux of vacancies, ${J}_{v}$, from the tensile region to the compressive region:
We approximate the the vacancy concentration gradient, $d{C}_{v}/dx$, as the difference in vacancy concentrations in the tensile and compressive regions of the grain, divided by the grain size,$d$:
$$\begin{array}{cc}{J}_{v}\approx {D}_{v}\left(\frac{{C}_{v}^{t}{C}_{v}^{c}}{d}\right)& (13.9)\end{array}$$
In order to relate the flux (units of vacancies/m${}^{2}\cdot $s) to the velocity of the grain edge (in m/s), we need to multiply ${J}_{v}$ by the atomic volume, $\Omega $:
Now we get the strain rate by dividing this velocity by $d:$
$$\begin{array}{cc}\frac{de}{dt}=\frac{\Omega {J}_{v}}{d}\approx \Omega {D}_{v}\left(\frac{{C}_{v}^{t}{C}_{v}^{c}}{{d}^{2}}\right)=\frac{2{D}_{v}}{{d}^{2}}sinh\left(\frac{\sigma \Omega}{{k}_{B}T}\right)exp\left(\frac{{G}_{v}}{{k}_{B}T}\right)& (13.11)\end{array}$$
For real systems the stresses are small enough so that $\sigma \Omega \ll {k}_{B}T$ , so we can use the fact that $sinh\left(x\right)\approx x$ for small x to obtain:
Vacancy diffusion is a thermally activated process, so we can express the vacancy diffusion coefficient in the following form:
$$\begin{array}{cc}{D}_{v}={D}_{0}^{v}exp\left(\frac{{Q}_{D}^{v}}{{k}_{B}T}\right)& (13.13)\end{array}$$
where ${Q}_{v}$ is the activation free energy for vacancy motion. We can define a combined quantity, ${Q}_{vol}$ that combines the the free energy of vacancy formation with the activation free energy for vacancy diffusion:
Here ${Q}_{vol}$ is the overall activation energy that accounts for the both both the formation and diffusive motion of vacancies diffusing through the interior of the grain.
13.3 Coble Creep: Grain Boundary Diffusion
$$\begin{array}{cc}\frac{de}{dt}\propto \frac{2{D}_{gb}{d}_{gb}}{{d}^{3}}\frac{\sigma \Omega}{{k}_{B}T}exp\left(\frac{{G}_{v}}{{k}_{B}T}\right)& (13.16)\end{array}$$
The grain boundary diffusion is also a thermally activated process, so we can combine fold the activation free energy for grain boundary diffusion into an overall activation energy, which we define as ${Q}_{gb}$:
$$\begin{array}{cc}{D}_{gb}={D}_{0}^{gb}exp\left(\frac{{Q}_{D}^{gb}}{{k}_{B}T}\right)& (13.17)\end{array}$$
$$\begin{array}{cc}{Q}_{gb}={Q}_{v}+{G}_{v}& (13.18)\end{array}$$
$$\begin{array}{cc}\frac{de}{dt}\propto \frac{2{D}_{0}^{gb}{d}_{gb}}{{d}^{3}}\frac{\sigma \Omega}{{k}_{B}T}exp\left(\frac{{Q}_{gb}}{{k}_{B}T}\right)& (13.19)\end{array}$$
Here ${Q}_{gb}$ is less than ${Q}_{vol}$, so Coble creep will be favored over NabarroHerring creep at low temperatures.
13.4 Weertman model of Creep by Dislocation Climb
In 1957 Johannes Weertman (one of the early and most important members of the Materials Science and Engineering Department at Northwestern University) developed a creep model that conforms to Eqs.
13.1 and
13.2, but with a specific exponent,
$n$, of 4.5.
$$\begin{array}{cc}\frac{de}{dt}=A{\sigma}^{3}sinh\left(B{\sigma}^{1.5}/{k}_{B}T\right)exp\left({Q}_{vol}/{k}_{B}T\right)& (13.20)\end{array}$$
13.5 Deformation mechanism maps
The different deformation mechanisms discussed above all have different dependencies on temperature and stress. We can put place them all on a deformation map as shown in Figure
13.9. The different regions correspond to combinations of stress and temperature where the indicated mechanism gives the highest deformation rate, and therefore dominates the deformation process. Note the following features:
 Coble creep and NabarroHerring (NH) creep both occur in polycrystalline materials, where the overall shape of the object mimics the change in shape of an individual grain. These mechanisms are suppressed by increasing the grain size, and do not occur at all in large, single crystals.
 The boundary between Coble creep and NabarroHerring creep is a vertical line, since the stress dependence for both measurements is the same.
 NabarroHerring creep occurs at higher temperatures than Coble creep, because diffusion through grains (the NabarroHerring mechanism) has a higher activation energy than diffusion at grain boundaries (the Coble mechanism).
 The Coble creep mechanism occurs at grain boundaries, it will be favored for materials with more grain boundaries, i.e., materials with a smaller grain size. As the grain size decreases, the line separating the Coble and NabarroHerring regions will move to the right.
 The line separating the NabarroHerring and Dislocation climb regions is horizontal because these mechanisms have the same temperature dependence, i.e., the activation energy is the same for both cases (see Figure 13.6).
14 Deformation of Polymers
14.0.1 Case Study: Fracture toughness of glassy polymers
Deformation is significant, but ${G}_{Ic}$ is still small compared to other engineering materials.
14.0.1.1 Deformation Mechanisms
Suppose we do a simple stress strain experiment on polystyrene. Polystyrene deforms by one of two different mechanisms:
 Shear bands due to strain softening (decrease in true stress after yield in shear).
 Crazing
 requires net dilation of sample (fracture mechanism for PS and PMMA).
Crazes are load bearing  but they break down to form cracks  failure of specimen.
14.0.1.2 Crazing
Fibrils are cold drawn polymer. Extension ratio remains constant as craze widens
Crack propagation:
 1) new fibrils are created at the craze tip
 2) fibrils break to form a true crack at the crack tip
 Crazing requires a stress field with a tensile hydrostatic component${\sigma}_{1}+{\sigma}_{2}+{\sigma}_{3}>0$ (crazes have voids between fibrils)
 Crazing occurs first for PMMA in uniaxial extension (${\sigma}_{2}=0$)
 ${G}_{Ic}$ is determined by energy required to form a craze ($\approx 1000\phantom{\rule{6px}{0ex}}J/{m}^{2}$)
 Crazing requires strain hardening of fibrils  material must be entangled ($M>{M}_{c}$),${M}_{c}$ typically$\approx $30,000 g/mol.
 In general, shear yielding competes with crazing at the crack tip
Meniscus instability mechanism
(fibril formation at craze tip)
Material near the craze tip is strain softened, and can flow like a fluid between two plates.
http://nervous.com/blog/?p=1556
Competition between Shear Deformation and Crazing
Shear deformation is preferable to crazing for producing high toughness.
Plane stress  shear yielding and crazing criteria (for PMMA)
14.0.1.3 Dugdale model
Assumptions:
1) Tensile stress throughout plastic zone is constant value,${\sigma}_{c}$
$$\begin{array}{cc}{G}_{IC}={\delta}_{c}{\sigma}_{c}& (14.1)\end{array}$$
The Dugdale zone (the craze in our polystyrene example) modifies the stress field so that it doesn't actually diverge to infinity, since infinite stresses are not really possible.
Tensile Behavior:
In Figure 14.15 we compare the tensile behavior of normal polystyrene (PS) and high impact polystyrene (HIPS). The rubber content in the material reduces the modulus but substantially increases the integrated area under the stress/strain curve up to the point of failure, which is a measure of the toughness of the material. We can summarize the differences between PS and HIPS as follows:
 PS is brittle, with $E$ =3GPa and relatively low fracture toughness
 HIPS is ductile, with a lower modulus, ($E=2.1$ GPa) and a much larger energy to fracture.
This area under the stress/strain curve is not a very quantitative measure of toughness because we don't have any information about the flaw size responsible for the eventual failure of the material.
deformation via crazing in vicinity of rubber particles (stress concentrators) throughout sample
Samples with Precrack:
(measurement of K_{IC} or G_{IC})
 Deformation limited to region around crack tip
 Much more deformation for HIPS  higher toughness
Impact Tests
15 Linear Viscoelasticity
Silly Putty is an obvious example of a viscoelastic material, with mechanical properties that depend on the time scale of the measurement. The simplest situation to start with is the small strain region, where the stresses are proportional to the strains, but depend on the time of the experiment. This linear viscoelastic regime is the focus of this Section.
15.1 Time and Frequency Dependent Moduli
As discussed in Section4.5, a variety of elastic constants can be defined for an isotropic material, although only two of these are independent. All of these elastic constants can be time dependent, although in this Section we focus on the specific case of a uniaxial tensile test, used to define the Young's modulus,
$E$, for an elastic material. Viscoelastic materials have mechanical properties that depend on the timescale of the measurement. The easiest way to think about this is to imagine a tensile experiment where a strain of
${e}_{0}$ is instantaneously applied to the sample (see the timedependent strain for the relaxation experiment in Figure
15.1). In a viscoelastic material, the resulting stress will decay with time while we maintain the strain at this fixed value. The time dependence of the resulting tensile stress,
$\sigma $, enables us to define a timedependent relaxation modulus
, $E\left(t\right)$:
$$\begin{array}{cc}E\left(t\right)=\frac{\sigma \left(t\right)}{{e}_{0}}& (15.1)\end{array}$$
$$\begin{array}{cc}e\left(t\right)={e}_{0}sin\left(\omega t\right)& (15.2)\end{array}$$
Note that the strain rate, $\frac{de}{dt}$, is also an oscillatory function, with the same angular frequency, but shifted with respect to the strain by 90${}^{\u25cb}$:
$$\begin{array}{cccc}\frac{de}{dt}& =& {e}_{0}\omega cos\left(\omega t\right)={\gamma}_{0}\omega sin\left(\omega t+\pi /2\right)& (15.3)\end{array}$$
The resulting stress is also an oscillatory function with an angular frequency of $\omega $, and is described by its amplitude and by the phase shift of the relative to the applied strain:
$$\begin{array}{cc}\sigma \left(t\right)={\sigma}_{0}sin\left(\omega t+\phi \right)& (15.4)\end{array}$$
Now we can define a complex modulus with real and imaginary components as follows:
There are a couple different ways to think about the complex modulus,${E}^{*}$. As a complex number we can express it either in terms of its real and imaginary components (${E}^{\prime}$ and${E}^{\prime \prime}$, respectively), or in terms of its magnitude,$\left{E}^{*}\right$ and phase,$\phi $. The magnitude of the complex modulus is simply the stress amplitude normalized by the strain amplitude:
$$\begin{array}{cc}\left{E}^{*}\right={\sigma}_{0}/{e}_{0}& (15.6)\end{array}$$
In order to understand the significance of the real and imaginary components of ${E}^{*}$ we begin with the Euler formula for the exponential of an imaginary number:
$$\begin{array}{cc}tan\phi =\frac{{E}^{\prime \prime}}{{E}^{\prime}}& (15.10)\end{array}$$
The loss tangent gives the ratio of the energy dissipated in one cycle of an oscillation to the maximum stored elastic energy during this cycle.
15.2 Torsional Resonator
The resonant angular frequency of the oscillator, ${\omega}_{n}$ is given by the following expression:
$$\begin{array}{cc}{\omega}_{n}=\sqrt{\frac{{K}_{t}}{I}}& (15.13)\end{array}$$
The solution to Eq. 15.11 for the case where
$T=0$ and we just let the fiber move (by twisting to a certain angle and letting it go, for example) is as follows:
Now all we need to do is to replace the modulus, $G,$ with the complex modulus, ${G}^{*}$, and everything works out fine. Here's what happens:
 The torsional stiffness, ${K}_{t}$ gets transformed to a complex torsional stiffness ${K}_{t}^{*}$
 The resonant frequency, ${\omega}_{n}$ gets transformed to a complex resonant frequency, ${\omega}_{n}^{*}$
 The imaginary part of the complex resonant frequency becomes an exponential damping function
To understand this last point, suppose write${\omega}_{n}^{*}$ in the following way:
$$\begin{array}{cc}{\omega}_{n}^{*}={\omega}_{n}^{\prime}+i{\omega}_{n}^{\prime \prime}& (15.16)\end{array}$$
$$\begin{array}{cc}\theta ={\theta}_{0}Real\left(exp\left(i{\omega}_{n}^{\prime}t\right)\right)exp\left({\omega}_{n}^{\prime \prime}t\right)& (15.17)\end{array}$$
Using the Euler formula (Eq.
15.7) to expand
$exp\left(i{\omega}_{n}^{\prime}t\right)$ and taking the real part gives:
$$\begin{array}{cc}\theta ={\theta}_{0}cos\left({\omega}_{n}^{\prime}t\right)exp\left({\omega}_{n}^{\prime \prime}t\right)& (15.18)\end{array}$$
So the imaginary part of the complex frequency describes the decay of oscillation with time.
15.3 Viscoelastic Models
$$\begin{array}{cc}\sigma =\eta \frac{de}{dt}& (15.19)\end{array}$$
Here ${\rho}_{s}$ is the density of the solid sphere, ${\rho}_{\ell}$ is the density of the surrounding liquid, $R$ is the radius of the sphere and $g$ is the gravitational acceleration (9.8 m/s${}^{2}$). This force is balanced by the viscous force exerted by the water on the sphere as the sphere moves the liquid with a velocity, $v$. For a viscous liquid this force is proportional to the velocity. It is also proportional to the viscosity of the liquid and to the radius of the sphere. The specific relationship is as follows:
Iron has a density of 7.87 g/cm${}^{3}$ (7870 kg/m${}^{3}$) and water has a density of 1000 kg/m${}^{3}$ and a viscosity of 10${}^{3}$ Pasec. From this we get a velocity of 16$\mu $m/s.
15.3.1 Maxwell model
$$\begin{array}{cc}\sigma =\eta \frac{d{e}_{1}}{dt}& (15.23)\end{array}$$
$$\begin{array}{cc}\sigma =E{e}_{2}& (15.24)\end{array}$$
The total strain, $e$ is ${e}_{1}+{e}_{2}$, so the time derivative of the total strain is:
$$\begin{array}{cc}\sigma \left(t\right)=E{e}_{0}exp\left(t/\tau \right)& (15.26)\end{array}$$
We can also obtain the solution for the stress in the case where we apply an oscillatory strain. The math is a bit more complicated, but we end up with the following complex modulus:
$$\begin{array}{cc}{E}^{\prime}\left(\omega \right)=\frac{E{\left(\omega \tau \right)}^{2}}{1+{\left(\omega \tau \right)}^{2}}& (15.29)\end{array}$$
$$\begin{array}{cc}{E}^{\prime \prime}\left(\omega \right)=\frac{E\omega \tau}{1+{\left(\omega \tau \right)}^{2}}& (15.30)\end{array}$$
15.3.2 Standard Linear Solid
$$\begin{array}{cc}E\left(t\right)={E}_{r}+{E}_{1}exp\left(t/{\tau}_{1}\right)& (15.31)\end{array}$$
$$\begin{array}{cc}{E}^{*}\left(\omega \right)={E}_{r}+\frac{{E}_{1}i\omega {\tau}_{1}}{1+i\omega {\tau}_{1}}& (15.32)\end{array}$$
So we see that the standard linear solid describes the relaxation of the modulus from an initial value of ${E}_{1}+{E}_{r}$ at very short times to a relaxed value of ${E}_{r}$ at long times. The standard linear solid is the simplest model for describing the transition of an amorphous, crosslinked polymer from glassy polymer behavior to rubbery behavior with typical values of ${E}_{1}$ and ${E}_{r}$ being $1{0}^{9}$ Pa and $1{0}^{6}$ Pa, respectively.
15.3.3 Generalized Maxwell Model
$$\begin{array}{cc}{E}^{*}\left(\omega \right)={E}_{r}+\sum _{j}\frac{{E}_{j}i\omega {\tau}_{j}}{1+i\omega {\tau}_{j}}& (15.34)\end{array}$$
with${\tau}_{j}={\eta}_{j}/{E}_{j}$.
15.3.4 KelvinVoigt Model
The KelvinVoigt model is the simplest model for the description of a creep experiment, where the stress jumps instantaneously from 0 to $\sigma ={\sigma}_{0}$ and we track the strain as a function of time. We have:
$$\begin{array}{cc}e\left(t\right)=\frac{{\sigma}_{0}}{E}\left(1exp\left(t/\tau \right)\right);\phantom{\rule{6px}{0ex}}\phantom{\rule{6px}{0ex}}\phantom{\rule{6px}{0ex}}\phantom{\rule{6px}{0ex}}\phantom{\rule{6px}{0ex}}\tau =\eta /E& (15.36)\end{array}$$The dynamic modulus for the KelvinVoigt element is given by the following expression:
$$\begin{array}{cc}{E}^{*}=E+i\omega \eta & (15.37)\end{array}$$
So that ${E}^{\prime}$ is simply equal to $E,$ and ${E}^{\prime \prime}$ is equal to $\omega \eta $.
16 Creep Behavior
In a creep experiment we apply a fixed stress to a material and monitor the strain as a function of time at this fixed stress. This strain is generally a function of both the applied stress, $\sigma $ and the time, $t$, since the load was applied.
16.1 Creep in the Linear Viscoelastic Regime
$$\begin{array}{cc}e\left(\sigma ,t\right)={e}_{1}\left(\sigma \right)+{e}_{2}\left(\sigma ,t\right)+{e}_{3}\left(\sigma ,t\right)& (16.2)\end{array}$$
In many cases we are interested in situations where the strain does not increase linearly with the applied stress. The yield point is one obvious example of nonlinear behavior. In ideally plastic system, we only have elastic strains for stresses below the yield stress, but above the yield stress the strains are much larger.
16.2 Nonlinear Creep: Potential Separability of Stress and Time Behaviors
The situation becomes much more complicated in the nonlinear regime, where it is nolonger possible to define a stressindependent creep compliance function. In general the strain in the nonlinear regime is a complex function of both the time and the applied stress. In some cases, however, we can separate the stress dependence from the time dependence and write the stress in the following way:
$$\begin{array}{cc}e\left(\sigma ,\phantom{\rule{6px}{0ex}}t\right)=f\left(\sigma \right)J\left(t\right)& (16.4)\end{array}$$
 Measure $e\left(t\right)$ at different different stresses (${\sigma}_{1}$ and ${\sigma}_{2}$) in Figure 16.3. If the ratio $e\left({\sigma}_{1},t\right)/e\left({\sigma}_{2},t\right)$ is constant for all value of $t$, then the separability into stress dependent and timedependent functions works (at least for these two stress levels). These experiment enable us to obtain the time dependent function $J\left(t\right)$.
 To get the stressdependent function, $f\left(\sigma \right)$, it is sufficient to make a series of measurements at a single, experimentally convenient time.
16.3 Use of empirical, analytic expressions
The procedure outlined in the previous Section only works if the ratio $e\left({\sigma}_{1},t\right)/e\left({\sigma}_{2},t\right)$ is independent of the time. This is not always the case. However, it is often possible to fit the data to relatively simple models. These models are similar to the spring and dashpot models of Section15, but a linear stress response is not necessarily assumed. One example corresponds to a nonlinear version of the linear model shown in Figure
16.1. As with the linear model, the strain components are assumed to consist of an elastic strain,
${e}_{1}$, a recoverable viscoelastic strain,
${e}_{2}$ and a plastic strain,
${e}_{3}$. However, we now use nonlinear elements to describe
${e}_{2}$ and
${e}_{3}$, with these two strain components given by the following expressions:
$$\begin{array}{cc}\begin{array}{c}E={E}_{1}\\ {C}_{1}=1/{\eta}_{2}\\ {C}_{2}={E}_{2}/{\eta}_{2}\\ {C}_{3}=1/{\eta}_{3}\end{array}& (16.6)\end{array}$$
This is just one possible nonlinear model that can be used. An additional nonlinear element is obtained from Eyring rate theory, and is described in the following Section.
16.4.1 Material Deformation as a Thermally Activated Process
16.4.1.1 Eyring Rate Law
We can develop an expression for the strain rate by recognizing that the net strain rate is given by the net frequency of hops in the forward direction. The frequency of hops in the forward and reverse directions, which we refer to as ${f}_{1}$ and ${f}_{2}$, respectively, are given as follows:
$$\begin{array}{cc}{f}_{1}={f}_{0}exp\left(\frac{{Q}^{*}v\sigma /2}{{k}_{B}T}\right)=exp\left(\frac{{Q}^{*}}{{k}_{B}T}\right)exp\left(\frac{v\sigma}{2{k}_{B}T}\right)& (16.7)\end{array}$$
$$\begin{array}{cc}{f}_{2}={f}_{0}exp\left(\frac{{Q}^{*}+v\sigma /2}{{k}_{B}T}\right)=exp\left(\frac{{Q}^{*}}{{k}_{B}T}\right)exp\left(\frac{v\sigma}{2{k}_{B}T}\right)& (16.8)\end{array}$$
The net rate of strain is proportional to ${f}_{1}{f}_{2}$, the net frequency of hops in the forward direction:
$$\begin{array}{cc}\frac{de}{dt}=A\left({f}_{1}{f}_{2}\right)=A{f}_{0}exp\left(\frac{{Q}^{*}}{{k}_{B}T}\right)\left[exp\left(\frac{v\sigma}{2{k}_{B}T}\right)exp\left(\frac{v\sigma}{2{k}_{B}T}\right)\right]& (16.9)\end{array}$$
$$\begin{array}{cc}sinh\left(x\right)=\left({e}^{x}{e}^{x}\right)/2& (16.10)\end{array}$$
we obtain the following expression for the strain rate:
16.4.1.2 Low Stress Regime
In the lowstress regime we can use the approximation $sinh\left(x\right)\approx x$ to get the following expression, valid for $v\sigma \ll {k}_{B}T$.
$$\begin{array}{cc}\eta =\frac{{k}_{B}T}{A{f}_{0}v}exp\left(\frac{{Q}^{*}}{{k}_{B}T}\right)& (16.14)\end{array}$$
So the Eyring theory reduces to and Arrhenius viscosity behavior in the linear, lowstress regime.
16.4.1.3 High Stress Regime
In the highstress regime, we use the fact that $sinh\left(x\right)\approx \frac{exp\left(x\right)}{2}$ for large $x$ to obtain the following expression for $v\sigma \gg {k}_{B}T$:
Equivalently, we can write the following:$$\begin{array}{cc}\frac{de}{dt}=A{f}_{0}exp\left(\frac{{Q}^{*}v\sigma /2}{{k}_{B}T}\right)& (16.16)\end{array}$$
The effective activation energy decreases linearly with increasing stress, giving a very nonlinear response. In practical terms the activation volume is obtained by plotting $ln\left(de/dt\right)$ vs. $\sigma ,$ with the slope being equal to $v/2{k}_{B}T$.
16.4.1.4 Physical significance of $v$
16.4.2 Additional Nonlinear Dashpot Elements
17 Materials Selection: A Case Study
The concept of a space elevator
is illustrated in Figure 17.1. The idea is that we run a cable directly from the earth out to a point in space. If the center of mass is at a geosynchronous orbit, the entire assembly orbits the earth at the same angular velocity at which the earth is rotating. To get into space, we no longer need to use a rocket. We can simply 'climb' up the cable at any velocity that we want. The concept is certainly appealing if we can get it to work. But could the concept ever actually work? That is determined by the availability (or lack thereof) of materials with sufficient strength for the cable. We'll need to start with some analysis to Figure out what sort of properties are needed.
$$\begin{array}{cc}F=m{\omega}^{2}r{G}_{gr}{M}_{e}m/{r}^{2}& (17.1)\end{array}$$
where ${G}_{gr}$ is the gravitational constant and ${M}_{e}$ is the mass of the earth. The angular velocity of the earth is 2$\pi $ radians per day, or in more useful units:
$$\omega =\frac{2\pi}{\left(24\phantom{\rule{6px}{0ex}}hr\right)\left(3600\phantom{\rule{6px}{0ex}}s/hr\right)}=7.3x1{0}^{5}\phantom{\rule{6px}{0ex}}{s}^{1}$$
It is convenient to rewrite the first term in terms of ${g}_{0}$, the gravitational acceleration at $r={r}_{0}$ (at the earth's surface):
$$\begin{array}{cc}{g}_{0}=\frac{{G}_{gr}{M}_{e}}{{r}_{0}^{2}}& (17.2)\end{array}$$
The net force on the mass at $r$ can be written as follows:
The net force is zero for $r={r}_{s}$ (geosynchronous orbit):
$$\begin{array}{cc}{r}_{s}={\left(\frac{{g}_{0}{r}_{0}^{2}}{{\omega}^{2}}\right)}^{1/3}=4.1x1{0}^{7}\phantom{\rule{6px}{0ex}}m\phantom{\rule{6px}{0ex}}(22,000\phantom{\rule{6px}{0ex}}miles)& (17.4)\end{array}$$
$$\begin{array}{cc}F=m{g}_{0}{r}_{0}^{2}\left(\frac{r}{{r}_{s}^{3}}\frac{1}{{r}^{2}}\right)& (17.5)\end{array}$$
$$\begin{array}{cc}{F}_{0}=\rho A{g}_{0}{r}_{0}^{2}{\int}_{{r}_{0}}^{{r}_{\ell}}\left(\frac{r}{{r}_{s}^{3}}\frac{1}{{r}^{2}}\right)dr& (17.6)\end{array}$$
After integration we get:
The maximum tension is at $r={r}_{s}$, and is obtained by integrating contributions to the force from${r}_{s}$ to ${r}_{\ell}$:
We have the following numbers:
 ${r}_{0}=6.4x1{0}^{6}$ m (4,000 miles)
 ${r}_{s}=4.1x1{0}^{7}$ m (22,000 miles)
 ${r}_{\ell}=1.5x1{0}^{8}$ m (90,000 miles)
 ${g}_{0}=9.8$ m/s${}^{2}$
From these numbers we get $\frac{{F}_{max}}{\rho A}=\frac{{\sigma}_{max}}{\rho}=4.8x1{0}^{7}\phantom{\rule{6px}{0ex}}\frac{N/{m}^{2}}{Kg/{m}^{3}}$
All is not lost yet, however, since we really haven't optimized the geometry. The design we considered above has a constant radius, which we would really not want to do. What if we optimize the geometry so that material has the largest cross section at $r={r}_{s}$ (where the load is maximized). We'll consider a design where the actual cross section varies in a way that keeps the stress (tensile force divided by cross sectional area) constant. The analysis is a bit more complicated than we want to bother with here, but we get a simple expression for the maximum cross sectional area, ${A}_{s}$ (at $r={r}_{s}$), to the cross sectional area at the earth's surface (${A}_{0},$ at $r={r}_{0}$) :
$$\begin{array}{cc}\frac{{A}_{s}}{{A}_{0}}=exp\left(\frac{0.77{r}_{0}\rho {g}_{0}}{\sigma}\right)& (17.9)\end{array}$$
18 Review Questions
18.1 Stress and Strain
 How are$E$,$G$,$K$,$\nu $ measured?
 How are these quantities related to one another for an isotropic material?
 What are principal stresses and strains?
 What does the strain ellipsoid lake for a state of pure/simple shear and for hydrostatic tension/compression?
 How are engineering shear strains related to principal extension ratios.
 When can the Mohr circle construction be used, and what is determined from it?
 How do we calculate the resolved shear stress and the resolved normal stress?
18.2 Linear Elasticity
 How do I use the stiffness and compliance matrices?
 What do these matrices look like for materials with different symmetries?
 What are typical moduli for metals, plastic and rubber.
 Under what conditions is the Poisson's ratio very close to 0.5? For what materials is this true?
 What is the relationship between elastic moduli and the sound velocity?
 What's the difference between a shear wave and a longitudinal wave? Which ones can travel in liquids and gases, and why?
18.3 Fracture Mechanics
 What is meant by a mode I, II or III crack? Which mode is most relevant for tensile failure of a homogeneous material, and why?
 What do the stress fields look like in front of a mode I crack?
 How is ${K}_{I}$ related to the crack length and crack tip radius of curvature, and applied tensile load?
 How is $K$ related to $G?$
 What is the significance of ${K}_{IC}$ and ${G}_{c}$?
 How is $G$ related to the compliance of a sample?
 How is $G$ determined in a double cantilever beam test?
 What is the difference between $E$ and the reduced modulus ${E}^{*}$ determined from a contact mechanics experiment?
 How is ${E}^{*}$ determined from an indentation test with flatended cylindrical punch on a thick material?
 What happens to the normal stress distribution under a flat punch as an indented, lowmodulus layer gets thinner and thinner? How is this related to the behavior of the compliance ($C),$ and of $G$ and${K}_{I}$.
 How does adding a rubber to polystyrene increase its toughness?
 How does transformation toughening work in$Zr{O}_{2}$?
 What does the loaddisplacement relationship look like for a rigid, curved indenter indenting a soft material? How does this depend on the materials stiffness and the indenter size?
 What is a Weibull distribution of failure probabilities look like? What is the significance of the Weibull modulus?
18.4 Yield Behavior
 What are the Tresca and Von Mises yield Criteria? How are they used to relate yield criteria for different experimental geometries (simple shear, uniaxial extension, etc.)?
 How is the yield stress obtained from a single crystal if you know the critical and the orientation of the slip system?
 How does this change for a polycrystal?
 How can you tell if a material forms a neck by looking at the true stress as a function of engineering strain or extension ratio? How do you determine the natural draw ratio of the material?
 What are the common strengthening mechanisms for metals, and how do they work?
 How does the strength of a metal depend on its grain size?
 How are the hardness and reduced modulus obtained from a nanoindentation curve with a Berkovich indenter?
18.5 Viscoelasticity
 How is the shear modulus obtained from a torsional resonator? What is the significance of the imaginary part of a complex frequency, and where does this complex frequency come from?
 What is the significance of the magnitude and phase of a complex modulus? What about the real and imaginary components?
 What do ${G}^{\prime}$ and ${G}^{\prime \prime}$ look like for a Maxwell model and for a KelvinVoigt model?
 How are springs and dashpots uses to describe viscoelastic behavior?
 How do Maxwell elements and a standard linear solid behave at high and low frequencies?
 What is a viscosity?
18.6 Creep
 What are the assumptions of the Eyring model?
 What is the significance of the activation volume, and how is it obtained from real experimental data?
 What are the creep mechanisms for dislocation creep, NabarroHerring creep and Coble creep? What determines the activation energy for each of these processes? How do they each depend on the applies stress, temperature and grain size?
 How is a creep deformation map used? How can you use this map to determine which processes have the larges stress dependence or activation energy?
19 Finite Element Analysis
Finite element analysis (FEA) or the finite element method (FEM), is a method for numerical solution of field problems. A field problem is one in which the spatial distribution of one or more unknown variables need to be determined. Thus, we may employ FEA to obtain the spatial distribution of displacements in a structure under load, or the distribution of temperature in an engine block. Mathematically, a field problem is described by a differential equation, or by an integral expression. Finite elements are formulated from such mathematical expressions of the problem.
Individual finite elements can be visualized as small pieces of the structure or domain of interest. Let's assume that the field quantity in question is a function $\phi \left(x,y,z\right)$. In each element $\phi $ is approximated to have only a simple spatial variation, e.g., linear, or quadratic. The actual variation in the region is in most cases more complex, so FEA provides an approximate solution. The elements are connected at points called nodes, and the particular arrangement of elements and nodes is called a mesh. Numerically, the FE mesh leads to a system of algebraic equations (a matrix equation) that is to be solved for the unknown field quantities at the nodes. Once the nodal values of $\phi $ are determined, in combination with the assumed spatial variation of $\phi $ within each element, the approximated field is completely determined within that element. Thus the function $\phi \left(x,y,z\right)$ is approximated element by element, in a piecewise manner. The essence of FEA, then, is the approximation by piecewise interpolation of a field quantity.
There are many advantages of FEA over other numerical methods:
 FEA is applicable to a diverse range of problems: heat transfer, stress analysis, magnetic fields, fluid flow and so on.
 There is no restriction on the geometry. The body or region to be analyzed may have very complex shapes.
 Boundary conditions and loading are not restricted. Any combination of distributed or concentrated forces may be applied to different parts of the body.
 Material properties may be linear or nonlinear, isotropic or anisotropic.
In this introductory class, we will consider only linear, steadystate problems that are displacementbased, that is, where the unknown field quantity is the displacement.
19.1 The direct stiffness method
In this Section we survey the entire computational process of linear steadystate FEA. The FEA process consists of the following steps:
 formulation of element matrices,
 their assembly into a structural matrix,
 application of loads and boundary conditions,
 solution of the structural equations,
 extraction of element stresses and strains.
The procedures are generally applicable regardless of element type, and therefore, we will make use of the simplest of elements  the twonode bar element. For this simple element, we can apply physical arguments to obtain the matrices that represent the element behavior, and for a simple structure consisting of a few bar elements, the computations can be done entirely by hand, allowing the various steps to be carefully examined.
19.1.1 Bar element
The simplest structural finite element is the twonode bar element. A bar, in structural terminology, can resist only axial loads. Consider a uniform elastic bar of length
$L$ and elastic modulus
$E$. Often a bar element is represented as a line, as in Figure
19.1, but the element has cross Section al area
$A$. A node is located at each end, labeled 1 and 2. The axial displacements at nodes are
${u}_{1}$ and
${u}_{2}$, which are the
degrees of freedom of the bar element. The forces acting on the nodes are
${f}_{1}$ and
${f}_{2}$, respectively. The internal axial stress
$\sigma $ can be related to nodal forces
${f}_{1}$ and
${f}_{2}$ by freebody diagrams,
$$\begin{array}{cc}A\sigma +{f}_{1}=0,\phantom{\rule{40px}{0ex}}A\sigma +{f}_{2}=0.& (19.1)\end{array}$$In turn,
$\sigma $ is related to elastic modulus
$E$ and axial strain
$\u03f5$,
$$\begin{array}{cc}\sigma =E\u03f5,\phantom{\rule{40px}{0ex}}\u03f5=\frac{{u}_{2}{u}_{1}}{L}\mathrm{.}& (19.2)\end{array}$$Thus, we obtain
$$\begin{array}{cccc}\frac{AE}{L}\left({u}_{1}{u}_{2}\right)& =& {f}_{1}& (19.3)\\ \frac{AE}{L}\left({u}_{2}{u}_{1}\right)& =& {f}_{2}& (19.4)\end{array}$$or, in matrix form,
$$\begin{array}{cc}\left[k\right]\overrightarrow{d}=\overrightarrow{f}& (19.5)\end{array}$$Here $\left[k\right]$ is the element stiffness matrix:$$\begin{array}{cc}\left[k\right]=\left[\begin{array}{cc}k& k\\ k& k\end{array}\right]\phantom{\rule{20px}{0ex}}where\phantom{\rule{20px}{0ex}}k=\frac{AE}{L},& (19.6)\end{array}$$and $\overrightarrow{d}$ and $\overrightarrow{f}$ are the element nodal displacement and the element nodal force vectors, respectively,$$\begin{array}{cc}\overrightarrow{d}=\left\{\begin{array}{c}{u}_{1}\\ {u}_{2}\end{array}\right\},\phantom{\rule{40px}{0ex}}\overrightarrow{f}=\left\{\begin{array}{c}{f}_{1}\\ {f}_{2}\end{array}\right\}\mathrm{.}& (19.7)\end{array}$$Note than $AE/L$ can be regarded as $k$, the stiffness of a linear spring. A bar and a spring therefore have the same behavior under axial load and are represented by the same stiffness matrix.
19.1.2 Structure of bar elements
Now consider a bar structure that consists of three segments attached endtoend as shown in Figure 19.2(a). Each of these segments has length, crossSectional area and elastic modulus of
${L}^{\left(1\right)}$,
${A}^{\left(1\right)}$,
${E}^{\left(1\right)}$,
${L}^{\left(2\right)}$,
${A}^{\left(2\right)}$,
${E}^{\left(2\right)}$, and
${L}^{\left(3\right)}$,
${A}^{\left(3\right)}$,
${E}^{\left(3\right)}$, respectively. Since abrupt discontinuities in crossSectional area will lead to stress concentrations at the segment junctions, we will let
${A}^{\left(1\right)}={A}^{\left(2\right)}={A}^{\left(3\right)}$. The structure is attached to a rigid wall at its left end, and a force
$\overrightarrow{F}$ is applied to its right end. In a timeindependent analysis, each segment of the structure can be modeled by a twonode bar element, resulting in a
finite element mesh consisting of three elements and four nodes. For a known applied force
$\overrightarrow{F}$, the objective is to determine the displacements in the structure, in particular, the displacement of the nodes
${u}_{1}$,
${u}_{2}$,
${u}_{3}$, and
${u}_{4}$.
The nodal displacement vector $\overrightarrow{D}$ and the nodal force vector $\overrightarrow{F}$ are:$$\begin{array}{cc}\overrightarrow{D}=\left\{\begin{array}{c}{u}_{1}\\ {u}_{2}\\ {u}_{3}\\ {u}_{4}\end{array}\right\},\phantom{\rule{40px}{0ex}}\overrightarrow{F}=\left\{\begin{array}{c}{f}_{1}\\ {f}_{2}\\ {f}_{3}\\ {f}_{4}\end{array}\right\}\mathrm{.}& (19.8)\end{array}$$The stiffness equation for the overall structure can be written as$$\begin{array}{cc}\overrightarrow{F}=\left[K\right]\overrightarrow{D}& (19.9)\end{array}$$where $\left[K\right]$ is the $4\times 4$ global stiffness matrix.
The global stiffness matrix $K$ is obtained through a process of breakdown and assembly. The goal of the breakdown step is to obtain the element stiffness equations. We start by ignoring all loads and supports, and disconnecting the elements in the structure as shown in Figure 19.2(b). The element stiffness relation of each twonode bar element then simply obeys the following equation:
$$\begin{array}{cc}{k}^{\left(e\right)}{d}^{\left(e\right)}={f}^{\left(e\right)},& (19.10)\end{array}$$so that
$$\begin{array}{cccc}\left[\begin{array}{cc}{k}^{\left(1\right)}& {k}^{\left(1\right)}\\ {k}^{\left(1\right)}& {k}^{\left(1\right)}\end{array}\right]\left\{\begin{array}{c}{u}_{1}^{\left(1\right)}\\ {u}_{2}^{\left(1\right)}\end{array}\right\}& =& \left\{\begin{array}{c}{f}_{1}^{\left(1\right)}\\ {f}_{2}^{\left(1\right)}\end{array}\right\},& \\ \left[\begin{array}{cc}{k}^{\left(2\right)}& {k}^{\left(2\right)}\\ {k}^{\left(2\right)}& {k}^{\left(2\right)}\end{array}\right]\left\{\begin{array}{c}{u}_{2}^{\left(2\right)}\\ {u}_{3}^{\left(2\right)}\end{array}\right\}& =& \left\{\begin{array}{c}{f}_{2}^{\left(2\right)}\\ {f}_{3}^{\left(2\right)}\end{array}\right\},& (19.11)\\ \left[\begin{array}{cc}{k}^{\left(3\right)}& {k}^{\left(3\right)}\\ {k}^{\left(3\right)}& {k}^{\left(3\right)}\end{array}\right]\left\{\begin{array}{c}{u}_{3}^{\left(3\right)}\\ {u}_{4}^{\left(3\right)}\end{array}\right\}& =& \left\{\begin{array}{c}{f}_{3}^{\left(3\right)}\\ {f}_{4}^{\left(3\right)}\end{array}\right\}\mathrm{.}& \end{array}$$Note that each element quantity is denoted by a superscript
$\left(e\right)$, where
$e=1,2,3$ is the element number, while the subscripts on the
$u$ and
$f$ indicate the node number. The element stiffnesses are given by the relationship
${k}^{\left(e\right)}={E}^{\left(e\right)}{A}^{\left(e\right)}/{L}^{\left(e\right)}$.
Now let's reconnect the elements and consider the structure as a whole. First, it is clear that for the structure to be physical, the displacements must be compatible. In other words, when elements meet at a node, the displacements of the common node for all the elements that share the node must be the same. Otherwise, there would be gaps in the structure at those nodes. In our structure example, it means ${u}_{2}^{\left(1\right)}={u}_{2}^{\left(2\right)}$, and ${u}_{3}^{\left(2\right)}={u}_{3}^{\left(3\right)}$. Therefore, the element superscripts can be dropped from the nodal displacements in equation (11) without any ambiguity. The element stiffness relations, equation (11) can now be expanded to include all the nodal degrees of freedom, using $0$'s as placeholders:$$\begin{array}{cccc}\left[\begin{array}{cccc}{k}^{\left(1\right)}& {k}^{\left(1\right)}& 0& 0\\ {k}^{\left(1\right)}& {k}^{\left(1\right)}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\left\{\begin{array}{c}{u}_{1}\\ {u}_{2}\\ {u}_{3}\\ {u}_{4}\end{array}\right\}& =& \left\{\begin{array}{c}{f}_{1}^{\left(1\right)}\\ {f}_{2}^{\left(1\right)}\\ 0\\ 0\end{array}\right\},& \\ \left[\begin{array}{cccc}0& 0& 0& 0\\ 0& {k}^{\left(2\right)}& {k}^{\left(2\right)}& 0\\ 0& {k}^{\left(2\right)}& {k}^{\left(2\right)}& 0\\ 0& 0& 0& 0\end{array}\right]\left\{\begin{array}{c}{u}_{1}\\ {u}_{2}\\ {u}_{3}\\ {u}_{4}\end{array}\right\}& =& \left\{\begin{array}{c}0\\ {f}_{2}^{\left(2\right)}\\ {f}_{3}^{\left(2\right)}\\ 0\end{array}\right\},& (19.12)\\ \left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& {k}^{\left(3\right)}& {k}^{\left(3\right)}\\ 0& 0& {k}^{\left(3\right)}& {k}^{\left(3\right)}\end{array}\right]\left\{\begin{array}{c}{u}_{1}\\ {u}_{2}\\ {u}_{3}\\ {u}_{4}\end{array}\right\}& =& \left\{\begin{array}{c}0\\ 0\\ {f}_{3}^{\left(3\right)}\\ {f}_{4}^{\left(3\right)}\end{array}\right\}\mathrm{.}& \end{array}$$Second, the total force acting on a node must be equal to the sum of all the element nodal forces. In our example, therefore${f}_{2}={f}_{2}^{\left(1\right)}+{f}_{2}^{\left(2\right)}$, and ${f}_{3}={f}_{3}^{\left(2\right)}+{f}^{\left(3\right)}$, where the absence of a superscript denoting the element number indicates a global or structure nodal force. The global force vector$F$ is therefore$$\begin{array}{cc}F={f}^{\left(1\right)}+{f}^{\left(2\right)}+{f}^{\left(3\right)}=\left({k}^{\left(1\right)}+{k}^{\left(2\right)}+{k}^{\left(3\right)}\right)D=K\phantom{\rule{6px}{0ex}}D& (19.13)\end{array}$$where the element stiffness matrices and force vectors refer to the expand ed matrices and vectors of equation (12). The global stiffness matrix$K$ is thus$$\begin{array}{cc}K=\left[\begin{array}{cccc}{k}^{\left(1\right)}& {k}^{\left(1\right)}& 0& 0\\ {k}^{\left(1\right)}& {k}^{\left(1\right)}+{k}^{\left(2\right)}& {k}^{\left(2\right)}& 0\\ 0& {k}^{\left(2\right)}& {k}^{\left(2\right)}+{k}^{\left(3\right)}& {k}^{\left(3\right)}\\ 0& 0& {k}^{\left(3\right)}& {k}^{\left(3\right)}\end{array}\right],& (19.14)\end{array}$$and the global stiffness equation is$$\begin{array}{cc}\left[\begin{array}{cccc}{k}^{\left(1\right)}& {k}^{\left(1\right)}& 0& 0\\ {k}^{\left(1\right)}& {k}^{\left(1\right)}+{k}^{\left(2\right)}& {k}^{\left(2\right)}& 0\\ 0& {k}^{\left(2\right)}& {k}^{\left(2\right)}+{k}^{\left(3\right)}& {k}^{\left(3\right)}\\ 0& 0& {k}^{\left(3\right)}& {k}^{\left(3\right)}\end{array}\right]\left\{\begin{array}{c}{u}_{1}\\ {u}_{2}\\ {u}_{3}\\ {u}_{4}\end{array}\right\}=\left\{\begin{array}{c}{f}_{1}\\ {f}_{2}\\ {f}_{3}\\ {f}_{3}\end{array}\right\}\mathrm{.}& (19.15)\end{array}$$This process of merging element stiffness equations to form the global equation is calledassembly.
19.1.3 Boundary conditions
However, the system as stated in equation (15) cannot be solved for the displacements because the stiffness matrix$K$ is singular. The singularity is due to unsuppressed rigid body motions. In other words, the structure is free to float around in space, and a unique solution for the nodal displacements cannot be determined. To eliminate rigid body motions the physical support conditions must be applied asboundary conditions. The support condition that prevents our structure from floating around in space is the attachment of the left end of the structure to a rigid wall. In other words, weknow that${u}_{1}=0$. Therefore,${u}_{1}$ ceases to be a degree of freedom and can be removed entirely from equation (15). Thereduced stiffness equation is then$$\begin{array}{cc}\left[\begin{array}{ccc}{k}^{\left(1\right)}+{k}^{\left(2\right)}& {k}^{\left(2\right)}& 0\\ {k}^{\left(2\right)}& {k}^{\left(2\right)}+{k}^{\left(3\right)}& {k}^{\left(3\right)}\\ 0& {k}^{\left(3\right)}& {k}^{\left(3\right)}\end{array}\right]\left\{\begin{array}{c}{u}_{2}\\ {u}_{3}\\ {u}_{4}\end{array}\right\}=\left\{\begin{array}{c}0\\ 0\\ F\end{array}\right\}& (19.16)\end{array}$$which can be solved to determine the unknown displacements${u}_{2}$,${u}_{3}$ and ${u}_{4}$ for a prescribed value of$F$.
19.1.4 Internal stresses
Once the nodal displacements are known, the average axial stress in each element can be found by ${\sigma}^{\left(e\right)}={E}^{\left(e\right)}{\u03f5}^{\left(e\right)}$. The element axial strain${\u03f5}^{\left(e\right)}={e}^{\left(e\right)}/{L}^{\left(e\right)}$, where${e}^{\left(e\right)}$ is the elongation of the element. For example, the average axial stress of element 2 can be determined as$$\begin{array}{cc}{\sigma}^{\left(2\right)}=\frac{{E}^{\left(2\right)}}{{L}^{\left(2\right)}}\left({u}_{3}{u}_{2}\right)\mathrm{.}& (19.17)\end{array}$$
19.1.5 Reaction forces
by solving equation (16) for the unknown displacements${u}_{2}$,${u}_{3}$ and ${u}_{4}$, the entire displacement vector is determined since${u}_{1}=0$ is prescribed by the boundary conditions. Multiplying the complete displacement solution by $K$ we get$$\begin{array}{cc}\left[\begin{array}{cccc}{k}^{\left(1\right)}& {k}^{\left(1\right)}& 0& 0\\ {k}^{\left(1\right)}& {k}^{\left(1\right)}+{k}^{\left(2\right)}& {k}^{\left(2\right)}& 0\\ 0& {k}^{\left(2\right)}& {k}^{\left(2\right)}+{k}^{\left(3\right)}& 0\\ 0& 0& {k}^{\left(3\right)}& {k}^{\left(3\right)}\end{array}\right]\left\{\begin{array}{c}{u}_{1}\\ {u}_{2}\\ {u}_{3}\\ {u}_{4}\end{array}\right\}=\left\{\begin{array}{c}{k}^{\left(1\right)}{u}_{2}\\ 0\\ 0\\ F\end{array}\right\}\mathrm{.}& (19.18)\end{array}$$Note that we have recovered thereaction force at the support:${f}_{1}={k}^{\left(1\right)}{u}_{2}$. It is easy to check that the complete system is in force equilibrium.
19.1.6 Plane Truss Problem
What if there are more than one degree of freedom per node rather than a single degree of freedom per node? Let's consider the 7bar plane truss of Figure
19.3. Each node of the truss has two degrees of freedom, the
$x$direction displacement
$u$, and the
$y$direction displacement
$v$. The nodal forces also consist of
$x$ and
$y$ components. Then, each displacement
${u}_{i}$ in the nodal displacement vector becomes a vector
${\left[\begin{array}{cc}{u}_{i}& {v}_{i}\end{array}\right]}^{T}$, each nodal force
${f}_{i}$ in the force vector becomes a vector
${\left[\begin{array}{cc}{f}_{ix}& {f}_{iy}\end{array}\right]}^{T}$, and each entry in the stiffness matrix becomes a 2 by 2 matrix. The global displacement vector
$D$ and the global force vector
$F$ for the plane truss of Figure
19.3 are therefore
$$\begin{array}{cc}\begin{array}{c}\overrightarrow{D}={\left[{u}_{1}\phantom{\rule{6px}{0ex}}{v}_{1}\phantom{\rule{6px}{0ex}}{u}_{2}\phantom{\rule{6px}{0ex}}{v}_{2}\phantom{\rule{6px}{0ex}}{u}_{3}\phantom{\rule{6px}{0ex}}{v}_{3}\phantom{\rule{6px}{0ex}}{u}_{4}\phantom{\rule{6px}{0ex}}{v}_{4}\phantom{\rule{6px}{0ex}}{u}_{5}\phantom{\rule{6px}{0ex}}{v}_{5}\right]}^{T}\\ \overrightarrow{F}={\left[{f}_{1x}\phantom{\rule{6px}{0ex}}{f}_{1y}\phantom{\rule{6px}{0ex}}{f}_{2x}\phantom{\rule{6px}{0ex}}{f}_{2y}\phantom{\rule{6px}{0ex}}{f}_{3x}\phantom{\rule{6px}{0ex}}{f}_{3y}\phantom{\rule{6px}{0ex}}{f}_{4x}\phantom{\rule{6px}{0ex}}{f}_{4y}\phantom{\rule{6px}{0ex}}{f}_{5x}\phantom{\rule{6px}{0ex}}{f}_{fy}\right]}^{T}\end{array}& (19.19)\end{array}$$
while the stiffness matrix $K$ that relates the displacements and forces in the relation$K\phantom{\rule{6px}{0ex}}D=F$ is a$10\times 10$ matrix.
19.2 Preliminaries
As we've seen in the previous Section , the basic concept in FEM is the subdivision of the domain (structure) into nonoverlapping elements. The response of each element is expressed in terms of a finite number of degrees of freedom characterized as the value of an unknown function, or functions, at a set of nodal points. The simple elements such as bar or beam elements allow direct formulation of the element matrices from physical arguments as we have seen. However, this is not the norm; the formulation of most elements for structural analysis rely on wellestablished tools of stress analysis, including stressstrain relations, straindisplacement relations, and energy considerations.
19.2.1 Equilibrium equations
Moving beyond the simple freebody diagram approach of Figure 19.1, the more general forces and stresses acting on a differential 2D element
$dx\times dy$ is shown in Figure
19.4. Body forces,
${F}_{x}$ and
${F}_{y}$, acting in the
$x$ and
$y$ directions respectively, are defined as forces per unit volume. The static balance of forces in the
$x$ and
$y$ directions lead to the differential equations for equilibrium (in 2D) as:
$$\begin{array}{cc}\begin{array}{c}{\sigma}_{x,x}+{\tau}_{xy,y}+{F}_{x}=0\\ {\tau}_{xy,y}+{\sigma}_{y,y}+{F}_{y}=0\end{array}& (19.20)\end{array}$$
or, in matrix form,$$\begin{array}{cc}{\partial}^{T}\sigma +F=0,& (19.21)\end{array}$$where$$\begin{array}{cc}\partial =\left[\begin{array}{cc}\frac{\partial}{\partial x}& 0\\ 0& \frac{\partial}{\partial y}\\ \frac{\partial}{\partial y}& \frac{\partial}{\partial x}\end{array}\right],\phantom{\rule{40px}{0ex}}\sigma =\left\{\begin{array}{c}{\sigma}_{x}\\ {\sigma}_{y}\\ {\tau}_{xy}\end{array}\right\},\phantom{\rule{40px}{0ex}}F=\left\{\begin{array}{c}{F}_{x}\\ {F}_{y}\end{array}\right\}\mathrm{.}& (19.22)\end{array}$$
The objective of most structural finite element problems is to find the displacement field that satisfies equation (21) for prescribed applied loads and boundary conditions.
The stresses and strain in an element are related through the constitutive matrix $E$ that contains the elastic constants,$$\begin{array}{cc}\sigma =E\u03f5\mathrm{.}& (19.23)\end{array}$$In twodimensions, the strain vector$\u03f5={\left[\begin{array}{ccc}{\u03f5}_{x}& {\u03f5}_{y}& {\gamma}_{xy}\end{array}\right]}^{T}$, and $$\begin{array}{cc}E=\left[\begin{array}{ccc}{E}_{11}& {E}_{12}& {E}_{13}\\ {E}_{21}& {E}_{22}& {E}_{23}\\ {E}_{31}& {E}_{32}& {E}_{33}\end{array}\right]\mathrm{.}& (19.24)\end{array}$$For example, in isotropic, plane stress conditions,$E$ takes the form$$\begin{array}{cc}E=\frac{E}{1{\nu}^{2}}\left[\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& \frac{1\nu}{2}\end{array}\right],& (19.25)\end{array}$$where$E$ is the elastic modulus and $\nu $ is the Poisson's ratio of the material.
Note that the unknown field quantity are the displacements. The straindisplacement relation,$\u03f5=\partial u$, together with the constitutive equation (23), relates the stresses in the element to the displacements, where$u={\left[\begin{array}{cc}u& v\end{array}\right]}^{T}$, and $u$ and $v$ are the displacements in the$x$ and the$y$ directions, respectively. The equilibrium equation (21) can thus be written in terms of, and solved for, the displacements.
19.2.2 Interpolation and Shape Functions
For the simple linear bar problem of the earlier example, it is clear that once the nodal displacements are determined, the entire displacement field of the structure is known. The displacement field within each element simply varies linearly between the two nodal values. For most problems, however, the relation between the values of the unknown field at the nodes and its variation within the element is not known. In the finite element method, therefore, weassume a simple (e.g., linear or quadratic) variation of the field within each element, by interpolation of the nodal values. To interpolate is to formulate a continuous function that satisfies prescribed conditions at a finite number of points. In FEA, the points are nodes of an element, and the prescribed conditions are the nodal values of a field quantity.
19.2.2.1 Linear Interpolation
Linear interpolation between two points$\left({x}_{1},{u}_{1}\right)$ and $\left({x}_{2},{u}_{2}\right)$. The linear interpolating function$\stackrel{\u02c6}{u}$ is$$\begin{array}{cc}\stackrel{\u02c6}{u}\left(x\right)={a}_{0}+{a}_{1}x,& (19.26)\end{array}$$so that
$$\begin{array}{cc}\begin{array}{c}{u}_{1}={a}_{0}+{a}_{1}{x}_{1}\\ {u}_{2}={a}_{0}+{a}_{1}{x}_{x}\end{array}& (19.27)\end{array}$$
Solving for${a}_{0}$ and ${a}_{1}$ and rearranging, we obtain$$\begin{array}{cc}\stackrel{\u02c6}{u}\left(x\right)={N}_{1}\left(x\right){u}_{1}+{N}_{2}\left(x\right){u}_{2}=Nu,& (19.28)\end{array}$$where$$\begin{array}{cc}N=\left[\begin{array}{cc}{N}_{1}& {N}_{2}\end{array}\right],\phantom{\rule{40px}{0ex}}u=\left\{\begin{array}{c}{u}_{1}\\ {u}_{2}\end{array}\right\},& (19.29)\end{array}$$and $$\begin{array}{cc}{N}_{1}\left(x\right)=\frac{{x}_{2}x}{{x}_{2}{x}_{1}},\phantom{\rule{40px}{0ex}}{N}_{2}\left(x\right)=\frac{x{x}_{1}}{{x}_{2}{x}_{1}},& (19.30)\end{array}$$are the two linear shape functions.
Figure 19.5(a) shows a twonode element with linear interpolation of the unknown field
$u\left(x\right)$ of the two nodal values
${u}_{1}$ and
${u}_{2}$. Note that the shape function
${N}_{1}=1$ at node 1, and
${N}_{1}=0$ at node 2, and
vice versa for
${N}_{2}$
19.2.2.2 Quadratic Interpolation
Quadratic interpolation fits a quadratic function to the three points$\left({x}_{1},{u}_{1}\right)$,$\left({x}_{2},{u}_{2}\right)$, and $\left({x}_{3},{u}_{3}\right)$. The interpolation function is$$\begin{array}{cc}\stackrel{\u02c6}{u}\left(x\right)={N}_{1}\left(x\right){u}_{1}+{N}_{2}\left(x\right){u}_{2}+{N}_{3}\left(x\right){u}_{3}=Nu,& (19.31)\end{array}$$where$$\begin{array}{cc}N=\left[\begin{array}{ccc}{N}_{1}& {N}_{2}& {N}_{3}\end{array}\right],\phantom{\rule{40px}{0ex}}u=\left\{\begin{array}{c}{u}_{1}\\ {u}_{2}\\ {u}_{3}\end{array}\right\},& (19.32)\end{array}$$and the three shape functions are$$\begin{array}{cc}{N}_{1}=\frac{\left({x}_{2}x\right)\left({x}_{3}x\right)}{\left({x}_{2}{x}_{1}\right)\left({x}_{3}{x}_{1}\right)},\phantom{\rule{10px}{0ex}}\phantom{\rule{10px}{0ex}}\phantom{\rule{10px}{0ex}}\phantom{\rule{10px}{0ex}}{N}_{2}=\frac{\left({x}_{1}x\right)\left({x}_{3}x\right)}{\left({x}_{1}{x}_{2}\right)\left({x}_{3}{x}_{2}\right)},\phantom{\rule{10px}{0ex}}\phantom{\rule{10px}{0ex}}\phantom{\rule{10px}{0ex}}\phantom{\rule{10px}{0ex}}{N}_{3}=\frac{\left({x}_{1}x\right)\left({x}_{2}x\right)}{\left({x}_{1}{x}_{3}\right)\left({x}_{2}{x}_{3}\right)}& (19.33)\end{array}$$Figure 19.5(b) shows a threenode element with quadratic interpolation of the unknown field
$u\left(x\right)$ of the three nodal values
${u}_{1}$,
${u}_{2}$ and
${u}_{3}$. Note again that the shape function
${N}_{1}=1$ at node 1 while
${N}_{1}=0$ at nodes 2 and 3, and similar behaviors for
${N}_{2}$ and
${N}_{3}$.
19.2.2.3 Linear Triangle
So far, we have only described interpolation of functions of a single variable,
$x$. In 2D and 3D problems, the field or fields are function of two (
$x,y$) or three (
$x,y,z$) independent variables. These interpolations are extensions of the onedimensional interpolations. A linear triangle is an example of a 2D element, and was the first element devised for plane stress analysis. It has 3 nodes and 6 degrees of freedom, as shown in Figure
19.6. The unknown fields are the displacement in the
$x$direction
$u\left(x,y\right)$, and in the
$y$direction
$v\left(x,y\right)$.
The linear triangle uses linear interpolations for both $u$ and $v$, we will call the interpolating functions$\stackrel{\u02c6}{u}$ and $\stackrel{\u02c6}{v}$:
$$\begin{array}{cc}\begin{array}{c}\stackrel{\u02c6}{u}={N}_{1}{u}_{1}+{N}_{2}{u}_{2}+{N}_{3}{u}_{3}\\ \stackrel{\u02c6}{v}={N}_{1}{v}_{1}+{N}_{2}{v}_{2}+{N}_{3}{v}_{3}\end{array}& (19.35)\end{array}$$
where the shape functions are
$$\begin{array}{cc}\begin{array}{c}{N}_{1}\left(x,y\right)=1\frac{x}{{x}_{2}}+\frac{\left({x}_{3}{x}_{2}\right)y}{{x}_{2}{y}_{3}}\\ {N}_{2}\left(x,y\right)=\frac{x}{{x}_{2}}\frac{{x}_{3}y}{{x}_{2}{y}_{3}}\\ {N}_{3}\left(x,y\right)=\frac{y}{{y}_{3}}\end{array}& (19.36)\end{array}$$
Notice that the same shape functions are used for both$\stackrel{\u02c6}{u}$ and $\stackrel{\u02c6}{v}$.
19.2.2.4 Linear Rectangle (Q4)
The linear rectangle is a 4 node element with 8 degrees of freedom as shown in Figure
19.7. The displacements
$u\left(x,y\right)$ and
$v\left(x,y\right)$ are again interpolated in both
$x$ and
$y$directions with the same linear variation, so we have:
$$\begin{array}{cc}\begin{array}{c}\stackrel{\u02c6}{u}\left(x,y\right)={N}_{1}\left(x,y\right){u}_{1}+{N}_{2}\left(x,y\right){u}_{2}+{N}_{3}\left(x,y\right){u}_{3}+{N}_{4}\left(x,y\right){u}_{4}\\ \stackrel{\u02c6}{v}\left(x,y\right)={N}_{1}\left(x,y\right){v}_{1}+{N}_{2}\left(x,y\right){v}_{2}+{N}_{3}\left(x,y\right){v}_{3}+{N}_{4}\left(x,y\right){v}_{4}\end{array}& (19.37)\end{array}$$
where the shape functions are given by
$$\begin{array}{cc}\begin{array}{cc}{N}_{1}\left(x,y\right)=\frac{\left(ax\right)\left(by\right)}{4ab},\phantom{\rule{6px}{0ex}}\phantom{\rule{6px}{0ex}}& N\left(x,y\right)=\frac{\left(a+x\right)\left(by\right)}{4ab}\\ {N}_{3}\left(x,y\right)=\frac{\left(a+x\right)\left(b+y\right)}{4ab},\phantom{\rule{6px}{0ex}}\phantom{\rule{6px}{0ex}}& N\left(x,y\right)=\frac{\left(ax\right)\left(b+y\right)}{4ab}\end{array}& (19.38)\end{array}$$
19.2.2.5 Number of Nodes and Order of Interpolation
19.2.3 Element Size
Imagine a onedimensional FE structure that has an exact solution$u\left(x\right)$ as shown in Figure 19.9. Assuming that the nodal values are exact (which is not always the case), the finite element solution obtained with linear elements is piecewise linear, and coincides with the exact solution only at the nodes. At all other points, there is an error between the exact and finite element solutions. When the mesh is refined, however, the piecewise linear solution approximates the exact solution more closely, and the error is reduced.
A correct FE formulation therefore converges to the exact solution of the mathematical model as the mesh is infinitely refined. Since an infinitely dense mesh is far from practical, the ideal case would be an acceptable accuracy at a reasonable mesh density. The rate of convergence, which describes how fast the solution converges to the exact one, varies for different element types, and can be obtained by analysis, or by a study of results from a sequence of successively refined meshes. In general, the convergence rate of a linear element is 1 for the displacement, which means that if the element length is halved, the error in the displacement is halved as well. A quadratic or secondorder element converges more rapidly, with a convergence rate of 2 for the displacement, meaning that if the element length is halved, the error in the displacement is quartered.
19.3 Example Problem