# 332: Mechanical Behavior of Materials

## 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:
1. Apply basic concepts of linear elasticity, including multiaxial stress-strain relationships through elastic constants for single and polycrystals.
2. 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).
3. Apply fracture mechanics concepts to determine quantitatively when existing cracks in a material will grow.
4. Describe how composite toughening mechanisms operate in ceramic matrix and polymer matrix composties.
5. Derive simple relationships for the composite stiffness and strength based on those of the constituent phases.
6. Exhibit a quantitative understanding of high temperature deformation in metals and ceramics, based on various creep mechanisms relate to diffusional and dislocation flow (Coble, Nabarro-Herring and Dislocation creep/climb mechanisms).
7. 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.
8. Describe the interplay between surface phenomena (environmental attack) and stresses leading to material embrittlement.
9. 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.
10. Use complex moduli to solve mechanics problems involving an oscillatory stress.
11. Prepare and characterize specimens for measurement of mechanical properties.
12. 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.
13. 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 so-called '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.
Figure 3.1: Available materials throughout history (from ref.[3]).

## 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 $\stackrel{\to }{P}$.
• Unit vectors are indicated with a caret above the symbol, like $\stackrel{ˆ}{n}.$
• 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]& \left(4.1\right)\end{array}$
The stress tensor must be symmetric, with ${\sigma }_{xy}={\sigma }_{yx}$. If this were not the case, the torques on the volume element shown above in Figure 4.1 would not balance, and the material would not be in static equilibrium. As a result the two dimensional stress state is specified by three components of the stress tensor:
• 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{ˆ}{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 $\stackrel{\to }{P}$, $\sigma$ and $\stackrel{\to }{n}$ is as follows:
$\begin{array}{cc}\left[\begin{array}{c}{P}_{1}\\ {P}_{2}\\ {P}_{3}\end{array}\right]=A\left[\begin{array}{ccc}{\sigma }_{11}& {\sigma }_{12}& {\sigma }_{13}\\ {\sigma }_{12}& {\sigma }_{22}& {\sigma }_{23}\\ {\sigma }_{13}& {\sigma }_{23}& {\sigma }_{33}\end{array}\right]\left[\begin{array}{c}{n}_{1}\\ {n}_{2}\\ {n}_{3}\end{array}\right]& \left(4.2\right)\end{array}$
or in more compact matrix notation:
$\begin{array}{cc}\stackrel{\to }{P}=A\left[\sigma \right]\stackrel{ˆ}{n}& \left(4.3\right)\end{array}$
Figure 4.2: 3 dimensional stress tensor
In graphical form the relationship is as shown in Figure 4.2. Like the 2-dimensional stress tensor mentioned above, the 3-dimensional stress tensor must also be symmetric in order for static equilibrium to be achieved. There are therefore 6 independent components of the three-dimensional 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 $\stackrel{\to }{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}& \left(4.4\right)\end{array}$Equivalent expressions can be obtained for the stresses acting on the 2 and 3 faces, by setting $\stackrel{\to }{n}=\left(0,1,0\right)$ and $\stackrel{\to }{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).

#### 4.2.1 Specification of the Transformation Matrix

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}^{○}.$ 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]& \left(4.5\right)\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:
Suppose we want to obtain the stress tensor in the transformed coordinate system obtained from a 45${}^{○}$ counterclockwise rotation around the z axis. The rotation matrix is given by Eq. 4.5, with$\phi =4{5}^{○}.$ The following MATLAB code solves for the full transformed tensor, with ${\sigma }_{ij}$ given by Eq. 4.8 and $\left[{\theta }_{ij}\right]$ given by Eq. 4.5 with $\phi =4{5}^{○}$:
```%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,90-phi,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
end
end
sigp %display the transformed tensor components

```
The Python version of the code looks like this:
```#!/usr/bin/env python3
# -*- coding: utf-8 -*-

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,90-phi,90], [90+phi,phi,90], [90,90,0]]
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
```
Figure 4.4: MATLAB output generated by rotate45.m.
The output generated by the MATLAB version of the code is shown in Figure 4.4, and corresponds to the following result:
$\begin{array}{cc}\left[{\sigma }^{\prime }\right]=\left[\begin{array}{ccc}2.5x1{0}^{6}& -2.5x1{0}^{6}& 0\\ -2.5x1{0}^{6}& 2.5x1{0}^{6}& 0\\ 0& 0& 0\end{array}\right]& \left(4.9\right)\end{array}$
Note the following:
• The normal stresses in the 1 and 2 directions are equal to one another.
• The transformed shear stress in the 1-2 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

A much easier way to do the transformation is to use a little bit of matrix math. The approach we use is described in a very nice web page put together by Bob McGinty[11]: A transformation matrix, ${Q}_{ij}$, is obtained by taking the cosines of all of these angles describing the relationship between the transformed and untransformed coordinate axes:
$\begin{array}{cc}\left[Q\right]=cos\left[\theta \right]& \left(4.10\right)\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]& \left(4.11\right)\end{array}$
The transformed stress is now obtained by the following simple matrix multiplication:
$\begin{array}{cc}\left[{\sigma }^{\prime }\right]=\left[Q\right]\left[\sigma \right]{\left[Q\right]}^{T}& \left(4.12\right)\end{array}$
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)& \left(4.13\right)\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]& \left(4.14\right)\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${}^{○}$ 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,90-phi,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

Figure 4.5: 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/2-phi,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]& \left(4.16\right)\end{array}$
Figure 4.6: MATLAB output generated by rotatesym_matrix.m.

#### 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:
Figure 4.7: MATLAB output generated by mohr_circle.m.
In the Mohr's circle construction the strain shear components of the stress tensor are plotted on the x axis and the shear component of the stress tensor is plotted on the y axis. An example is shown in Figure 4.8, where each of the components of the stress tensor is assumed to be a positive number. In this particular example, the coordinate axes need to be rotated by $2\phi$ to get obtain the principal axes. Whether this rotation is clockwise or counterclockwise depends on the sign convention in the definition of the shear stress. We're not going to worry about it here, but you can refer to the Mohr's Circle Wikipedia article[17] for the details (see the Section on the sign convention).
Figure 4.8: Mohr's circle construction.
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)& \left(4.17\right)\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.
Figure 4.9: Three-dimensional Mohr's Circle.
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 non-zero shear stress. Consider, for example, the following stress state:
$\begin{array}{cc}\left[\sigma \right]=\left[\begin{array}{ccc}3& 0& 2\\ 0& 3& 0\\ 2& 0& 5\end{array}\right]\phantom{\rule{6px}{0ex}}MPa& \left(4.18\right)\end{array}$
In this case there is only one non-zero shear stress (${\sigma }_{13}$), so we can determine the principals stresses in the following manner:
1. 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 non-zero shear stress in our example is ${\sigma }_{13},$ one of the principal stresses is ${\sigma }_{22}$=3 MPa.
2. Draw a Mohr's circle construction using the two normal stresses and the non-zero 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.
3. 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(5-4\right)}^{2}+{\sigma }_{13}^{2}}=2.24\phantom{\rule{6px}{0ex}}MPa$
4. The principal stresses are given by the intersections of the circle with the horizontal axis:
${\sigma }^{p}={\sigma }_{c}±{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.
5. 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.
Figure 4.10: Mohr's circle in the 13 plane for the stress state specified by Eq. 4.18.

#### 4.3.2 Critical Resolved Shear Stress for Uniaxial Tension

As an example of the Mohr circle construction we can consider the calculation of the resolved shear stress on a sample in a state of uniaxial tension. The Mohr's circle representation of the stress state is shown in Figure 4.11. The resolved shear stress, ${\tau }_{rss}$, for sample in uniaxial tension is given by the following expression:
$\begin{array}{cc}{\tau }_{rss}={\sigma }_{1}^{p}cos\phi cos\lambda & \left(4.19\right)\end{array}$
$\begin{array}{cc}{\tau }_{rss}={\sigma }_{1}^{p}cos\phi cos\left(9{0}^{○}-\phi \right)& \left(4.20\right)\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)& \left(4.21\right)\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}& \left(4.22\right)\end{array}$
The untransformed 2-dimensional stress tensor looks like this:
$\begin{array}{cc}\left[\sigma \right]=\left[\begin{array}{cc}{\sigma }_{1}^{p}& 0\\ 0& 0\end{array}\right]& \left(4.23\right)\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]& \left(4.24\right)\end{array}$
Figure 4.11: Mohr's circle construction and calculation of the resolved shear stress for a 2-dimensional sample in uniaxial extension.

#### 4.3.3 Principal Stress Calculation

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${}^{○}$ about the 3 axis. The MATLAB script to do this is very simple and is as follows:
```clear all % clear the workspace of previously defined variables
sig=1e6*[2.5,2.5,0;2.5,2.5,0;0,0,0];
[directions,principalstresses]=eigs(sig);
% the columns in 'directions' correspond to the dot product of the
% principal axes with the orignal coordinate system
% The rotation angles are obtained by calculating the inverse cosines
theta=acosd(directions)
principalstresses
```
Here's the output generated by this script:
Figure 4.12: MATLAB output generated by principal_stress_calc.m.
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)& \left(4.25\right)\end{array}$
$\begin{array}{cc}{I}_{1}={\sigma }_{11}+{\sigma }_{22}+{\sigma }_{33}& \left(4.26\right)\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}& \left(4.27\right)\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}& \left(4.28\right)\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$.

### 4.4 Strain

There are 3 related definitions of the strain:
1. Engineering strain
2. Tensor strain
3. 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)$
Figure 4.13: Location of two points,${P}_{1}$ and${P}_{2}$, before and after an applied deformation.

#### 4.4.1 Small Strains

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}& \left(4.29\right)\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}& \left(4.32\right)\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}& \left(4.33\right)\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}& \left(4.35\right)\end{array}$
$\begin{array}{cc}{e}_{1}^{t}=ln\left({\lambda }_{1}\right)& \left(4.36\right)\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.

### 4.5 Deformation Modes

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

Figure 4.15: Schematic representation of 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]& \left(4.38\right)\end{array}$
$\begin{array}{cc}{e}_{xy}=\frac{u}{d}& \left(4.39\right)\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}\epsilon =\left[\begin{array}{ccc}0& {e}_{xy}/2& 0\\ {e}_{xy}/2& 0& 0\\ 0& 0& 0\end{array}\right]& \left(4.40\right)\end{array}$
We're generally going to use engineering strains and not tensor strains, so we generally don't need to worry about the factor of two. The exception is when we want to use a coordinate transformation to find the principal strains. To do this we use a procedure exactly analogous to the procedure described in Section 4.3, but we need to make sure we are using the tensor strains when we do the calculation.
The shear modulus is simply the ratio of the shear stress to the shear strain.
$\begin{array}{cc}G\phantom{\rule{6px}{0ex}}\left(shear\phantom{\rule{6px}{0ex}}modulus\right)=\frac{{\sigma }_{xy}}{{e}_{xy}}& \left(4.41\right)\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

Mohr's cirlcle for strain looks just like Mohr's circle for stress, provided that we use the appropriate tensor components. That means that we need to plot ${e}_{xy}/2$ on the vertical axis and the normal strains on the vertical axis, as shown in Figure 4.16 (where we have used the common notation for the simple shear geometry, with ${e}_{xy}=\gamma$). One thing that we notice from this plot is that $\gamma$ is simply given by the difference between the two principal strains:
$\begin{array}{cc}\gamma ={e}_{1}^{p}-{e}_{2}^{p}={\lambda }_{1}-{\lambda }_{2}& \left(4.42\right)\end{array}$
I
Figure 4.16: Mohr's circle construction for SMALL strains.

#### 4.5.3 Hydrostatic Compression

$\begin{array}{cc}{K}_{b}=-V\frac{dp}{dV}& \left(4.43\right)\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]& \left(4.44\right)\end{array}$
Figure 4.17: Hydrostatic Compression.

#### 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]& \left(4.45\right)\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 1-2 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}}& \left(4.46\right)\end{array}$
$\begin{array}{cc}{e}_{22}=\frac{\Delta w}{w}& \left(4.47\right)\end{array}$
$\begin{array}{cc}\nu =-{e}_{22}/{e}_{33}& \left(4.49\right)\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 non-zero 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& \left(4.50\right)\end{array}$
Figure 4.19: Longitudinal Compression

### 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 }}& \left(4.51\right)\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& \left(4.52\right)\end{array}$
If the compression is applied slowly enough so that the temperature of the gas can equilibrate, we have:
So we expect that for a gas, the compressive modulus just equal to the pressure. The situation is a bit more complicated for gas, since we need to use the adiabatic modulus, which is about 40% higher than the pressure itself. (For a detailed explanation, see the Wikipedia article on the speed of sound (http://en.wikipedia.org/wiki/Speed_of_sound)[20]. The brief explanation is that for sound propagation, the derivative in Eq. 4.53 needs to be evaluated at constant entropy and not constant temperature, because the sound oscillation is so fast that the heat does not have time to escape). With $\rho$=1.2 kg/m${}^{3}$ and $K=1.4x1{0}^{5}$ Pa, we end up with a sound velocity of 344 m/s.

## 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 1-6, 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:
1. Start with a column vector consisting of the 6 elements of an applied stress, and use the compliance matrix to calculate the strains.
2. 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×6$ matrix to relate two 6-element 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]& \left(5.1\right)\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]& \left(5.2\right)\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]& \left(5.3\right)\end{array}$

### 5.3 Symmetry requirements on the compliance (or stiffness) matrix.

#### 5.3.1 Orthorhombic symmetry

Extruded polymer sheets, like the one shown schematically in Figure 5.1 have orthorhombic symmetry, with different elastic properties in the extrusion, thickness and width directions. These materials have orthorhombic symmetry.
Figure 5.1: Schematic representation of an extruded sheet.
For materials with orthorhombic symmetry, the principal axes of stress and strain are identical, and all compliance components relating a shear strain (e4, e5 or e6) to normal stresses (${\sigma }_{1}$,${\sigma }_{2}$ or${\sigma }_{3}$) or to another shear stress must be zero. The stiffness matrix is as shown in Eq. 5.4 below, and there are 9 independent elastic constants. These 9 elastic constants can be identified as follows:
• ${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.
$\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}& 0& 0& 0\\ {s}_{12}& {s}_{22}& {s}_{23}& 0& 0& 0\\ {s}_{13}& {s}_{23}& {s}_{33}& 0& 0& 0\\ 0& 0& 0& {s}_{44}& 0& 0\\ 0& 0& 0& 0& {s}_{55}& 0\\ 0& 0& 0& 0& 0& {s}_{66}\end{array}\right]\left[\begin{array}{c}{\sigma }_{1}\\ {\sigma }_{2}\\ {\sigma }_{3}\\ {\sigma }_{4}\\ {\sigma }_{5}\\ {\sigma }_{6}\end{array}\right]& \left(5.4\right)\end{array}$

#### 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:
1. Many liquid crystalline polymers (e.g. Kevlar).
2. 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 obtain ${s}_{33}$ and ${s}_{13}$ by performing a tensile test along the fiber axis (the 3 direction) as shown in Figure 5.2. The strain in the 3 direction is given by the fractional change in the length of the fiber after application of the load, and the strains in the 1 and 2 directions are given by the fractional changes in the diameter of the fiber:$\begin{array}{cc}{e}_{3}=\Delta \ell /\ell ;\phantom{\rule{6px}{0ex}}\phantom{\rule{6px}{0ex}}{e}_{1}={e}_{2}=\Delta d/d& \left(5.6\right)\end{array}$
We then obtain ${s}_{33}$ and ${s}_{13}$ from Eq. 5.5, recalling that ${\sigma }_{3}$ is the only non-zero 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}& \left(5.7\right)\end{array}$
Figure 5.2: Fiber extension.
##### 5.3.2.2 Fiber Torsion: Measurement of ${s}_{44}$
Figure 5.3: Fiber torsion.
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 }& \left(5.8\right)\end{array}$
The corresponding shear stress is obtained by multiplying by the shear modulus, $G$ characterizing deformation in the 1-2 and 2-3 planes:
$\begin{array}{cc}{\sigma }_{\theta z}=Gr{\theta }_{0}/\ell & \left(5.9\right)\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}$$\right)$ we know the shear modulus, $G$. This shear modulus is simply the inverse of${s}_{44}$:
$\begin{array}{cc}G=\frac{1}{{s}_{44}}& \left(5.11\right)\end{array}$

#### 5.3.3 Fiber compression in 1-2 plane: determination of ${s}_{11}$ and ${s}_{12}$

The last two elastic constants for a material with fiber symmetry can be determined from an experiment where the fiber is confined between two surfaces and compressed as shown in Figure 5.4. The elastic constants can be determined by measuring the width of the contact region between the fiber and the confining surface. If the confining surfaces are much stiffer than the fiber itself, than this contact width, $2b$, is determined by the elastic deformation of the material in the 1-2 plane. If there is no friction between the fiber and the confining surfaces the fiber is allowed to extend in the 3 direction as it is compressed and the contact width is given by the following expression:
$\begin{array}{cc}{b}^{2}=\frac{2F{d}_{0}{s}_{11}}{\ell }& \left(5.12\right)\end{array}$
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}& \left(5.13\right)\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)& \left(5.14\right)\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)& \left(5.15\right)\end{array}$
Figure 5.4: Transverse deformation of a fiber.
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)& \left(5.16\right)\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}& \left(5.19\right)\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}& \left(5.20\right)\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}& \left(5.22\right)\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)& \left(5.23\right)\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)}& \left(5.24\right)\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.
$\begin{array}{cc}E=\frac{9G}{3+G/{K}_{b}}=2G\left(1+\nu \right)& \left(5.25\right)\end{array}$
We can also equate the two expressions for $G$ in Eq. 5.25 to get the following expression for$\nu$:
$\begin{array}{cc}\nu =\frac{3-2G/{K}_{b}}{6+2G/{K}_{b}}& \left(5.26\right)\end{array}$
Note that if ${K}_{b}\gg G$, $E\approx 3G$ and $\nu \approx 0.5$.

## 6 Other Symmetry-Related Constitutive Relationships

Elasticity is just one of many properties that relate some sort of field (like stress) to a response (like the strain). A range of other properties also exist, with different properties defined by the same sort of symmetry relationships described above. Some of these are listed below in Table 3. A symmetry map that relates the appropriate property coefficients for these different cases is shown in Figure 6.1. Here we show 3 of these symmetry maps corresponding to 3 of 32 symmetry classes for materials. These are discussed in more detail in 361, but the concept is introduced here because it relates to the linear elastic properties.
 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 symmetry-related properties.

## 7 Contact Mechanics

Figure 7.1: Indentation if a soft surface with a rigid, flat-ended cylindrical punch of radius ${a}_{0}$.

### 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}& \left(7.1\right)\end{array}$

### 7.2 Flat Punch Indentation

Consider a flat-ended 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.
Figure 7.2: Flat punch contact geometry. For an elastic half space,$h\to \infty .$

#### 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& \left(7.2\right)\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}& \left(7.3\right)\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}& \left(7.4\right)\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}}& \left(7.6\right)\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 non-zero 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}& \left(7.9\right)\end{array}$
Now we can put this value back into Eq. 7.7 and solve for${e}_{1}$:
$\begin{array}{cc}{e}_{1}=\frac{1}{E}\left({\sigma }_{1}-{\nu }^{2}{\sigma }_{1}\right)& \left(7.10\right)\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}}& \left(7.11\right)\end{array}$

### 7.3 Flat Punch Detachment and the Energy Release Rate

If adhesive forces cause the punch to stick to the substrate, we can use fracture mechanics to understand the force required for detachment to occur. The situation is as shown in Figure 7.3 for a flat-ended cylindrical punch with a radius of ${a}_{0}$. The surface profile of the substrate (assumed in this case to be an elastic half space, i.e., $h=\infty$) is given by the following expression[6]:
• $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.
Figure 7.3: Surface profile under a flat punch, from Eq. 7.12.
The decrease in $a$ from ${a}_{0}$ to ${a}_{0}/2$ is accompanied by a decrease in the stored elastic strain energy, and this strain energy is what drives the decrease in the contact area. While it may not be immediately obvious from Figure 7.3, the detachment problem is actually a fracture mechanics problem. This is because the edge of the contact can be viewed as a crack, which grows as the contact area shrinks. In the following section we describe a generalized energy based approach to for quantifying the driving force for the contact area to decrease before applying this approach to the specific problem of a flat cylindrical punch.

#### 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.
The energy release rate, $G$, describes the amount of energy that is used to move a crack forward by some incremental distance. Formally it is described in the following way:
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}& \left(7.14\right)\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/m2 = 1 erg/cm2 = 1 dyne/cm).
• Polymers: 20-50 mJ/m2 Van der Waals bonding between molecules
• Water: 72 mJ/m2 Hydrogen bonding between molecules
• Metals:$\approx$1000 mJ/m2 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 displacement-force curve:
$\begin{array}{cc}C={.\frac{d{\delta }_{t}}{d{P}_{t}}|}_{{A}_{c}}& \left(7.15\right)\end{array}$
(a)
(b)
Figure 7.4: Derivation of the compliance expression for$G$.
Now suppose that the crack area is increased by an amount $d{A}_{c}$ while the load remains fixed at ${P}_{t}$, i.e. the system moves from point 1 to point two on Figure 7.4b. This increases the compliance by an amount $dC$, resulting in corresponding increase in the displacement of ${P}_{t}\delta C$. If we now unload the sample from point 2 back to the origin, the slope of this unloading curve is given by the enhanced compliance, $C+\delta C$. At the end of this loading cycle be have put energy into the sample equal to the shaded area in Figure 7.4b. This is the total work done on the system by the external stresses ($W$ in Eq. 7.13), and given by the following expression:
$\begin{array}{cc}\delta W=\frac{1}{2}{P}_{t}^{2}\delta C& \left(7.16\right)\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:
$\begin{array}{cc}G={lim}_{\delta {A}_{c}\to 0}\frac{\delta W}{\delta {A}_{c}}=\frac{{P}_{t}^{2}}{2}\frac{dC}{d{A}_{c}}& \left(7.17\right)\end{array}$

#### 7.3.2 Stable and Unstable Contact

Two different behaviors are obtained as two contacting materials are separated, depending on the relationship between $G$ and the contact area $A$. These behaviors are referred to as stable contact ($dG/dA>0\right)$ and unstable contact ($dG/dA<0\right)$. The difference between these behaviors is illustrated in Figure 7.5, and can be understood by considering two surfaces that are initially brought into contact to establish a contact area, ${A}_{0}$. We then increase the tensile load, ${P}_{t}$, thereby increasing the applied energy release rate. As the tensile load and the corresponding tensile displacement are increased, $G$ increases until it reaches the critical value, ${G}_{c}$. The tensile load at this point is defined as the critical load ${P}_{c}$, and is the load at which $A$ begins to decrease. We fix the tensile load at ${P}_{c}$ and observe one of two possible behaviors:
Figure 7.5: Illustration of stable contact, where ${P}_{t}$ must increase continuously in order for the contact area to continue to decrease, and unstable contact, where the contact area reduces rapidly to zero as soon as a critical tensile load is attained.
##### Unstable Detachment :
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.
##### Stable Detachment :
If $dG/dA>0,$ a decrease in $a$ corresponds to a decrease in $G.$ 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

The edge of the contact is a crack, which advances as $a$ decreases. We can use Eq. 7.17 for the energy release rate to obtain the following:
$\begin{array}{cc}G=-\frac{{P}_{t}^{2}}{2}\frac{dC}{dA}=-\frac{{P}_{t}^{2}}{4\pi a}\frac{dC}{da}& \left(7.18\right)\end{array}$
where we have assumed that the contact area remains circular, with $A=\pi {a}^{2}$. Not that $A$ in this expression is the contact area between the indenter and the substrate, and NOT the crack area. The negative sign in Eq. 7.18 emerges from the fact an decrease in contact area corresponds to an equivalent increase in the crack area, so we have:
$\begin{array}{cc}\frac{dC}{dA}=-\frac{dC}{d{A}_{c}}& \left(7.19\right)\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{{P}_{t}^{2}}{8\pi {E}_{r}{a}^{3}}& \left(7.20\right)\end{array}$
$\begin{array}{cc}G=-\frac{{\delta }_{t}^{2}}{2{C}^{2}}\frac{dC}{dA}=-\frac{{\delta }_{t}^{2}}{4\pi a{C}^{2}}\frac{dC}{dA}& \left(7.21\right)\end{array}$
If the compliance is the value for an elastic half space (Eq. 7.5), then we obtain the following expression for the energy release rate in terms of the displacement:
$\begin{array}{cc}G=\frac{{E}_{r}{\delta }_{t}^{2}}{2\pi a}& \left(7.22\right)\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 pull-off 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 pull-off 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}& \left(7.25\right)\end{array}$
Note that the detachment stress increases with decreasing punch size.
Figure 7.6: Electron micrograph of Gecko setae.
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 pull-off 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.
Figure 7.7: Schematic representation of an array of pillars in contact with a flat surface.

#### 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}& \left(7.26\right)\end{array}$
Figure 7.8: A thin, compliant layer being indented with a rigid, cylindrical punch.
$\begin{array}{cc}{f}_{C}={\left[1+1.33\left(a/h\right)+1.33{\left(a/h\right)}^{3}\right]}^{-1}& \left(7.27\right)\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}.$ 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:
$\begin{array}{cc}G=-\frac{{P}_{t}^{2}}{4\pi a}\frac{dC}{da}=\frac{{P}_{t}^{2}}{8\pi {E}_{r}{a}^{3}}{f}_{Gp}& \left(7.28\right)\end{array}$
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}& \left(7.29\right)\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 }& \left(7.30\right)\end{array}$
In this case ${f}_{G\delta }$ includes the dependence on $a/h$ of both the compliance and it's derivative with respect to $a$. This dependence is evident from Eq. 7.21, where $G\left({\delta }_{t}\right)$ is seen to be proportional to $dC/da$ and is inversely proportional to ${C}^{2}$. The $a/h$ dependence of $dC/da$ is accounted for by ${f}_{Gp}$, and the $a/h$ dependence of $C$ is accounted for by ${f}_{C}$, so we obtain the following for ${f}_{G\delta }$:
${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

Suppose that the indenter is not flat, but has a parabolic profile that can be described by the following expression:
$\begin{array}{cc}z={A}_{p}{r}^{2}& \left(7.31\right)\end{array}$
$\begin{array}{cc}{r}^{2}+{\left(z-R\right)}^{2}={R}^{2}& \left(7.32\right)\end{array}$
Solving Eq. 7.32 for$z$ gives:
$\begin{array}{cc}z=R\left(1±\sqrt{1-{\left(r/R\right)}^{2}}\right)& \left(7.33\right)\end{array}$
For small $x,$ $\sqrt{1-x}\approx 1-x/2$, so for $r\ll R$ we have:
$\begin{array}{cc}z=\frac{{r}^{2}}{2R}& \left(7.34\right)\end{array}$where we have taken the solution with the smaller value of $z$, corresponding to the bottom of the sphere. From a comparison of Eqs.7.32 and 7.34, we see the paraboloid is a good approximation for the shape of a sphere, with the sphere radius given by $1/2{A}_{p}$. For this reason we use $R$ instead of $A$ to characterize the parabolic shape, since the results can be applied to contact of spheres, provided that the the contact dimensions are much smaller than $R$. Generally everything works well as long as $r/R<4.$

Figure 7.10: Non-adhesive contact of a rigid, parabolic indenter into an elastic material.
The compressive a rigid parabolic indenter into the surface of the material (${\delta }_{h}$ in Figure 7.10) is given by the following expression:
$\begin{array}{cc}{\delta }_{h}={a}^{2}/R& \left(7.35\right)\end{array}$Note that this is a completely geometric relationship that does not depend on the modulus of the material that is being indented. The compressive force required to establish a contact radius of $a$ is referred to as ${P}_{h}$, and is given by the following expression:
We can use Eq. 7.35 to substitute for $a$ and obtain a relationship between ${P}_{h}$ and${\delta }_{h}$:
$\begin{array}{cc}{P}_{h}=\frac{4{E}_{r}}{3}{R}^{1/2}{\delta }_{h}^{3/2}& \left(7.37\right)\end{array}$
$\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\}& \left(7.38\right)\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}\frac{\delta -{\delta }_{h}}{P-{P}_{h}}={C}_{0}=\frac{1}{2{E}_{r}a}& \left(7.39\right)\end{array}$
Combination of Eqs. 7.5 and 7.39 gives the following relationship between $\delta$, $P$ and ${E}_{r}$:
$\begin{array}{cc}\delta =\frac{{a}^{2}}{3R}+\frac{P}{2{E}_{r}a}& \left(7.40\right)\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.
(b)
Figure 7.11: An example of the surface profile for adhesive contact for the case where ${\delta }_{t}=0$.
Once we know the reduced modulus of the system, we can obtain the energy release rate. The expression for the energy release rate for curved object in contact with surface in a way that is very similar to what we did for the flat punch in Section7.3. The only difference is that in the absence of adhesion we need to apply a compressive load, ${P}_{h}$ (given by Eq. 7.36):
$\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}}& \left(7.41\right)\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)& \left(7.42\right)\end{array}$

### 7.5 Indentation with Berkovich Trips

Figure 7.12: Geometry of a Berkovich tip commonly used in indentation experiments. The angle, $a$, is $65.3{5}^{○}$ for a standard Berkovich tip.
The hardness, $H$, of a material is given by the ratio of the load to the projected contact area of the non-recoverable indent made in the material by the indenter. In our case we obtain the hardness from the maximum load,${P}_{max}$ (illustrated in Figure 7.13), and from the corresponding projected area, $A$, of the hardness impression: $\begin{array}{cc}H=\frac{{P}_{max}}{A}& \left(7.43\right)\end{array}$
The projected area is related to the contact depth, ${\delta }_{c}$, by a relationship that depends on the shape of the indenter[9]. For a Berkovich tip the appropriate relationship is:
$\begin{array}{cc}A=24.5{\delta }_{c}^{2}& \left(7.44\right)\end{array}$
where ${\delta }_{max}$ is the maximum penetration depth of the indenter tip and $S$ is the contact stiffness, determined experimentally as the initial slope of the linear portion of unloading curve (see Figure 7.13). From the measured values of $S$, ${P}_{max}$ and ${\delta }_{max}$, we use Equations 7.44 and 7.45 to determine $A$. The reduced modulus is then obtained from the following expression for the contact stiffness, assuming a value for the contact stiffness that is the same for a circular contact of the same area:$\begin{array}{cc}S=\frac{2}{\sqrt{\pi }}{E}_{r}\sqrt{A}& \left(7.46\right)\end{array}$
Figure 7.13: Typical load-displacement curve for indentation of the polyester resin used to embed the paint samples, labeled to illustrate the values of${P}_{max}$,${h}_{max}$ and$S$.

## 8 Fracture

The stress-strain behavior for a many material can exhibit a range of phenomena, depending on the temperature. This is particularly true of many polymers, which can show the range of behaviors in a uniaxial tensile test shown in Figure 8.1. While not all of these behaviors are necessarily observed in the same material, the following general regimes can often be identified, based on 4 different temperature regimes (${T}_{1}$, ${T}_{2}$, ${T}_{3}$ and ${T}_{4}$).
• T1: Brittle behavior. This is generally observed at sufficiently low temperatures.
• T2: Ductile behavior (yield before fracture)
• T3: cold drawing (stable neck)
• T4: uniform deformation
Figure 8.1: Typical generic temperature behavior at different temperatures.
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

Different fracture modes are defined by the relationship between the applied stress and the crack geometry. These are illustrated schematically in Figure 8.2 Fracture of a homogeneous material fracture generally occurs under Mode I conditions, and this is the most important condition. Mode II conditions, where a shear stress is applied in the direction perpendicular to the crack front, is often important for interfacial fracture, including the adhesive bonding of materials with different properties. Mode III is generally not important for our purposes.
Figure 8.2: Fracture Modes .

### 8.2 Stress Concentrations

In the previous section on contact mechanics we introduced the concept of the energy release rate, $G$, which can be viewed as the driving force for crack propagation. Failure occurs when $G$ exceeds a critical value, ${G}_{c}$. This energy-based approach was originally formulated by Griffith, and is referred to as the
Griffith model
for this reason. We can also describe the driving force for crack propagation in terms of the detailed stress field in the vicinity of the tip of a propagating crack. This approach was developed by Irwin, and is referred to here as the
Irwin model
. The key concept here is that stresses are enhanced, or 'concentrated' in the vicinity of a defect like a crack. The easiest way to start thinking about this is to look at the nature of the stress distribution around a circular hole in a two-dimensional plate (Figure 8.3). A stress is a force per unit area, so we can imagine dividing up the stress into individual force lines, which are equidistant when the stress is uniform. Near a defect the lines of force are closer to one another, indicating that the stress is higher in this area. The maximum tensile stress at the sides of the hole is three times the average applied stress.
Figure 8.3: Force lines around a circular defect.
For an ellipse of with axis ${a}_{c}$ perpendicular to the applied stress and axis ${b}_{c}$ parallel to the applied stress (see Figure 8.4), the maximum stress in this case is given by the following expression:
$\begin{array}{cc}{\sigma }_{max}={\sigma }_{0}\left(1+2\frac{{a}_{c}}{{b}_{c}}\right)& \left(8.1\right)\end{array}$
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:
Combination of Eqs. 8.1 and 8.2 gives:
$\begin{array}{cc}{\sigma }_{max}={\sigma }_{0}\left(1+2\sqrt{{a}_{c}/{\rho }_{c}}\right)& \left(8.3\right)\end{array}$
$\begin{array}{cc}{\sigma }_{max}\propto {\sigma }_{0}\sqrt{{a}_{c}}& \left(8.4\right)\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.
Figure 8.4: Elliptical crack with a crack tip radius of curvature,${\rho }_{c}$.

### 8.3 Stress Intensity Factor

Consider a planar crack in the x-z plane, as shown conceptually Figure 8.5. The stress in the vicinity of the crack tip can be expressed in the following form:
$\begin{array}{cc}\sigma =\frac{K}{\sqrt{2\pi d}}f\left(\theta \right)& \left(8.5\right)\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).
(a)
(b)
Figure 8.5: Cartesian (a) and polar (b) coordinate axes use d to define stresses in the vicinity of a crack tip.

$\begin{array}{cc}\left(\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{12}\end{array}\right)=\frac{{K}_{I}}{\sqrt{2\pi d}}cos\frac{\theta }{2}\left(\begin{array}{c}1-sin\frac{\theta }{2}sin\frac{3\theta }{2}\\ 1+sin\frac{\theta }{2}sin\frac{3\theta }{2}\\ cos\frac{3\theta }{2}sin\frac{\theta }{2}\end{array}\right)& \left(8.6\right)\end{array}$
This compact notation is used to specify the three relevant values of $f\left(\theta \right).$ 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(1-sin\frac{\theta }{2}sin\frac{3\theta }{2}\right)& \left(8.7\right)\end{array}$
These expressions assume that the crack tip is very sharp, with a very small radius of curvature,${\rho }_{c}$. If $d$ is comparable to ${\rho }_{c}$, these equations no longer apply. Consider for example, the presence of an internal crack of length ${a}_{c}$ and radius of curvature ${\rho }_{c}$ in a thin sheet of material, shown schematically in Figure 8.6. In this case the stress at the crack edge is ${\sigma }_{max}$ as given by Eq. 8.3. An assumption in the use of Eq. 8.6 is that the stresses are substantially less than ${\sigma }_{max}$. In other words, $K$ describes the stress field close to the crack tip, but still at distances away from the crack tip that are larger than the crack trip radius of curvature, ${\rho }_{c}$.
The mode I stress intensity factor for this geometry is given by the applied stress, ${\sigma }_{0}$ and the crack length ${a}_{c}$:
Figure 8.6: An internal crack in a homogeneous solid.

$\begin{array}{cc}\left(\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{12}\end{array}\right)=\frac{{K}_{II}}{\sqrt{2\pi d}}\left(\begin{array}{c}-sin\frac{\theta }{2}\left(2+cos\frac{\theta }{2}cos\frac{3\theta }{2}\right)\\ sin\frac{\theta }{2}cos\frac{\theta }{2}cos\frac{3\theta }{2}\\ cos\frac{\theta }{2}\left(1-sin\frac{\theta }{2}sin\frac{3\theta }{2}\right)\end{array}\right)& \left(8.9\right)\end{array}$

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 stress-based 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:
$\begin{array}{cc}{K}_{IC}={\sigma }_{f}\sqrt{\pi {a}_{c}}& \left(8.10\right)\end{array}$Rearranging gives:
The fracture toughness, ${K}_{IC}$ has strange units - a stress multiplied by the square root of a length. In order to understand where this characteristic stress and the characteristic length actually come from, we need to consider the actual shape of the crack tip. Using Eq. 8.3 we see that the maximum stress in front of the crack tip, ${\sigma }_{max}^{f}$ , at the point of fracture is: $\begin{array}{cc}{\sigma }_{max}^{f}\approx 2{\sigma }_{f}\sqrt{{a}_{c}/{\rho }_{c}}& \left(8.12\right)\end{array}$
where we have assumed that $\sqrt{{a}_{c}/{\rho }_{c}}\gg 1$, so that we can ignore the extra factor of 1 in Eq. 8.3. Now we can use Eq. 8.10 to substitute ${K}_{IC}$ for ${\sigma }_{f}$. After rearranging we get:
$\begin{array}{cc}{K}_{IC}\approx {\sigma }_{max}^{f}\frac{\sqrt{\pi }}{2}\sqrt{{\rho }_{c}}\approx {\sigma }_{max}^{f}\sqrt{{\rho }_{c}}& \left(8.13\right)\end{array}$
This expression is really only valid for a crack tip with a well-defined 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 stress-based 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}& \left(8.15\right)\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.
The fact that $G\propto {K}_{I}^{2}$ for a mode I fracture experiment is illustrated in Figure 8.7, which we use to show the relationship between stress and stored elastic energy for an elastically deformed sample. The energy input to the sample up to the point of fracture, which we refer to as${U}_{f}$, is the area under the stress strain curve:
$\begin{array}{cc}{U}_{f}=\frac{1}{2}{\sigma }_{f}{e}_{f}=\frac{1}{2}\frac{{\sigma }_{f}^{2}}{{E}_{r}}& \left(8.16\right)\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}}& \left(8.17\right)\end{array}$
This is consistent with Eq. 8.14, but we need to do a more detailed analysis to get the prefactor exactly right.
Figure 8.7: Schematic stress/strain curve for a brittle material in the presence of a crack.

### 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.
Figure 8.8: Double cantilever beam geometry.
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}& \left(8.19\right)\end{array}$
So we have:
We now combine Eqs. 8.18 and 8.20 to obtain the following for the energy release rate:
$\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}}& \left(8.21\right)\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}}& \left(8.22\right)\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):
$\begin{array}{cc}\frac{{\sigma }_{zz}}{{\sigma }_{avg}}=0.5{\left(1-{\left(r/a\right)}^{2}\right)}^{-1/2}& \left(8.23\right)\end{array}$
Here the average stress,${\sigma }_{avg}$ is defined as follows:
$\begin{array}{cc}{\sigma }_{avg}\equiv \frac{P}{\pi {a}^{2}}& \left(8.24\right)\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 a-r& \left(8.25\right)\end{array}$
$\begin{array}{cc}\frac{{\sigma }_{zz}}{{\sigma }_{avg}}=\frac{1}{2}{\left[2d/a-{\left(d/a\right)}^{2}\right]}^{-1/2}& \left(8.26\right)\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}& \left(8.27\right)\end{array}$
by comparing to Eq. 8.6 for${K}_{I}$ we obtain the following: $\begin{array}{cc}{K}_{I}=\frac{1}{2}{\sigma }_{avg}{\left(\pi a\right)}^{1/2}& \left(8.28\right)\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).
Figure 8.9: Comparison of the full solution for the flat punch contact stresses (Eq. 8.23) with the${d}^{-1/2}$ singularity obtained from${K}_{I}$ (Eq. 8.28).

#### 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.
Figure 8.10: Dependence of the stress distribution under a flat punch for different values of the confinement ration,$a/h$.
• Max. stress in center for thin, incompressible layers ($\nu =0.5\right)$.
• 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:
Note that the prefactor in this expression is slightly different than what is given in Eq. 8.8 because of the different crack geometries. Eq. 8.8 is for a rectangular crack and Eq. 8.30 is for a circular crack. The energy release rate for the circular crack is given by using the relationship between ${K}_{I}$ and $G$ valid for a mode I crack at the interface between compliant and rigid materials (Eq. 8.29):
$\begin{array}{cc}G=\frac{{K}_{I}^{2}}{2{E}_{r}}=\frac{2{a}_{c}{\sigma }_{zz}^{2}}{\pi {E}_{r}}& \left(8.31\right)\end{array}$

### 8.7 Fracture Toughness of Materials

In the Griffith (energy-based) 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}}& \left(8.32\right)\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 1-2\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}}}& \left(8.33\right)\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}}& \left(8.34\right)\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 (non-crystalline 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}\right)$? 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}}}^{○}$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.
Figure 8.11: a) Schematic of a solid oxide fuel cell utilizing a zirconia solid electrolyte [1]. b) Image of a fuel cell stack, with several fuel cells stacked in series with one another.
An additional problem with the use of zirconia is that different phases are present at equilibrium at different temperatures. In its pure form, zirconia has three different crystal phases: a monoclinic phase at equilibrium at room temperature, a tetragonal phase at equilibrium at between 1170 ${}^{○}$C and 2370 ${}^{○}$C, and a cubic phase at equilibrium above 2370 ${}^{○}$C. These structures are shown schematically in Figure 8.12. The phases have different densities with the cubic phase having the highest density (smallest volume for a given mass of material) and with the monoclinic phase having the lowest density (highest volume per a given mass of material). Volume changes occurring during processing, when the sample is cooled from high temperatures where the cubic phase is stable, give rise to substantial cracking of the material, making pure zirconia unprocessible. A common solution to this problem is to add an alloying element that stabilizes the cubic phase to lower temperatures. Yttrium is a common alloying element, and can be added in its oxidized form (${Y}_{2}{O}_{3},$or Yttria), to form yttria stabilized zirconia (YSZ ).
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 higher-volume, 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.
Figure 8.13: $Zr{O}_{2}/{Y}_{2}{O}_{3}$ Phase Diagram (from ref. [2]).

#### 8.8.2 Tempered Glass

Figure 8.15: Comparison of tempered glass and float glass samples that have been fractured. (from http://www.bharatsafety.com/tempered_glass.html).
Figure 8.16: Laminated safety glass.

#### 8.8.3 Fiber Composites

Figure 8.17: Schematic representation of crack propagation through a fiber-reinforced composite. (from http://www.aml.engineering.columbia.edu/ntm/level2/ch05/html/l2c05s02.html)
Figure 8.18: use of Composites in the Boeing 787 Dreamliner
Figure 8.19: Views of a carbon fiber/epoxy composite perpendicular (left) and parallel (right) to the fiber direction.(from http://www.scielo.br/scielo.php?pid=S1516-14392010000300022&script=sci_arttext).
Figure 8.20: Laminate plies in a carbon fiber composite.
Figure 8.21: SEM image of failure of crack propagation in a carbon fiber composite (left), and a corresponding schematic representation of the crack path (right).

#### 8.8.4 Nondestructive Testing of a Fiber Laminate Composite

Figure 8.22: Nondestructive evaluation of damage in a carbon fiber composite.

## 9 Weibull Analysis of Failure

Failure of brittle materials is determined by the largest flaw size, since the largest flaw size will have the largest value of $K$ for a given applied stress. As an example of the use of Weibull statistics, consider the adhesive transfer of a thin layer of a material from a flat, flexible substrate to a rigid, curved indenter as illustrated in Figure 9.1. The basic geometry of the experiment is illustrated in Figure 9.1a, and consists of a thin, viscoelastic film that is coated on an elastomeric substrate. A hemispherical glass indenter is brought into contact with the film and is then pulled away from the surface. The system is designed so that the adhesion of the film to the glass indenter is stronger than the adhesion to the elastomeric substrate, the film will be transfered from the substrate to the indenter. The process occurs by the sequence of steps illustrated in parts b-e of Figure 9.1:
• 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.
Figure 9.1: Adhesive transfer of a thin, viscoelastic film.
Details of this experiment can be found in reference[12]. The most important thing to keep in mind is that the whole transfer process is controlled by the initial appearance of a cavity at the indenter/substrate interface (Figure 9.1b). This happens when ${p}_{max}$, the maximum hydrostatic tension at the film/substrate interface, reaches a critical value that we refer to as ${p}_{cav}$. Qualitatively we find that ${p}_{cav}$ is close to ${E}_{sub}$, the elastic modulus of the substrate, but cavitation occurs for different values of ${p}_{cav}$. The Weibull distribution for the survival probability, ${P}_{s}$ (the probability that cavitation has NOT occurred) is as follows:
$\begin{array}{cc}ln\left[ln\left(1/{P}_{s}\right)\right]=Mln{p}_{max}-Mln{p}_{cav}& \left(9.2\right)\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:
1. 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.
2. Organize this list from the lowest value of ${p}_{max}$ to the highest value.
3. 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}.$ 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%:
$\begin{array}{cc}{P}_{s}=exp\left[-{\left(\frac{\sigma }{{\sigma }_{avg}}\right)}^{M}\right]& \left(9.3\right)\end{array}$
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 1-x$ 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

Factors that increase the hardness or tensile strength (${\sigma }_{ts}\right)$ of a material generally decrease the toughness, ${K}_{IC}$, so that the possible values of these two materials lie along a curve like the one shown in Figure 10.1. This type of relationship exists a certain degree of yielding is necessary in order to increase the characteristic crack tip radius of curvature, ${\rho }_{c}$ (recall that ${K}_{IC}\propto \sqrt{{\rho }_{c}}$ as described by Eq. 8.13. The mode I toughness of a material can be expressed in terms of a critical stress intensity factor, ${K}_{IC}$, or a critical energy release rate, ${G}_{C}$, which are related to one another through the following version Eq. 8.14:
$\begin{array}{cc}{G}_{C}=\frac{{K}_{IC}}{{E}_{r}}& \left(10.1\right)\end{array}$
It is common to consider the intrinsic toughness and the extrinsic toughness of a material. Extrinsic toughening is based on the modification of the stress field in the region of a growing crack. We can describe it by dividing the macroscopic stress intensity factor,${K}_{I}$, into two parts,${K}_{tip}$ and${K}_{s}$:
$\begin{array}{cc}{K}_{I}={\sigma }_{0}\sqrt{\pi a}={K}_{tip}+{K}_{s}& \left(10.2\right)\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}}& \left(10.3\right)\end{array}$

### 10.1 Hydraulic Fracturing (Fracking)

Figure 10.2: The hydraulic fracturing (fracking) process.

## 11 Yield Behavior

The yield point of a material corresponds to the onset of permanent deformation, originating for example from the movement of dislocations. The stress/strain curve for a simple tensile test is shown in Figure 11.1, with the tensile yield stress, ${\sigma }_{y}$, corresponding to the onset of permanent, irreversible deformation in the material. For most materials this corresponds to the onset of nonlinearity in the stress-strain curve (rubber is the exception, and that case is discussed in more detail in 331). What we need is a generalized criterion that can be used to determine the onset of yield for any stress state. These yield criteria all focus on the importance of the shear stress, and we begin with the discussion of yield in single crystals in the following Section .
Figure 11.1: Tensile testing sample and representative data for a ductile sample.

### 11.1 Critical Resolved Shear Stress

The simplest criterion for the yield of a material is that the resolved shear stress , ${\tau }_{RSS}$, exceeds a critical value referred to simply as the critical resolved shear stress , ${\tau }_{CRSS}$. The resolved shear stress is obtained from the applied stress and the orientation of the slip plane as illustrated in Figure 11.2. For a single crystal the relevant resolved shear stress is the shear stress on acting on the slip plane, in the direction of the Burgers vector. In a tensile experiment it is given by the applied tensile stress, $\sigma$, the angle $\phi$ between the slip plane normal and the tensile axis, and the angle $\lambda$ between the tensile axis and the slip direction:
$\begin{array}{cc}{\tau }_{RSS}=\sigma cos\phi cos\lambda & \left(11.1\right)\end{array}$
$\begin{array}{cc}{\sigma }_{y}=\frac{{\tau }_{CRSS}}{cos\phi cos\lambda }& \left(11.2\right)\end{array}$
Figure 11.2: Resolved shear stress.
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=⟨\frac{1}{cos\phi cos\lambda }⟩& \left(11.4\right)\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}& \left(11.5\right)\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 3-dimensionsal 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{{|{\sigma }_{i}^{p}-{\sigma }_{j}^{p}|}_{max}}{2}>{\tau }_{crit}& \left(11.6\right)\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}}& \left(11.7\right)\end{array}$
Yielding occurs when ${\sigma }_{e}>{\sigma }_{y}$. The Tresca and Von Mises yield surfaces for a two dimensional stress state (${\sigma }_{3}^{p}$=0) are shown in Figure 11.3. Yielding does not occur inside the surface, but does occur outside the surface.
Figure 11.3: Cross sections of the Tresca and Von Mises yield surfaces at the${\sigma }_{3}^{p}=0$ plane.
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:
Figure 11.4: Uniaxial tensile test.
$\begin{array}{cc}{\sigma }_{eng}=P/{A}_{0}& \left(11.8\right)\end{array}$$\begin{array}{cc}{\sigma }_{true}=P/A& \left(11.9\right)\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}& \left(11.10\right)\end{array}$
The relationship between the true stress and the engineering stress is therefore as follows:

#### 11.3.1 Considére Construction

The Considére construction is a simple construction that can be used to determine the stability of regions in a sample at large tensile deformations. It can be used in to distinguish between unstable and unstable necking of a sample, illustrated schematically in Figure 11.5. We begin by considering a region of the sample that has a slightly thinner cross section than the rest of the sample. The true stress in this region of the sample will be higher than the rest of the sample because we are dividing the applied load by a lower cross section . Two things can happen at this point, the first possibility is that the larger stress in this region of the sample leads to greater deformation, and the sample breaks as the necked region begins to thin down. This is the unstable necking condition illustrated on the left side of Figure 11.5. In this case the maximum force, $P$, applied to the sample is the force where the neck begins to form.
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.
Figure 11.5: Schematic representation of stable and unstable necking of a sample under tensile loading conditions.
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& \left(11.13\right)\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& \left(11.14\right)\end{array}$
which we rearrange to the following:
$\begin{array}{cc}\frac{d{\sigma }_{t}}{d\lambda }=\frac{{\sigma }_{t}}{\lambda }& \left(11.15\right)\end{array}$
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 stress-strain curve (point A in Figure 11.6). A necking instability forms when the true stress reaches this value, resulting in a thinned-down 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 stress-strain curve lies below the tangent line) and the tangent at point B represents a minimum in the applied force (the stress-strain 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.
$\begin{array}{cc}\frac{d{\sigma }_{t}}{d\lambda }=\frac{{\sigma }_{t}}{\lambda }& \left(11.16\right)\end{array}$

#### 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 Table6 for a variety of metals.
Table 6: Strain Hardening Coefficients for Various Materials (From Hertzberg, Table 2.8).
 Material Strain Hardening Coefficient Stainless Steel 0.45-0.55 Brass 0.35-0.4 Copper 0.3-0.35 Aluminum 0.15-0.25 Iron 0.05-0.15

#### 11.3.4 Necking Instability in a Power-Law Strain Hardening Material

If we write the true stress appearing in Eq. 11.17 as $ln\lambda$ (see Eq. 4.37), we have:
$\begin{array}{cc}{\sigma }_{t}=K{\left(ln\left(\lambda \right)\right)}^{n}& \left(11.18\right)\end{array}$
Differentiating with respect to $\lambda$ gives:
$\frac{d{\sigma }_{t}}{d\lambda }=\frac{nK{\left(ln\left(\lambda \right)\right)}^{n-1}}{\lambda }=\frac{nK{e}_{t}^{n-1}}{\lambda }$
and then using Eq. 11.15 to equate $d{\sigma }_{\lambda }/d\lambda$ to$\sigma /\lambda$ gives the following:
$\frac{nK{e}_{t}^{n-1}}{\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}& \left(11.19\right)\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

Figure 11.7: Example of a series of twin boundaries in an FCC crystal (from ref.[10], adapted with permission from Macmillan Publishers Ltd.).
Figure 11.8: Undeformed material (a) and deformed material after slip (b) and deformation twinning (c).
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:
1. We can make a material with a very low dislocation density, thereby eliminating dislocation as a yield mechanism.
2. 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 Frank-Read dislocation source in Figure 12.1. (See the 316-1 text to review Frank-Read 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 316-2. 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.
Figure 12.1: Schematic illustration of a Frank-Read source.

### 12.1 Review of some Dislocation Basics

The basics of dislocations are covered in 316-1. Here are the most important concepts to remember and use:
1. Dislocations correspond to a discontinuity of the displacement field. This continuity is quantified by the Burgers vector, $\stackrel{\to }{b}$.
2. Dislocations move along slip planes that contain both $\stackrel{\to }{b}$ and the dislocation line.
3. When a dislocation exits the sample, portions of the sample on either side of the slip plane are displaced by $\stackrel{\to }{b}$.
4. The energy per length of a dislocation, ${T}_{s}$, is proportional to $G{b}^{2}$, where $b$ is the magnitude of $\stackrel{\to }{b}$. This energy per unit length has the dimensions of a force, and can be viewed as a 'line tension' acting along the dislocation.
5. The resolved shear stress required to bend a dislocation through a radius of curvature of $r$ is ${T}_{s}/rb$.
6. 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 $\stackrel{\to }{{b}_{1}}$ and ${\stackrel{\to }{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}\stackrel{\to }{\tau }=\frac{G\stackrel{\to }{b}}{2\pi d}& \left(12.1\right)\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 $\stackrel{\to }{b}$.
This stress induces a force on the dislocation, given by the following expression:
If ${\stackrel{\to }{b}}_{1}$ and ${\stackrel{\to }{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.
Figure 12.2: Interactions between screw dislocations of the same sign (solid line) and opposite sign(dashed line).
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{\stackrel{\to }{b}}_{2}\cdot {\stackrel{\to }{b}}_{1}}{2\pi \left(1-\nu \right)d}& \left(12.3\right)\end{array}$
The expression is the same as the expression for the screw dislocations, with the extra factor of $1-\nu$.
If the two edge dislocations do not lie on the same slip plane the situation is more complicated, but we can still make some qualitative arguments based on the nature of the strain fields around the dislocations. Consider two edge dislocations on glide planes that are separated by a distance $h$, as illustrated in Figure 12.3. In this case the Dislocations begin to line up due to cancellation of the strains in the regions where these strain field overlap. The tensile strain induced below the upper dislocation is partially compensated for by the compressive strain above the lower dislocation. As a result the dislocations move along their respective glide planes so that they lie on top of one another.
Figure 12.3: Overlapping stress fields for two adjacent dislocations of the same sign.
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.
Figure 12.4: Map showing the regimes where two edge dislocations move toward one another with attractive interaction (-) or apart from one another with a repulsive interaction (+). One of the dislocations is at the origin, with $\stackrel{ˆ}{s}$ pointing into the plane of the figure. The other dislocation is assumed to have the same value for $\stackrel{ˆ}{s}$ and $\stackrel{\to }{b}$. Both dislocations move in glide planes that are perpendicular to the y axis. These dislocations move toward each other if the second dislocation is located within a region that is attractive (-), and apart from one another if the second dislocation is located in a region that is repulsive (+).
Dislocations can also interact with point defects or other atoms. The reason for this is that point defects also distort the lattice, giving strain fields that overlap with the strain field of a dislocation. This effect is easiest to visualize for edge dislocations, as illustrated in Figure 12.5. Substitutional impurities with different sizes than the majority component atoms can reduce the strain fields surrounding the dislocation by moving to the appropriate region near the dislocation. Large atoms will tend to substitute for atoms in regions where the local strain is tensile, and small atoms will tend to substitute for atoms in regions where the local strain is compressive. Interstitial impurities will segregate to regions of tensile stress as well.
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

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.
Figure 12.6: Qualitative relationship between strength and dislocation density.

### 12.4 Grain Boundary Strengthening

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 Hall-Petch equation is commonly used to describe the relationship between the yield strength, ${\sigma }_{y}$ and the grain size, $d$:
$\begin{array}{cc}{\sigma }_{y}={\sigma }_{i}+{k}_{y}{d}^{-1/2}& \left(12.4\right)\end{array}$
Here${\sigma }_{i}$ and ${k}_{y}$ are factors determined by fitting to experimental data. Note that the Hall-Petch equation is based on experimental observations, and is not an expression derived originally from theoretical considerations.

### 12.5 Precipitation Hardening

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 Frank-Read source in Figure 12.1. We know from 316-1 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

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.
Figure 13.1: Absolute melting temperatures (K) for the elements.
Typical creep response for a metal is shown in Figure 13.2. The easiest way to do the experiment is under constant load conditions. Because the cross sectional area of the sample decreases as the material deforms, the true stress increases during a constant load creep experiment. If the experiment is done under conditions where the load is reduced to maintain a constant true stress, the creep rate does not increase as fast at the longer times. The increase in the creep rate at long times for the constant load test can be attributed to the increased true stress in the sample at the later stages of the creep experiment. Failure of the sample during a constant load test is referred to as the creep rupture point. The creep rupture time decreases with increasing temperature and increasing applied engineering stress, as illustrated by the data for an iron alloy in Figure 13.3.
Figure 13.2: Qualitative difference between creep behavior under constant load and constant stress conditions.
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, ${\stackrel{˙}{\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, steady-state creep regime in the following discussion.
Figure 13.3: Creep rupture data for an iron-based alloy.
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:
$\begin{array}{cc}{\stackrel{˙}{\epsilon }}_{II}=f\left(\sigma \right)exp\left(-\frac{{Q}_{c}}{{k}_{B}T}\right)& \left(13.1\right)\end{array}$
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.
In the high temperature regime, the activation energy for creep, ${Q}_{c}$, is very similar to the activation energy for self-diffusion in the metal (see Figure 13.6), indicating that high temperature creep is controlled by the same processes that control self diffusion. As we describe below, this common process is vacancy diffusion. Vacancy diffusion can contribute to the creep behavior by enabling dislocation climb, or by changing the shapes of individual grains in a manner that does not require the existence of dislocations. These mechanisms are Nabarro-Herring creep and Coble Creep, as described below.
Figure 13.6: Activation energies for creep (for $T/{T}_{m}>0.5\right)$ and atomic diffusion (from ref.[14]).

### 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& \left(13.3\right)\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)& \left(13.4\right)\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 (Nabarro-Herring 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 Nabarro-Herring Creep: Bulk Diffusion Within a Grain

Nabarro-Herring 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.
Figure 13.7: Conceptual illustration of Nabarro-Herring creep.
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)& \left(13.5\right)\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}& \left(13.6\right)\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$:
Combining Eqs.13.7 and13.8 gives:
$\begin{array}{cc}{J}_{v}\approx -{D}_{v}\left(\frac{{C}_{v}^{t}-{C}_{v}^{c}}{d}\right)& \left(13.9\right)\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$:
$\begin{array}{cc}v=\Omega {J}_{v}\approx \Omega {D}_{v}\left(\frac{{C}_{v}^{t}-{C}_{v}^{c}}{d}\right)& \left(13.10\right)\end{array}$
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)& \left(13.11\right)\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)& \left(13.13\right)\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:
Combination of Eqs. 13.12-13.14 gives the following for the strain rate:
$\begin{array}{cc}\frac{de}{dt}\approx \frac{2{D}_{0}}{{d}^{2}}\frac{\sigma \Omega }{{k}_{B}T}exp\left(-\frac{{Q}_{vol}}{{k}_{B}T}\right)& \left(13.15\right)\end{array}$
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

If diffusion occurs predominantly along the grain boundaries, the mechanisms is slightly different and is referred to as Coble Creep. The expression for the strain rate is very similar to the expression for Nabarro-Herring creep, with two difference. First, we replace the bulk vacancy diffusion coefficient, ${D}_{v}$ in Eq. 13.7 with the grain boundary diffusion coefficient, ${D}_{gb}$, which is generally larger than ${D}_{v}$. In addition, we modify the relationship between ${J}_{v}$ to account for the fact that vacancies are not diffusing throughout the entire grain, but only along a narrow grain of width ${d}_{gb}$. As a result we need to multiply the velocity in Eq. 13.10 by the overall volume fraction of the grain boundary, which is proportional to ${d}_{gb}/d$. The net result is the following expression for the strain rate for Coble creep:
$\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)& \left(13.16\right)\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)& \left(13.17\right)\end{array}$
$\begin{array}{cc}{Q}_{gb}={Q}_{v}+{G}_{v}& \left(13.18\right)\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)& \left(13.19\right)\end{array}$
Here ${Q}_{gb}$ is less than ${Q}_{vol}$, so Coble creep will be favored over Nabarro-Herring 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 and13.2, but with a specific exponent, $n$, of 4.5.
Figure 13.8: Schematic representation of dislocation climb.
In the Weertman model the creep rate is controlled by dislocation climb, illustrated schematically in Figure 13.8. The following relationship for the strain rate was obtained from theoretical considerations:
$\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)& \left(13.20\right)\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 Nabarro-Herring (N-H) 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 Nabarro-Herring creep is a vertical line, since the stress dependence for both measurements is the same.
• Nabarro-Herring creep occurs at higher temperatures than Coble creep, because diffusion through grains (the Nabarro-Herring 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 Nabarro-Herring regions will move to the right.
• The line separating the Nabarro-Herring 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).
Figure 13.9: Deformation mechanism map.

## 14 Deformation of Polymers

The deformation mechanisms in polymeric materials are completely different from those in metals and ceramics, and (almost) never have anything to do with the motion of dislocations. To begin with, we can separate into the polymers that are partially crystalline and those that are not. As one would expect, the structures of these non-crystalline (amorphous) and semicrystalline polymers are very different, as shown schematically in Figure Crystallization and glass formation are the two most important concepts underlying the physical properties of polymers. The ways in which the molecules are organized in non-crystalline (amorphous) polymers and semicrystalline polymers are very different, as illustrated in Figure s Polymers crystallize at temperatures below ${T}_{m}$ (melting temperature) and form glasses at temperatures below ${T}_{g}$ (glass transition temperature). All polymers will form glasses under the appropriate conditions, but not all polymers are able to crystallize. The classification scheme shown in Figure 14.3 divides polymeric materials based on the locations of ${T}_{g}$ and ${T}_{m}$ (relative to the use temperature, $T$) and is a good place to start when understanding different types of polymers.
Figure 14.1: Structure of an amorphous polymer.
Figure 14.4: Load-Displacement and true stress-true strain curves for PVC (polyvinyl chloride) and HDPE (high density polyethylene).
Figure 14.5: Stress-strain curves for PMMA (polymethylmethacrylate) under tensile and compressive loading conditions.
Figure 14.6: Yield surface for PMMA at 20${}^{○}$C and at 90${}^{○}$C. For comparison a map of the Tresca yield criterion (where normal forces do not matter) shown as the dashed line.

#### 14.0.1 Case Study: Fracture toughness of glassy polymers

Figure 14.7: Molecular weight dependence of fracture toughness for polystyrene (PS) and poly(methyl methacrylate) (PMMA).
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:
(a)
(b)
Figure 14.8:
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. 1) new fibrils are created at the craze tip
2. 2) fibrils break to form a true crack at the crack tip
3. Crazing requires a stress field with a tensile hydrostatic component${\sigma }_{1}+{\sigma }_{2}+{\sigma }_{3}>0$ (crazes have voids between fibrils)
4. Crazing occurs first for PMMA in uniaxial extension (${\sigma }_{2}=0$)
5. ${G}_{Ic}$ is determined by energy required to form a craze ($\approx 1000\phantom{\rule{6px}{0ex}}J/{m}^{2}$)
6. Crazing requires strain hardening of fibrils - material must be entangled ($M>{M}_{c}$),${M}_{c}$ typically$\approx$30,000 g/mol.
7. 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.
Figure 14.10: Meniscus instability mechanism for the formation of craze fibrils.
http://n-e-r-v-o-u-s.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)
Figure 14.11: Deformation map for the shear yielding and crazing for plane stress conditions (${\sigma }_{3}^{p}=0\right)$
##### 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}& \left(14.1\right)\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.
Figure 14.12: Crack tip stresses in the Dugdale model.
###### High Impact Polystyrene:
Polystyrene (PS) is a big business - how do we make it tougher? High impact polystyrene (HIPS) is a toughened version of polystyrene produced by incorporating small, micron-sized rubber particles in the material. The morphology is shown in Figure 14.13, and consists small PS inclusions embedded in rubber particles that are in turn embedded in the PS matrix materials. The rubber particles act as stress concentrators that act as nucleation points for crazes.
Figure 14.13: Morphology of high impact polystyrene.
Figure 14.14: Multiple crazes in high impact polystyrene.
###### 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.
Figure 14.15: Schematic stress-strain curves for polystyrene (PS) and high impact polystyrene (HIPS) in the absence of a crack.
deformation via crazing in vicinity of rubber particles (stress concentrators) throughout sample
###### Samples with Precrack:
(measurement of KIC or GIC)
Figure 14.16: Schematic stress-strain curves for polystyrene (PS) and high impact polystyrene (HIPS) in the presence of a crack.
• Deformation limited to region around crack tip
• Much more deformation for HIPS - higher toughness
###### Impact Tests
Impact tests are designed to investigate the failure of materials at high rates. Two common standardized are the Izod impact test and the Charpy impact test . They both involve measuring the loss in kinetic energy of a swinging pendulum as it fractures a sample. The geometry of a Charpy impact test is illustrated in Figure 14.17 Decrease in pendulum velocity after breaking sample gives impact toughness. For a useful discussion of a Charpy impact test, see https://www.youtube.com/watch?v=tpGhqQvftAo.
Figure 14.17: Charpy impact test.
Figure 14.18: Izod impact geometry.

## 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 time-dependent 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 time-dependent relaxation modulus , $E\left(t\right)$:
$\begin{array}{cc}E\left(t\right)=\frac{\sigma \left(t\right)}{{e}_{0}}& \left(15.1\right)\end{array}$
We are often interested in the application of an oscillatory strain to a material. Examples include the propagation of sound waves, where wave propagation is determined by the response of the material at the relevant frequency of the acoustic wave that is propagating through the material. In an oscillatory experiment, referred to as a dynamic mechanical experiment in Figure 15.1, the applied shear strain is an oscillatory function with an angular frequency, $\omega$, and an amplitude, ${e}_{0}$:
$\begin{array}{cc}e\left(t\right)={e}_{0}sin\left(\omega t\right)& \left(15.2\right)\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${}^{○}$:
$\begin{array}{cccc}\frac{de}{dt}& =& {e}_{0}\omega cos\left(\omega t\right)={\gamma }_{0}\omega sin\left(\omega t+\pi /2\right)& \left(15.3\right)\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)& \left(15.4\right)\end{array}$
Now we can define a complex modulus with real and imaginary components as follows:
$\begin{array}{cc}{E}^{*}={E}^{\prime }+i{E}^{\prime \prime }=|{E}^{*}|{e}^{i\phi }& \left(15.5\right)\end{array}$
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,$|{E}^{*}|$ and phase,$\phi$. The magnitude of the complex modulus is simply the stress amplitude normalized by the strain amplitude:
$\begin{array}{cc}|{E}^{*}|={\sigma }_{0}/{e}_{0}& \left(15.6\right)\end{array}$
Figure 15.2: Time dependent stress and strain in a dynamic mechanical experiment.
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:
We can combine Eqs. 15.8 and 15.9 to get the following expression for $tan\phi$, commonly referred to simply as the loss tangent:
$\begin{array}{cc}tan\phi =\frac{{E}^{\prime \prime }}{{E}^{\prime }}& \left(15.10\right)\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

We can now consider the specific example of a torsional resonator, shown below in Figure 15.3. The equation of motion for the spring as it is being twisted is[21]:
$\begin{array}{cc}T={K}_{t}\theta +I\frac{{d}^{2}\theta }{d{t}^{2}}& \left(15.11\right)\end{array}$
$\begin{array}{cc}{K}_{t}\equiv \frac{T}{\theta }=\frac{\pi G{d}^{4}}{32\ell }& \left(15.12\right)\end{array}$
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}}& \left(15.13\right)\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:
$\begin{array}{cc}\theta ={\theta }_{0}cos\left({\omega }_{n}t\right)& \left(15.14\right)\end{array}$
$\begin{array}{cc}\theta ={\theta }_{0}Real\left(exp\left(i{\omega }_{n}t\right)\right)& \left(15.15\right)\end{array}$
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 }& \left(15.16\right)\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)& \left(15.17\right)\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)& \left(15.18\right)\end{array}$
So the imaginary part of the complex frequency describes the decay of oscillation with time.

### 15.3 Viscoelastic Models

Models of viscoelasticity can be represented visually by connecting solid-like and liquid-like elements together. In materials science we're used to the solid-like elements. They are simply elastic springs where the force is proportional to the extension of the spring. In terms of stress and strain, the stress, $\sigma$ is proportional to the strain, $e$. There is not time-dependence to the behavior of an ideal spring. The stress is proportional to the strain, regardless of how fast or slow the strain is applied. In many real materials the rate at which the strain is applied matters as well, and this is where the dashpots come in. A dashpot is a liquid-like element where the stress is proportional to the rate at which the strain is changing. The proportionality constant is the viscosity , $\eta$, as defined in the following expression:
$\begin{array}{cc}\sigma =\eta \frac{de}{dt}& \left(15.19\right)\end{array}$
Figure 15.4: Schematic representations of springs and dashpots used to represent the time-dependent response of viscoelastic materials.
As an simple illustration of the significance of the viscosity, we can calculate the time it takes for a small sphere of metal to drop to the bottom of pool of water. Suppose we use an iron particle with a radius of 1$\mu$m. The situation is as shown in Figure 15.5. The gravitational force, F, causing the ball to sink is proportional to the volume of the sphere, and the difference in densities between the solid and liquid:
Figure 15.5: Terminal velocity of a sphere descending in a viscous fluid.
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:
At steady state, the velocity reaches a value that is obtained by equating the forces in Eqs.15.20 and15.21. In this way we obtain:
$\begin{array}{cc}V=\frac{2}{9}\frac{g{R}^{2}}{\eta }\left({\rho }_{s}-{\rho }_{\ell }\right)& \left(15.22\right)\end{array}$
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}$ Pa-sec. From this we get a velocity of 16$\mu$m/s.

#### 15.3.1 Maxwell model

The simplest viscoelastic model that contains both liquid-like and solid-like elements is the Maxwell model consisting of a linear dashpot with viscosity $\eta$ in series with a linear spring with a modulus, $E$ (see Figure 15.6). Because the elements are in series the strain is the same in each one of them. This stress can be related to the strain ${e}_{1}$ in the dashpot and the strain ${e}_{2}$ in the spring through the following expressions:
$\begin{array}{cc}\sigma =\eta \frac{d{e}_{1}}{dt}& \left(15.23\right)\end{array}$
$\begin{array}{cc}\sigma =E{e}_{2}& \left(15.24\right)\end{array}$
Figure 15.6: Maxwell model.
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)& \left(15.26\right)\end{array}$
with $\tau =\eta /E$. We divide by ${e}_{0}$ to obtain the relaxation modulus: $\begin{array}{cc}E\left(t\right)\equiv \sigma \left(t\right)/{e}_{0}=Eexp\left(-t/\tau \right)& \left(15.27\right)\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:
We can obtain the real and imaginary parts of ${E}^{*}$ by multiplying the top and bottom of Eq. 15.28 by $1+i/\left(\omega \tau \right)$, remembering that ${i}^{2}=-1$. Taking the real and imaginary components of${E}^{*}$ gives${E}^{\prime }$ and${E}^{\prime \prime }$:
$\begin{array}{cc}{E}^{\prime }\left(\omega \right)=\frac{E{\left(\omega \tau \right)}^{2}}{1+{\left(\omega \tau \right)}^{2}}& \left(15.29\right)\end{array}$
$\begin{array}{cc}{E}^{\prime \prime }\left(\omega \right)=\frac{E\omega \tau }{1+{\left(\omega \tau \right)}^{2}}& \left(15.30\right)\end{array}$

#### 15.3.2 Standard Linear Solid

For a single maxwell element the relaxation modulus decays to zero at long times ($t\gg \tau$). In many real systems we are interested in describing the day of the relaxation modulus from a large value at short times to a much smaller value at larger times. The standard linear solid accounts for this by putting an elastic element with a modulus of ${E}_{r}$ (the 'relaxed' modulus). In parallel with a single Maxwell element as shown in Figure 15.7. Extra spring adds directly to stress response, so all we have to do is add ${E}_{r}$ to the expressions for $E\left(t\right)$ and ${E}^{*}\left(w\right)$ that were obtained from the Maxwell model:
$\begin{array}{cc}E\left(t\right)={E}_{r}+{E}_{1}exp\left(-t/{\tau }_{1}\right)& \left(15.31\right)\end{array}$
$\begin{array}{cc}{E}^{*}\left(\omega \right)={E}_{r}+\frac{{E}_{1}i\omega {\tau }_{1}}{1+i\omega {\tau }_{1}}& \left(15.32\right)\end{array}$
Figure 15.7: Standard linear solid.
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

For real systems, the behaviors of $E\left(t\right)$ and ${E}^{*}\left(\omega \right)$ are usually much more complicated than given by the predictions of the Maxwell model or standard linear solid. To describe the behavior for real materials, we can add an arbitrary number of Maxwell elements in parallel, resulting in the generalized Maxwell model shown in Figure 15.8. The stresses for each of the parallel elements are additive, so $E\left(t\right)$ and ${E}^{*}\left(w\right)$ are given by the following expressions:$\begin{array}{cc}E\left(t\right)={E}_{r}+\sum _{j}{E}_{j}exp\left(-t/{\tau }_{j}\right)& \left(15.33\right)\end{array}$
$\begin{array}{cc}{E}^{*}\left(\omega \right)={E}_{r}+\sum _{j}\frac{{E}_{j}i\omega {\tau }_{j}}{1+i\omega {\tau }_{j}}& \left(15.34\right)\end{array}$
with${\tau }_{j}={\eta }_{j}/{E}_{j}$.
Figure 15.8: Generalized Maxwell model.

#### 15.3.4 Kelvin-Voigt Model

The Kelvin-Voigt model consists of a spring and dashpot in parallel, as shown in Figure 15.9. In this case the stresses in two elements are additive and the strains in the two elements are the same:$\begin{array}{cc}\sigma =\eta \frac{de}{dt}+Ee& \left(15.35\right)\end{array}$
The Kelvin-Voigt 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(1-exp\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& \left(15.36\right)\end{array}$The dynamic modulus for the Kelvin-Voigt element is given by the following expression:
$\begin{array}{cc}{E}^{*}=E+i\omega \eta & \left(15.37\right)\end{array}$
So that ${E}^{\prime }$ is simply equal to $E,$ and ${E}^{\prime \prime }$ is equal to $\omega \eta$.
Figure 15.9: Kelvin-Voigt Model.

## 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

An example of the sort of spring/dashpot model used to describe creep behavior of a linear viscoelastic material is shown in Figure 16.1. The model consists of a single Kelvin-Voigt element in series with a spring and a dashpot. In this example the strain obtained in response to a jump in the stress from 0 to $\sigma$ at $t=0$ can be represented as the sum of three contributions: ${e}_{1}$, ${e}_{2}$, and ${e}_{3}$:
$\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)& \left(16.2\right)\end{array}$
Figure 16.1: Linear creep model.
$\begin{array}{cc}\begin{array}{c}{e}_{1}=\frac{\sigma }{{E}_{1}}\\ {e}_{2}=\frac{\sigma }{{E}_{2}}\left[1-exp\left(-t/{\tau }_{2}\right)\right]\\ {e}_{3}=\frac{\sigma }{{\eta }_{3}}t\end{array}& \left(16.3\right)\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 no-longer possible to define a stress-independent 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)& \left(16.4\right)\end{array}$
1. 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 time-dependent functions works (at least for these two stress levels). These experiment enable us to obtain the time dependent function $J\left(t\right)$.
2. To get the stress-dependent function, $f\left(\sigma \right)$, it is sufficient to make a series of measurements at a single, experimentally convenient time.
Figure 16.3: Creep response of a material with separable dependencies on stress and 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}_{1}=\frac{\sigma }{{E}_{1}}\\ {e}_{2}={C}_{1}{\sigma }^{n}\left[1-exp\left(-{C}_{2}t\right)\right]\\ {e}_{3}={C}_{3}{\sigma }^{n}t\end{array}& \left(16.5\right)\end{array}$
$\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}& \left(16.6\right)\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.
Figure 16.4: Nonlinear creep model.

### 16.4 Eyring Model of Steady State Creep

#### 16.4.1 Material Deformation as a Thermally Activated Process

Our starting point is to realize that material deformation is a thermally activated process, meaning that there is some energy barrier that needs to be overcome in order for deformation to occur. The general idea is illustrated in Figure 16.5. The stress does an amount of work on the system equal to $\sigma v$, where $\sigma$ is the applied stress and $v$ is the volume of the element that moves in response to this applied stress. The quantity $v$ is typically referred to as an activation volume . The net result of the application of the stress is to reduce the activation barrier for motion in the stress direction by an amount equal to $v\sigma /2$ and to increase the activation barrier in the opposite direction by this same amount.
##### 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)& \left(16.7\right)\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)& \left(16.8\right)\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]& \left(16.9\right)\end{array}$
$\begin{array}{cc}sinh\left(x\right)=\left({e}^{x}-{e}^{-x}\right)/2& \left(16.10\right)\end{array}$
we obtain the following expression for the strain rate:
Before considering the behavior of the Eyring rate equation for high and low stresses, it is useful to consider the overall behavior of the sinh function, which is illustrated in Figure 16.6. Note that for small $x$, $sinhx\approx x$, and for large x, $sinh\left(x\right)\approx 0.5exp\left(x\right)$.
Figure 16.6: Behavior of the hyperbolic sine function.
##### 16.4.1.2 Low Stress Regime
In the low-stress 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$.
Comparing Eqs.16.12 and16.13 gives the following for the viscosity:
$\begin{array}{cc}\eta =\frac{{k}_{B}T}{A{f}_{0}v}exp\left(\frac{{Q}^{*}}{{k}_{B}T}\right)& \left(16.14\right)\end{array}$
So the Eyring theory reduces to and Arrhenius viscosity behavior in the linear, low-stress regime.
##### 16.4.1.3 High Stress Regime
In the high-stress 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)& \left(16.16\right)\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$
Eyring rate models are most often used for polymeric systems. In this case the activation volume can be viewed as the volume swept out the portions of a polymer molecule which move during a fundamental creep event, as illustrated schematically in Figure 16.7. A large activation volume means that cooperative deformation of a large region of the material is required in order for the material to flow. Low values of the activation volume indicate that the deformation is controlled by a very localized event, corresponding, for example, to the rupture of a single covalent bond.
Figure 16.7: Illustration of the activation volume.

#### 16.4.2 Additional Nonlinear Dashpot Elements

Nonlinear elements based on the Eyring rate model can also be included in our spring and dashpot models of viscoelastic behavior. For example, we can use an Eyring rate model to describe the stress dependence of the steady state creep of a nonlinear model, in which case Figure 16.4 gets modified to Figure 16.8. The steady state component is specified by the prefactor $A{f}_{0}$ the activation energy ${Q}^{*},$ and the activation volume, $v$.
Figure 16.8: Nonlinear creep model, including an Eyring rate element for the steady-state creep component.

## 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.
Figure 17.1: Schematic representation of the space elevator.
Consider a mass, $m$, that is located a distance$r$ from the center of the earth, as illustrated in Figure 17.2. The net force on the object is a centripetal force acting outward (positive in our sign convention) and a gravitational force acting inward (negative in our sign convention):
$\begin{array}{cc}F=m{\omega }^{2}r-{G}_{gr}{M}_{e}m/{r}^{2}& \left(17.1\right)\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}}& \left(17.2\right)\end{array}$
The net force on the mass at $r$ can be written as follows:
$\begin{array}{cc}F=m{\omega }^{2}r-m{g}_{0}{\left(\frac{{r}_{0}}{r}\right)}^{2}& \left(17.3\right)\end{array}$
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}}\left(22,000\phantom{\rule{6px}{0ex}}miles\right)& \left(17.4\right)\end{array}$
$\begin{array}{cc}F=m{g}_{0}{r}_{0}^{2}\left(\frac{r}{{r}_{s}^{3}}-\frac{1}{{r}^{2}}\right)& \left(17.5\right)\end{array}$
Figure 17.2: Radial forces acting on an orbiting mass.
Now consider a cable that extends from the earths surface $\left(r={r}_{0}\right)$ to a distance ${r}_{\ell }$ from the earth's surface, as shown in Figure 17.3. We need the cable to be in tension everywhere so that it doesn't buckle. If we get the tension at the earth's surface ($r={r}_{0}\right)$ we're in good shape. The mass increment for a cable of length $dr$ is $\rho Adr$, where A is the cross sectional area of the cable and $\rho$ is the density of the material from which it was made. We obtain the total force, ${F}_{0}$ at the earth's surface by integrating contributions to the force from the the whole length of the cable:
$\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& \left(17.6\right)\end{array}$
After integration we get:
$\begin{array}{cc}{F}_{0}=\rho A{g}_{0}{r}_{0}^{2}\left[\frac{1}{2{r}_{s}^{3}}\left({r}_{\ell }^{2}-{r}_{0}^{2}\right)+\left(\frac{1}{{r}_{\ell }}-\frac{1}{{r}_{0}}\right)\right]& \left(17.7\right)\end{array}$
Figure 17.4: ${F}_{0}$ as given by Eq. 17.7.
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}}$
Let's compare that to the best materials that are actually available. An Ashby plot of tensile strength (${\sigma }_{f}$) and density is shown in Figure 17.5. A line with a slope of 1 on this double logarithmic plot corresponds to a range of materials with a constant value of ${\sigma }_{f}/\rho$. The line drawn on Figure 17.5 corresponds to ${\sigma }_{f}/\rho$ is$\approx$2.8x10${}^{6}$ (in Si units), which is the most optimistic value possible for any known material - corresponding to the best attainable properties for diamond. Forgetting about any issues of cost, fracture toughness, etc., we could imagine that we see that we are a factor of 20 below the design requirement. So there's no way this is ever going to work, no matter how good your team of materials scientists is.
Figure 17.5: Ashby plot of tensile strength and density.
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)& \left(17.9\right)\end{array}$
If we assume $\sigma /\rho =2.8x1{0}^{6}\phantom{\rule{6px}{0ex}}\frac{N/{m}^{2}}{Kg/{m}^{3}}$, so that the system is operating at the value of $\sigma /f$ corresponding to the solid line in Figure 17.5 gives $\frac{{A}_{s}}{{A}_{0}}=110$. So in this case cable with a diameter of 1 cm at $r={r}_{0}$ needs to have a diameter of $\approx 10$ cm at $r={r}_{s}$. We don't have much leeway in decreasing $\sigma /\rho$, however. If the best we can do is $\sigma /\rho =1.0x1{0}^{6}\phantom{\rule{6px}{0ex}}\frac{N/{m}^{2}}{Kg/{m}^{3}}$, we get $\frac{{A}_{s}}{{A}_{0}}=1{0}^{21}$, which is clearly not workable. So we are stuck with the requirement that the cable have Good luck with that.

## 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 flat-ended cylindrical punch on a thick material?
• What happens to the normal stress distribution under a flat punch as an indented, low-modulus layer gets thinner and thinner? How is this related to the behavior of the compliance ($C\right),$ 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 load-displacement 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 Kelvin-Voigt 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, Nabarro-Herring 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, steady-state problems that are displacement-based, 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 steady-state FEA. The FEA process consists of the following steps:
1. formulation of element matrices,
2. their assembly into a structural matrix,
3. application of loads and boundary conditions,
4. solution of the structural equations,
5. 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 two-node 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 two-node 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 free-body diagrams,$\begin{array}{cc}A\sigma +{f}_{1}=0,\phantom{\rule{40px}{0ex}}-A\sigma +{f}_{2}=0.& \left(19.1\right)\end{array}$In turn, $\sigma$ is related to elastic modulus $E$ and axial strain $ϵ$,$\begin{array}{cc}\sigma =Eϵ,\phantom{\rule{40px}{0ex}}ϵ=\frac{{u}_{2}-{u}_{1}}{L}.& \left(19.2\right)\end{array}$Thus, we obtain$\begin{array}{cccc}\frac{AE}{L}\left({u}_{1}-{u}_{2}\right)& =& {f}_{1}& \left(19.3\right)\\ \frac{AE}{L}\left({u}_{2}-{u}_{1}\right)& =& {f}_{2}& \left(19.4\right)\end{array}$or, in matrix form,
$\begin{array}{cc}\left[k\right]\stackrel{\to }{d}=\stackrel{\to }{f}& \left(19.5\right)\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},& \left(19.6\right)\end{array}$and $\stackrel{\to }{d}$ and $\stackrel{\to }{f}$ are the element nodal displacement and the element nodal force vectors, respectively,$\begin{array}{cc}\stackrel{\to }{d}=\left\{\begin{array}{c}{u}_{1}\\ {u}_{2}\end{array}\right\},\phantom{\rule{40px}{0ex}}\stackrel{\to }{f}=\left\{\begin{array}{c}{f}_{1}\\ {f}_{2}\end{array}\right\}.& \left(19.7\right)\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.
Figure 19.1: A two-node bar element of length $L$ and cross-Section area $A$, with nodal forces ${f}_{1}$ and ${f}_{2}$, showing internal stress$\sigma$ and nodal d.o.f. ${u}_{1}$ and ${u}_{2}$.

#### 19.1.2 Structure of bar elements

Now consider a bar structure that consists of three segments attached end-to-end as shown in Figure 19.2(a). Each of these segments has length, cross-Sectional 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 cross-Sectional 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 $\stackrel{\to }{F}$ is applied to its right end. In a time-independent analysis, each segment of the structure can be modeled by a two-node bar element, resulting in a finite element mesh consisting of three elements and four nodes. For a known applied force $\stackrel{\to }{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 $\stackrel{\to }{D}$ and the nodal force vector $\stackrel{\to }{F}$ are:$\begin{array}{cc}\stackrel{\to }{D}=\left\{\begin{array}{c}{u}_{1}\\ {u}_{2}\\ {u}_{3}\\ {u}_{4}\end{array}\right\},\phantom{\rule{40px}{0ex}}\stackrel{\to }{F}=\left\{\begin{array}{c}{f}_{1}\\ {f}_{2}\\ {f}_{3}\\ {f}_{4}\end{array}\right\}.& \left(19.8\right)\end{array}$The stiffness equation for the overall structure can be written as$\begin{array}{cc}\stackrel{\to }{F}=\left[K\right]\stackrel{\to }{D}& \left(19.9\right)\end{array}$where $\left[K\right]$ is the $4×4$ global stiffness matrix.
Figure 19.2: (a) A bar consisting of three segments, with stiffnesses ${k}^{\left(1\right)}={E}^{\left(1\right)}{A}^{\left(1\right)}/{L}^{\left(1\right)}$, ${k}^{\left(2\right)}={E}^{\left(2\right)}{A}^{\left(2\right)}/{L}^{\left(2\right)}$ and ${k}^{\left(3\right)}={E}^{\left(3\right)}{A}^{\left(3\right)}/{L}^{\left(3\right)}$, respectively, has a force $F$ applied on its right end while the left end of the structure is fixed to a rigid wall. The structure is discretized with each segment represented by a two-node bar element. (b) The three elements and their element nodal displacements and forces.
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 two-node bar element then simply obeys the following equation:$\begin{array}{cc}{k}^{\left(e\right)}{d}^{\left(e\right)}={f}^{\left(e\right)},& \left(19.10\right)\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\},& \left(19.11\right)\\ \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\}.& \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\},& \left(19.12\right)\\ \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\}.& \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& \left(19.13\right)\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],& \left(19.14\right)\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\}.& \left(19.15\right)\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\}& \left(19.16\right)\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)}{ϵ}^{\left(e\right)}$. The element axial strain${ϵ}^{\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).& \left(19.17\right)\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\}.& \left(19.18\right)\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 7-bar 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}\stackrel{\to }{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}\\ \stackrel{\to }{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}& \left(19.19\right)\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×10$ matrix.
Figure 19.3: A plane truss loaded by a force$F$.

### 19.2 Preliminaries

As we've seen in the previous Section , the basic concept in FEM is the subdivision of the domain (structure) into non-overlapping 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 well-established tools of stress analysis, including stress-strain relations, strain-displacement relations, and energy considerations.

#### 19.2.1 Equilibrium equations

Moving beyond the simple free-body diagram approach of Figure 19.1, the more general forces and stresses acting on a differential 2D element$dx×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}& \left(19.20\right)\end{array}$
or, in matrix form,$\begin{array}{cc}{\partial }^{T}\sigma +F=0,& \left(19.21\right)\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\}.& \left(19.22\right)\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.
Figure 19.4: Stresses and body forces that act on a 2D differential element. The equilibrium equations (20) state that the differential element is in equilibrium.
The stresses and strain in an element are related through the constitutive matrix $E$ that contains the elastic constants,$\begin{array}{cc}\sigma =Eϵ.& \left(19.23\right)\end{array}$In two-dimensions, the strain vector$ϵ={\left[\begin{array}{ccc}{ϵ}_{x}& {ϵ}_{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].& \left(19.24\right)\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],& \left(19.25\right)\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 strain-displacement relation,$ϵ=\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{ˆ}{u}$ is$\begin{array}{cc}\stackrel{ˆ}{u}\left(x\right)={a}_{0}+{a}_{1}x,& \left(19.26\right)\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}& \left(19.27\right)\end{array}$
Solving for${a}_{0}$ and ${a}_{1}$ and rearranging, we obtain$\begin{array}{cc}\stackrel{ˆ}{u}\left(x\right)={N}_{1}\left(x\right){u}_{1}+{N}_{2}\left(x\right){u}_{2}=Nu,& \left(19.28\right)\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\},& \left(19.29\right)\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}},& \left(19.30\right)\end{array}$are the two linear shape functions.
Figure 19.5(a) shows a two-node 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}$
Figure 19.5: (a) Linear interpolation and shape functions. (b) Quadratic interpolation and shape functions.
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{ˆ}{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,& \left(19.31\right)\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\},& \left(19.32\right)\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)}& \left(19.33\right)\end{array}$Figure 19.5(b) shows a three-node 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 one-dimensional 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)$.
Figure 19.6: A linear triangle element.
The linear triangle uses linear interpolations for both $u$ and $v$, we will call the interpolating functions$\stackrel{ˆ}{u}$ and $\stackrel{ˆ}{v}$:
$\begin{array}{cc}\begin{array}{c}\stackrel{ˆ}{u}={N}_{1}{u}_{1}+{N}_{2}{u}_{2}+{N}_{3}{u}_{3}\\ \stackrel{ˆ}{v}={N}_{1}{v}_{1}+{N}_{2}{v}_{2}+{N}_{3}{v}_{3}\end{array}& \left(19.35\right)\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}& \left(19.36\right)\end{array}$
Notice that the same shape functions are used for both$\stackrel{ˆ}{u}$ and $\stackrel{ˆ}{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{ˆ}{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{ˆ}{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}& \left(19.37\right)\end{array}$
where the shape functions are given by
$\begin{array}{cc}\begin{array}{cc}{N}_{1}\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(a+x\right)\left(b-y\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(a-x\right)\left(b+y\right)}{4ab}\end{array}& \left(19.38\right)\end{array}$
Figure 19.7: A bilinear quadrilateral element.
##### 19.2.2.5 Number of Nodes and Order of Interpolation
As noted previously, displacements (or other degrees of freedom) are calculated at the nodes of the element. At any other point in the element, the displacements are obtained by interpolating the nodal displacements. In most cases, the order of interpolation is determined by the number of nodes in the element. Elements that have nodes only at their corners, such as the linear triangle, linear rectangle, and the 8-node brick element shown in Figure 19.8(a) use linear interpolation in each direction and are often called linear elements or first-order elements. Elements with midside nodes, and possibly interior nodes, such as the 20-node brick element shown in Figure 19.8(b) use quadratic interpolation and are often called quadratic elements or second-order elements.
Figure 19.8: (a) Linear and (b) quadratic 3D elements.

#### 19.2.3 Element Size

Imagine a one-dimensional 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 second-order 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.
Figure 19.9: The piecewise linear finite element solution converges to the exact solution as the mesh is refined. Note that the error, which is the difference between the exact and finite element solutions, is much smaller in the refined mesh of (b).

### 19.3 Example Problem

As an example of finite element analysis, consider the problem shown in Figure 19.10. The beam is of uniform cross Section , and the material is linearly elastic, isotropic and homogeneous. We wish to find the stresses and deflections due to the loading of the cantilever beam.
Figure