332: Mechanical Behavior of Materials

Department of Materials Science and EngineeringNorthwestern University

Table of Contents

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).

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, σ 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.
image: 8_home_ken_Mydocs_MSEcore_332_figures_propmaps50kbc.png image: 9_home_ken_Mydocs_MSEcore_332_figures_propmaps1900.png
image: 10_home_ken_Mydocs_MSEcore_332_figures_propmaps1945.png image: 11_home_ken_Mydocs_MSEcore_332_figures_propmaps2010.png
Figure 3.1: Available materials throughout history (from ref.[3]).

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 , σ 12 , etc.
  • Vectors are indicated with an arrow over the symbol, like P .
  • Unit vectors are indicated with a caret above the symbol, like n ˆ .
  • Matrices are enclosed in square brackets, like [ σ ]

4.1 Tensor Representation of Stress

image: 12_home_ken_Mydocs_MSEcore_332_figures_stress_state_2d.svg
Figure 4.1: 2-dimensional stress tensor
The stress applied to an object, which we denote as σ ij or [ σ ] is the force acting over an area of an object, divided by the area over which this force is acting. Note note that [ σ ] is a matrix with individual components, σ ij specified by the indices i and j . These indices have the following significance:
To obtain the Engineering stress , [ σ eng ] , we use the undeformed areas of the stress-free object to obtain the stress tensor, whereas the true stress (which is what we generally mean when we write [ σ ] ) we use the actual areas in the as-stressed state.
The stress matrix is a tensor , which means that it obeys the coordinate transformation laws describe below. In two dimensions it has the following form:
[ σ ] =[ σ xx σ xy σ yx σ yy ] (4.1)
The stress tensor must be symmetric, with σ xy = σ 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:
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 , P 2 and P 3 ) acting on a given plane. The plane is specified by the orientation of the unit vector, 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. ( n 1 2 + n 2 2 + n 3 2 ) 1/2 =1. The relationship between P , σ and n is as follows:
[ P 1 P 2 P 3 ]=A[ σ 11 σ 12 σ 13 σ 12 σ 22 σ 23 σ 13 σ 23 σ 33 ][ n 1 n 2 n 3 ] (4.2)
or in more compact matrix notation:
P =A[ σ ] n ˆ (4.3)
Here A is the total cross sectional area of the plane that we are interested in. (If you need a refresher on matrix multiplication, the Wikipedia page on Matrix Multiplication (https://en.wikipedia.org/wiki/Matrix_multiplication)[16] is very helpful).
image: 13_home_ken_Mydocs_MSEcore_332_figures_stress_tensor.svg
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:
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 n =( 1,0,0 ) in Eq. 4.2, we get the following for the stresses acting on the 1 face of the volume element:
P 1 /A= σ 11 P 2 /A= σ 12 P 3 /A= σ 13 (4.4) Equivalent expressions can be obtained for the stresses acting on the 2 and 3 faces, by setting n =( 0,1,0 ) and n =( 0,0,1 ) , 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 2 and 3 . As an example, consider the simple counterclockwise rotation around the 3 axis by an angle φ , 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 θ ij , where i specifies the transformed axes ( 1 , 2 or 3 ) and j specifies the untransformed axis (1, 2 or 3). In our simple example, the angle between the 1 and 1 axes is φ , so θ 11 = φ . The angle between the 2 and 2 axes is also φ , so θ 22 = φ . The 3 axis remains unchanged in our rotation example, so θ 33 =0. The 3/ 3 axis remains perpendicular to the 1, 1 , 2, 2 axes, so we have θ 31 = θ 32 = θ 13 = θ 23 =9 0 . Finally, we see that the angle between the 1 and the 2 axis is 90- φ ( θ 12 =90- φ ) and the angle between the 2 and 1 axis is 90+ φ ( θ 21 =90+ φ ). The full [ θ ] matrix in this case is as follows: [ θ ] =[ φ 90- φ 90 90+ φ φ 90 90 90 0 ] (4.5)
Note that the [ θ ] matrix is NOT symmetric ( θ ij θ ji ), so you always need to make sure the first index, i, (denoting the row in the [ θ ] matrix) corresponds to the transformed axes, and the second index, j (denoting the column in the [ θ ] matrix) corresponds to the original, untransformed axes.
image: 14_home_ken_Mydocs_MSEcore_332_figures_Stress_transformation_2D.svg
Figure 4.3: Rotation of coordinate system. The coordinate system is rotated by θ about the 3 axis, transforming the 1 axis to 1 and the 2 axis to 2 .

4.2.2 Expressions for the Stress Components

Once we specify all the different components of [ θ ] , 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:
σ ij = k,l cos θ jk cos θ il σ kl (4.6)
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, σ 12 is given as follows:
σ 12 = cos θ 21 cos θ 11 σ 11 + cos θ 21 cos θ 12 σ 12 + cos θ 21 cos θ 13 σ 13 + cos θ 22 cos θ 11 σ 21 + cos θ 22 cos θ 12 σ 22 + cos θ 22 cos θ 13 σ 23 + cos θ 23 cos θ 11 σ 31 + cos θ 23 cos θ 12 σ 32 + cos θ 23 cos θ 13 σ 33 (4.7)
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:
[ σ ] =[ 5x1 0 6 0 0 0 0 0 0 0 0 ] (4.8)
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 φ =4 5 . The following MATLAB code solves for the full transformed tensor, with σ ij given by Eq. 4.8 and [ θ ij ] given by Eq. 4.5 with φ =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
for i=1:3
    for j=1:3
        for k=1:3
            for l=1:3
            end  % cosd takes degrees instead of radians as the unit
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]]
theta = np.deg2rad(theta)  # trig functions need angles in radians
for i in [0, 1, 2]:
    for j in [0, 1, 2]:
        for k in [0, 1, 2]:
            for l in [0, 1, 2]:

print(sigp)  # display the transformed tensor components
Note that Python is not at all necessary for anything in 332. It is, however, rapidly replacing MATLAB and will almost certainly be more useful to you than MATLAB once you graduate. It is free, powerful, and quite easy to learn if you already know MATLAB. Python examples of the MATLAB codes are include in Section 20. These are presented as examples of how to do some useful things in Python.
image: 15_home_ken_Mydocs_MSEcore_332_figures_rotate45-result.svg
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:
[ σ ] =[ 2.5x1 0 6 -2.5x1 0 6 0 -2.5x1 0 6 2.5x1 0 6 0 0 0 0 ] (4.9)
Note the following:

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:
[ Q ] = cos [ θ ] (4.10)
For the simple case of rotation about the z axis, the angles are given by Eq. 4.5, so that [ Q ] is given as:
[ Q ] =[ cos φ cos ( 90- φ ) cos 90 cos ( 90+ φ ) cos φ cos 90 cos 90 cos 90 cos 0 ]=[ cos φ sin φ 0 - sin φ cos φ 0 0 0 1 ] (4.11)
The transformed stress is now obtained by the following simple matrix multiplication:
[ σ ] =[ Q ] [ σ ] [ Q ] T (4.12)
where the [ Q ] T is the transpose of [ Q ] :
Q T ( i,j ) =Q( j,i ) (4.13)
For the rotation by φ around the z axis, [ Q ] T is given by the following:
[ Q ] T =[ cos φ - sin φ 0 sin φ cos φ 0 0 0 1 ] (4.14)
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:
sig=zeros(3); % create stress tensor and set to zero
sig(1,1)=5e6; % this is the only nonzero component

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

Any stress state (true stress) can be written in terms of three principal stresses σ 1 p , σ 2 p and σ 3 p , applied in three perpendicular directions as illustrated in Figure 4.5.
image: 16_home_ken_Mydocs_MSEcore_332_figures_principal_stresses.svg
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:
[ σ ] =[ σ 1 p 0 0 0 σ 2 p 0 0 0 σ 3 p ] (4.15)
In order to gain some insight into the three points mentioned above, it is useful to consider a range of rotation angles, and not just a singe rotation angle of 4 5 . One way to do this is to use the symbolic math capability of MATLAB to obtain the full stress tensor as a function of the rotation angle. We'll use a slightly more complicated stress tensor as our input, which has principal stresses in all three untransformed directions, but without shear stresses, as shown below in Eq. 4.15.
The following MATLAB script results in the output shown in Figure 4.6.
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
The components of the stress tensor in the rotated coordinate system are:
[ σ ] =[ σ 1 p - σ 1 p sin 2 φ + σ 2 p sin 2 φ sin ( 2 φ ) ( σ 1 p - σ 2 p ) /2 0 sin ( 2 φ ) ( σ 1 p - σ 2 p ) /2 σ 2 p + σ 1 p sin 2 φ - σ 2 p sin 2 φ 0 0 0 0 ] (4.16)
image: 17_home_ken_Mydocs_MSEcore_332_figures_matlab_symbolic.png
Figure 4.6: MATLAB output generated by rotatesym_matrix.m.
This is not yet a very illuminating result, but it is the basis for the Mohr circle construction , which provides a very useful way to visualize two dimensional stress states. This construction is described in more detail in the following Section .

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, σ 1 p and σ 2 p , and by the orientation of the principal axes. The Mohr circle is drawn with a diameter of σ 1 p - σ 2 p , centered at ( σ 1 p + σ 2 p ) . We can use MATLAB's symbolic math capability to derive the Mohr circle construction as follows:
% Mohr's circle derivation - beginning from previous calculation of sigp
center=(sig1p+sig2p)/2;  % the center of Mohr's circle
radius=(sig1p-sig2p)/2;  % radius of Mohr's circle
% now we shift the origin to the center of Mohr's circle, and 
% normalize by the radius of Mohrs circle
% shear stresses are plotted on the y axis, so we don't need to shift the
% origin, we just need to normlized by the radius

The output of Mohr_circle.m is shown In Figure 4.7.
image: 18_home_ken_Mydocs_MSEcore_332_figures_mohr_circle_output.png
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 φ 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).
image: 19_home_ken_Mydocs_MSEcore_332_figures_Mohr_Circle.svg
Figure 4.8: Mohr's circle construction.
In general, there are three principal stresses, σ 1 p , σ 2 p and σ 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 σ 1 p is the largest principal stress and that σ 3 p is the smallest principal stress, i.e. σ 1 p > σ 2 p > σ 3 p . An important result is that the largest shear stress, τ max , is given by the difference between the largest principal stress and the smallest one:
τ max = 1 2 ( σ 1 p - σ 3 p ) (4.17)
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.
image: 20_home_ken_Mydocs_MSEcore_332_figures_Mohr_Circle3d.svg
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 σ 1 p is directed, C 2 corresponds to rotation around the direction in which σ 2 p is directed, and C 3 corresponds to rotation around the direction in which σ 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:
[ σ ] =[ 3 0 2 0 3 0 2 0 5 ] MPa (4.18)
In this case there is only one non-zero shear stress ( σ 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 σ 13 , one of the principal stresses is σ 22 =3 MPa.
  2. Draw a Mohr's circle construction using the two normal stresses and the non-zero shear stress, in this case σ 11 , σ 33 and σ 13 . Mohr's circle is centered at the the average of these two normal stresses, in our case at ( σ 11 + σ 33 ) /2 = 4 MPa. We refer to this stress as σ c , as illustrated below in Figure 4.10.
  3. Determine the radius of the circle, R M , is given as:
    R m = ( σ 33 - σ c ) 2 + σ 13 2 = ( 5-4 ) 2 + σ 13 2 =2.24 MPa
  4. The principal stresses are given by the intersections of the circle with the horizontal axis:
    σ p = σ c ± R m =6.24 MPa, 1.76 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.
image: 21_home_ken_Mydocs_MSEcore_332_figures_Mohr_Circle_radius_and_center.svg
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, τ rss , for sample in uniaxial tension is given by the following expression:
τ rss = σ 1 p cos φ cos λ (4.19)
where σ 1 p is the applied tensile stress, φ is the angle between the tensile axis and a vector normal to the plane of interest, and λ is the angle between the tensile axis and the direction of the shear stress. This shear stress has to be in the plane itself, so for a 2-dimensional sample λ + φ =9 0 . This means we can rewrite Eq. 4.19 in the following way:
τ rss = σ 1 p cos φ cos ( 9 0 - φ ) (4.20)
We can use the identities cos ( 90- φ ) = sin φ and sin ( 2 φ ) =2 sin φ cos φ to obtain the following:
τ rss = σ 1 p 2 sin ( 2 φ ) (4.21)
We can get the same thing from the Mohr's circle construction to redefine the axes by a rotation of φ . The shear stress is simply the radius of the circle ( φ 1 p /2 in this case) multiplied by sin ( 2 φ ) . Mohr's circle also gives us the normal stresses:
σ 11 = σ 1 p 2 + σ 1 p 2 cos ( 2 φ ) σ 22 = σ 1 p 2 - σ 1 p 2 cos ( 2 φ ) (4.22)
The untransformed 2-dimensional stress tensor looks like this:
[ σ ] =[ σ 1 p 0 0 0 ] (4.23)
The transformed stress tensor (after rotation by φ to give the resolved shear stress) looks like this:
[ σ ] =[ σ 1 p 2 + σ 1 p 2 cos ( 2 φ ) σ 1 p 2 sin ( 2 φ ) σ 1 p 2 sin ( 2 φ ) σ 1 p 2 - σ 1 p 2 cos ( 2 φ ) ] (4.24)
image: 22_home_ken_Mydocs_MSEcore_332_figures_Mohr_Circle_uniaxial.svg
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

Principal stresses can easily by calculated for any stress state just by obtaining the eigenvalues of the stress tensor. In addition, the orientation of the principal axes (the coordinate system for which there are no off-diagonal components in the stress tensor). If you need a refresher on what eigenvalues and eigenvectors actually are, take a look at the appropriate Wikipedia page (http://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors). We'll use MATLAB to do this, using the 'eigs' command .
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
% 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
Here's the output generated by this script:
image: 23_home_ken_Mydocs_MSEcore_332_figures_principal_stress_result.png
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:

4.3.4 Stress Invariants

Some quantities are invariant to choice of axes. The most important one is the
hydrostatic pressure
, p , given by summing the diagonal components of the stress tensor and dividing by 3:
p=- 1 3 ( σ 11 + σ 22 + σ 33 ) =- 1 3 ( σ 1 p + σ 2 p + σ 3 p ) (4.25)
The negative sign appears because a positive pressure is compressive, but positive stresses are tensile. The hydrostatic pressure is closely related to a quantity referred to as the 'first stress invariant', I 1 :
I 1 = σ 11 + σ 22 + σ 33 (4.26)
The second and third stress invariants, I 2 and I 3 , are also independent of the way the axes are defined:
I 2 = σ 11 σ 22 + σ 22 σ 33 + σ 33 σ 11 - σ 12 2 - σ 13 2 - σ 23 2 (4.27)
I 3 = σ 11 σ 22 σ 33 - σ 11 σ 23 2 - σ 22 σ 13 2 - σ 33 σ 12 2 +2 σ 12 σ 13 σ 23 (4.28)
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')
Each of these definitions of strain describe the way different points an object move relative to one another when the material is deformed. Consider two points P 1 and P 2 , separated initially by the increments d x , d y and d z along the x , y and z directions. After the deformation is applied, these points move by the following amounts, as illustrated in Figure 4.13:
image: 24_home_ken_Mydocs_MSEcore_332_figures_straindef.svg
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:
du= u x dx+ u y dy+ u z dz dv= v x dx+ v y dy+ v z dz dw= w x dx+ w y dy+ w z dz (4.29)
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:
e xx = u x e yy = v y e zz = w z (4.30)
The engineering shear strains are defined as follows:
e yz = w y + v z e xz = u z + w x e xy = v x + u y (4.31)
Note: shear strains are often represented by the lower case Greek gamma to distinguish them from normal strains:
γ xy e xy ; γ yz e yz ; γ xz = e xz (4.32)

4.4.2 Tensor Shear Strains

Engineering strains relate two vectors to one another, ( dx,dy,dz ) and ( du,dv,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 ε to indicate tensor strains. The tensor normal strains are exactly the same as the engineering normal strains:
ε xx = e xx ε yy = e yy ε zz = e zz (4.33)
Engineering shear strains ( e yz , e xz , e zy ) are divided by two to give tensor shear strains:
ε yz = e yz /2 ε xz = e xz /2 ε xy = e xy /2 (4.34)
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

We can also define the strain by considering a cube of side that is deformed into a parallelepiped with dimensions of (along principal strain axes). After deformation, the cube has dimensions of λ 1 , λ 2 , λ 3 . Alternatively, a sphere of radius r 0 is deformed into and ellipsoid with principal axes of λ 1 r 0 , λ 2 r 0 and λ 3 r 0 , as shown in Figure 4.14. The quantities λ 1 , λ 2 , λ 3 are extension ratios , and are related to the principal strains as follows:
λ 1 =1+ e 1 P λ 2 =1+ e 2 P λ 3 =1+ e 3 P (4.35)
The true strains , e 1 t , are obtained as by taking the natural log of the relevant extension ratio. For example, for a uniaxial tensile test, the true strain in the tensile direction (assumed to be the 1 direction here) is:
e 1 t = ln ( λ 1 ) (4.36)
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/ , where is the current length of the material as it is being deformed. If the initial length is 0 and the final, deformed length is f , then the total true strain, e 1 true is obtained by integrating the incremental strains accumulated throughout the entire deformation history:
e 1 t = 0 f d = . ln | 0 f = ln ( f 0 ) = ln λ 1 (4.37)
image: 25_home_ken_Mydocs_MSEcore_332_figures_strain_ellipsoid.svg
Figure 4.14: Unit sphere deformed into a strain ellipsoid with dimensions of λ 1 , λ 2 , λ 3 .
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.
A more detailed description of generalized strain, with a lot of relevant matrix math, is provided in the Wikipedia article on finite strain theory (https://en.wikipedia.org/wiki/Finite_strain_theory). This information is not needed for this class, but the reference is provided here as a useful resource for large strain. If you come across concepts like the
Cauchy-Green deformation tensor
or the
Finger deformation tensor
, this article provides a useful introduction (but prepared for a lot of matrix math). These concepts appear in a range of useful description of mechanical response, including many in the biomedical field (muscle actuation, deformation of skin, etc).

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

Simple shear is a two dimensional strain state, which means that one of the principal strains is zero (or one of the principal extension ratios is 1).
image: 26_home_ken_Mydocs_MSEcore_332_figures_shear_strain.svg
Figure 4.15: Schematic representation of simple shear.
The stress tensor looks like this:
σ =[ 0 σ xy 0 σ xy 0 0 0 0 0 ] (4.38)
From the definition of the engineering shear strain (Eq. 4.31) we have:
e xy = u d (4.39)
We need to divide the engineering shear strains by 2 to get the tensor strains, so we get the following:
ε =[ 0 e xy /2 0 e xy /2 0 0 0 0 0 ] (4.40)
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.
G ( shearmodulus ) = σ xy e xy (4.41)
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 = γ ). One thing that we notice from this plot is that γ is simply given by the difference between the two principal strains:
γ = e 1 p - e 2 p = λ 1 - λ 2 (4.42)
For simple shear this relationship is valid, even for large strains, even though there are other aspects of the Mohr's circle construction that no longer work at large strains. The primary difficulty is that the frame of reference for the strained and unstrained case are not the same. In general, strains rotate the frame of reference by an amount that we don't want to worry about for the purposes of this course. For small strains typically obtained in metals or ceramics (with strain amplitudes of a few percent or less) we don't need to worry about this rotation, but it can become important for polymeric systems that undergo very large strains prior to failure. However, if all we want a measure of the shear strain in the material, we can still use Eq. 4.42, regardless of how large the principal strains (and corresponding extension ratios) actually are.
image: 27_home_ken_Mydocs_MSEcore_332_figures_Mohr_Circle_strain.svg
Figure 4.16: Mohr's circle construction for SMALL strains.

4.5.3 Hydrostatic Compression

bulk compressive modulus
, K b , of a material describes its resistance to a change in density. Formally, it is defined in terms of the dependence of the volume of the material on the hydrostatic pressure, p :
K b =-V dp dV (4.43)
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:
σ =[ -p 0 0 0 -p 0 0 0 -p ] (4.44)
where p is the
hydrostatic pressure
image: 28_home_ken_Mydocs_MSEcore_332_figures_hydrostatic.svg
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:
σ =[ 0 0 0 0 0 0 0 0 σ 33 ] (4.45)
image: 29_home_ken_Mydocs_MSEcore_332_figures_uniaxial_extension.svg
Figure 4.18: Uniaxial tensile deformation.
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:
e 33 = Δ 0 (4.46)
e 22 = Δ w w (4.47)
From these strains we can define
Young's modulus
, E , and
Poisson's ratio
, ν : E= σ 33 / e 33 (4.48)
ν =- e 22 / e 33 (4.49)

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=[ 0 0 0 0 0 0 0 0 e 33 ]
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 σ 11 and σ 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 which is the ratio of σ 33 to e 33 for this strain state. Note that this deformation state changes both the shape and volume of the material, so E involves both G and K b :
E = σ 33 e 33 = K b + 4 3 G (4.50)
image: 30_home_ken_Mydocs_MSEcore_332_figures_longitudinal.svg
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.
Table 1: Representative elastic moduli for different materials.
G (Pa)
K b (Pa)
1.0x10 5
2.2x10 9
1 0 4
2.2x10 9
1 0 9
2x10 9
8x10 10
1.6x10 11

4.7 Case Study: Speed of Sound

The speed of sound, or sound velocity , V sound , is actually a mechanical property. It is related to a modulus, M , in the following way:
V sound = M ρ (4.51)
Here ρ 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:
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:
P= n V RT (4.52)
If the compression is applied slowly enough so that the temperature of the gas can equilibrate, we have:
K b =-V dP dV = n V RT=P (4.53)
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 ρ =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: x1 , y2 , z3 . 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.
Table 2: Definition of the matrix components of stress and strain.
Engineering stress
Matrix Stress
Engineering Strain
Matrix Strain
σ 11
σ 1
e 11
e 1
σ 22
σ 2
e 22
e 2
σ 33
σ 3
e 33
e 3
σ 23
σ 4
e 23
e 4
σ 13
σ 5
e 13
e 5
σ 12
σ 6
e 12
e 6
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

[ e 1 e 2 e 3 e 4 e 5 e 6 ]=[ 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 ][ σ 1 σ 2 σ 3 σ 4 σ 5 σ 6 ] (5.1)
The matrix must be symmetric, with s ij = s ji , so there are a maximum of 21 independent compliance coefficients:
[ e 1 e 2 e 3 e 4 e 5 e 6 ]=[ 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 ][ σ 1 σ 2 σ 3 σ 4 σ 5 σ 6 ] (5.2)
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. [ σ 1 σ 2 σ 3 σ 4 σ 5 σ 6 ]=[ 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 ][ e 1 e 2 e 3 e 4 e 5 e 6 ] (5.3)

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.
image: 31_home_ken_Mydocs_MSEcore_332_figures_orthorhombic_symmetry.svg
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 ( σ 1 , σ 2 or σ 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 e 2 e 3 e 4 e 5 e 6 ]=[ 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 ][ σ 1 σ 2 σ 3 σ 4 σ 5 σ 6 ] (5.4)

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 :
[ e 1 e 2 e 3 e 4 e 5 e 6 ]=[ s 11 s 12 s 13 0 0 0 s 12 s 11 s 13 0 0 0 s 13 s 13 s 33 0 0 0 0 0 0 s 44 0 0 0 0 0 0 s 44 0 0 0 0 0 0 2( s 11 - s 12 ) ][ σ 1 σ 2 σ 3 σ 4 σ 5 σ 6 ] (5.5)
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. 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: e 3 = Δ /; e 1 = e 2 = Δ d/d (5.6)
We then obtain s 33 and s 13 from Eq. 5.5, recalling that σ 3 is the only non-zero stress component in this situation:
s 33 = e 33 / σ 3 s 13 = e 1 / σ 3 (5.7)
image: 32_home_ken_Mydocs_MSEcore_332_figures_fiber_extension.svg
Figure 5.2: Fiber extension. Fiber Torsion: Measurement of s 44
We obtain the shear modulus by looking at the torsional stiffness of the fiber, i.e. , the torque, T , required to rotate the top and bottom of the fiber by an angle θ 0 , as illustrated in Figure 5.3
image: 33_home_ken_Mydocs_MSEcore_332_figures_fiber_torsion.svg 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 θ around the z axis. The shear strain in the θ -z plane depends only on r , and is given by :
e θ z =r d θ dz =r θ 0 (5.8)
The corresponding shear stress is obtained by multiplying by the shear modulus, G characterizing deformation in the 1-2 and 2-3 planes:
σ θ z =Gr θ 0 / (5.9)
We integrate the shear stress to give the torque, T :
T= 0 d/2 r σ θ z 2 π r r = π G θ 0 d 4 32 (5.10)
So once we know the torsional stiffness of the fiber (the relationship between the applied T and θ 0 ) we know the shear modulus, G . This shear modulus is simply the inverse of s 44 :
G= 1 s 44 (5.11)

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:
b 2 = 2F d 0 s 11 (5.12)
where P is the force applied to the fiber, d 0 is its undeformed diameter and is its length.
A practical situation that is often observed is that friction between the fiber and confining surfaces keeps the fiber length from increasing, so the strain in the three direction must be zero. If we express the stress state in terms of the principle stresses in 1 and 2 directions, we have the following from Eq. 5.5:
e 3 = s 13 ( σ 1 p + σ 2 p ) + s 33 σ 3 (5.13)
Setting e 3 to zero in this equation gives the following for σ 3 :
σ 3 = - s 13 s 33 ( σ 1 p + σ 2 p ) (5.14)
A consequence of this stress is that the frictionless expression for b gets modified to the following:
b 2 = 2F d 0 ( s 11 - s 13 2 s 33 ) (5.15)
image: 34_home_ken_Mydocs_MSEcore_332_figures_fiber_transverse_compression.svg
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.
d d 0 =f( P, s 11 , s 13 , s 33 , s 12 ) (5.16)

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:
[ e 1 e 2 e 3 e 4 e 5 e 6 ]=[ s 11 s 12 s 12 0 0 0 s 12 s 11 s 12 0 0 0 s 12 s 12 s 11 0 0 0 0 0 0 s 44 0 0 0 0 0 0 s 44 0 0 0 0 0 0 s 44 ][ σ 1 σ 2 σ 3 σ 4 σ 5 σ 6 ] (5.17)

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:
[ e 1 e 2 e 3 e 4 e 5 e 6 ]=[ s 11 s 12 s 12 0 0 0 s 12 s 11 s 12 0 0 0 s 12 s 12 s 11 0 0 0 0 0 0 2( s 11 - s 12 ) 0 0 0 0 0 0 2( s 11 - s 12 ) 0 0 0 0 0 0 2( s 11 - s 12 ) ][ σ 1 σ 2 σ 3 σ 4 σ 5 σ 6 ] (5.18)
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( s 11 - s 12 ) , so there are two independent elastic constants. The shear modulus, G , Young's modulus E and Poisson's ratio, ν are given as follows:
G=1/2( s 11 - s 12 ) ;E=1/ s 11 ; ν =- s 12 / s 11 (5.19) 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:
K b =-V dP dV (5.20)
The stress state in this case is as follows:
σ =[ -p 0 0 0 -p 0 0 0 -p ] (5.21)
From the compliance matrix (Eq. 5.18) we get e 1 = e 2 = e 3 =-p( s 11 +2 s 12 ) . The change in volume, Δ V can be written in terms of the three principal extension ratios, λ 1 , λ 2 and λ 3 :
Δ V V 0 = V V 0 -1= λ 1 λ 2 λ 3 -1=( 1+ e 1 ) ( 1+ e 2 ) ( 1+ e 3 ) -1 e 1 + e 2 + e 3 (5.22)
Now we use the fact that for small x , ( 1+x ) 3 1+3x , so we have:
Δ V V 0 = e 1 + e 2 + e 3 =-3p( s 11 +2 s 12 ) (5.23)
Recognizing that the derivative dP/dV in the definition of K b can be written as the limit of p/ Δ V for very small p allows us to obtain the expression we want for K b : K b = lim p0 -p Δ V/ V 0 = 1 3( s 11 +2 s 12 ) (5.24) Relationship between the Isotropic Elastic Constants:
Because there are only two independent elastic constants for an isotropic system E and ν can be expressed in terms of some combination of K b and G . For E the relevant relationship is as follows.
E= 9G 3+G/ K b =2G( 1+ ν ) (5.25)
We can also equate the two expressions for G in Eq. 5.25 to get the following expression for ν :
ν = 3-2G/ K b 6+2G/ K b (5.26)
Note that if K b G , E3G and ν 0.5 .

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.
Field (row)
Response (column)
tensor properties of rank 0 (scalars)
specific heat
C (1x1)
tensor properties of rank 1 (vectors)
p i (3x1)
D i (electric displacement)
tensor properties of rank 2
Thermal expansion
α i (6x1)
e i
Dielectric permittivity
κ ij (3x3)
E j (elec. field)
D i (electric displacement)
Electrical conductivity
σ ij (3x3)
E j
j i (current density)
Σ ij (3x3)
T x j (Temp. gradient)
E i (electric field)
tensor properties of rank 3
d ij (3x6)
σ j (stress)
D i (electric displacement)
Converse piezeoelectricty
d ij (6x3)
E j (elec. field)
e i (strain)
tensor properties of rank 4
s ij (6x6)
σ j
e i
Table 3: Some symmetry-related properties.
image: 35_home_ken_Mydocs_MSEcore_332_figures_symmetry_props_isotropic.svg
cubic (point group 432)
image: 36_home_ken_Mydocs_MSEcore_332_figures_symmetry_props_432.svg
quartz (point group 32)
image: 37_home_ken_Mydocs_MSEcore_332_figures_symmetry_props_32.svg
Figure 6.1: Symmetry-property maps for 3 of the 32 point groups.

Contact Mechanics

In a simple tensile test involving a sample with a uniform cross section, the stresses and strains are both uniform throughout the entire sample. In almost any real application where we care about mechanical properties, this is not the case however. A simple example of this is the case where we press a rigid, cylinder into s soft, compliant material as shown in Figure 6.1.
image: 38_home_ken_Mydocs_MSEcore_332_figures_flatpunch.svg
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:
In order not to get too hung up in issues related to the sign, we define P t and δ t as the tensile loads and displacements:
P t =-P δ t =- δ (7.1)

7.2 Flat Punch Indentation

(Note: Many of issues presented here are discussed in more detail in a published review article: see ref.[13]).
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, δ , applied to the indenter.
image: 39_home_ken_Mydocs_MSEcore_332_figures_flat_punch_geometry.svg
Figure 7.2: Flat punch contact geometry. For an elastic half space, h.

7.2.1 Flat Punch: Approximate Result for an elastic half space.

For an elastic half space ( h ), 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, δ from the following approximate concepts:
e avg =- δ /a (7.2)
C 0 δ P 1 π Ea (7.4)

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 ν 1 and ν 2 , then the expression for C 0 is:
C 0 δ P = δ t P t = 1 2 E r a (7.5) where E r is the following reduced modulus :
1 E r = 1- ν 1 2 E 1 + 1- ν 2 2 E 2 (7.6)
Note that for a stiff indenter, ( E 2 E 1 ) we have E r = E 1 1- ν 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;
[ e 1 e 2 e 3 e 4 e 5 e 6 ]= 1 E [ 1 - ν - ν 0 0 0 - ν 1 - ν 0 0 0 - ν - ν 1 0 0 0 0 0 0 2( 1+ ν ) 0 0 0 0 0 0 2( 1+ ν ) 0 0 0 0 0 0 2( 1+ ν ) ][ σ 1 σ 2 σ 3 σ 4 σ 5 σ 6 ] (7.7)
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 =- ν /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 σ 3 =0 , we have:
e 2 = - ν E σ 1 + σ 2 E (7.8)
If e 2 is constrained to be zero, we have:
σ 2 = ν σ 1 (7.9)
Now we can put this value back into Eq. 7.7 and solve for e 1 :
e 1 = 1 E ( σ 1 - ν 2 σ 1 ) (7.10)
The plane strain modulus relates σ 1 to e 1 , which for the case assumed above ( e 2 =0 ) gives:
E r = σ 1 e 1 = E 1- ν 2 (7.11)

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= ) is given by the following expression[6]:
u z =( 2 δ t / π ) arcsin ( a/r ) (7.12) where δ t is the applied tensile displacement and a is the actual radius of the contact area between the punch and the substrate. In Figure 7.3 we compare the shapes of the surface for the following two cases:
image: 40_home_ken_Mydocs_MSEcore_332_figures_flatpunchtensile.svg
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:
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:
G = d d A c ( W- U el ) (7.13)
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:
G = G c (7.14)
The lowest possible value of G c is 2 γ where γ 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).
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, δ 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:
C= . d δ t d P t | A c (7.15)
image: 41_home_ken_Mydocs_MSEcore_332_figures_crack_extension_G_derivation.svg
image: 42_home_ken_Mydocs_MSEcore_332_figures_g_derivation.svg
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 δ 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+ δ 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:
δ W= 1 2 P t 2 δ C (7.16)
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 δ A c becomes very small to replace δ W/ δ C with the derivative, dW/dC to obtain:
G = lim δ A c 0 δ W δ A c = P t 2 2 dC d A c (7.17)

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 ( d G /dA>0) and unstable contact ( d G /dA<0) . 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:
image: 43_home_ken_Mydocs_MSEcore_332_figures_stable_unstble_crack_.svg
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 d G /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 d G /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:
G =- P t 2 2 dC dA =- P t 2 4 π a dC da (7.18)
where we have assumed that the contact area remains circular, with A= π 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:
dC dA =- dC d A c (7.19)
With C= C 0 =1/2 E r a (Eq. 7.5) we obtain the following expression for G :
G = P t 2 8 π E r a 3 (7.20)
In some situations it is more convenient to express the energy release rate in terms of the tensile displacement, δ t . The most general expression is used by using C= δ t / P t to substitute P t with δ t /C in Eq. 7.17:
G =- δ t 2 2 C 2 dC dA =- δ t 2 4 π a C 2 dC dA (7.21)
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:
G = E r δ t 2 2 π a (7.22)
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= π 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, σ c :
σ c = P c π a 2 = ( 8 E r G c π a ) 1/2 (7.25)
Note that the detachment stress increases with decreasing punch size.
image: 44_home_ken_Mydocs_MSEcore_332_figures_geckosatae.jpg
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:
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 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.
image: 45_home_ken_Mydocs_MSEcore_332_figures_manypillars.svg
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:
C= 1 2 E r a f C (7.26)
image: 39_home_ken_Mydocs_MSEcore_332_figures_flat_punch_geometry.svg
Figure 7.8: A thin, compliant layer being indented with a rigid, cylindrical punch.
For an elastic half space ( h) f C =1. The factor f C accounts for changes in the compliance due to the decreased thickness of the layer. In general it depends on Poisson's ratio for the compliant layer and the confinement ratio, a/h (the ratio of the punch radius to the thickness of the layer). For an incompressible compliant layer with ν =0.5 the following expression for f C provides an excellent approximation to the behavior of the compliance on the aspect ratio, a/h :[13]
f C = [ 1+1.33( a/h ) +1.33 ( a/h ) 3 ] -1 (7.27)
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:
G =- P t 2 4 π a dC da = P t 2 8 π E r a 3 f G p (7.28)
Here f G p 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 :
f G p = dC da / . dC da | a/h=0 (7.29)
image: 46_home_ken_Mydocs_MSEcore_332_figures_flat_punch_finite_size_corrections.svg
Figure 7.9: Geometric correction factors for the flat punch geometry (generated with python code given as example 3. in the appendix).
Finally, we can use the the fact that P t = δ t /C to get a similar expression for G in terms of the tensile displacement, δ t :
G = E r δ t 2 2 π a f G δ (7.30)
In this case f G δ 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 ( δ t ) 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 G p , and the a/h dependence of C is accounted for by f C , so we obtain the following for f G δ :
f G δ = f G p f c 2
The confinement functions f G p and f G δ are both equal to one for a/h=0 and are plotted as a function of a/h for ν =0.5 in Figure 7.9. A practical consequence of the decrease in f G p 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 G p , a larger tensile load needs to be applied in order for G to exceed the critical energy release rate, G c .

7.4 Contact of Paraboloids

7.4.1 Non-Adhesive Case

Suppose that the indenter is not flat, but has a parabolic profile that can be described by the following expression:
z= A p r 2 (7.31)
Here z is the vertical distance from the apex of the parabola, r is the radial distance from symmetry axis for the paraboloid and A p is a constant that defines the shape of the paraboloid. A sphere has a parabolic shape near the apex, which can be seen by considering the equation for a sphere of Radius R that has it's center at r=0 , z=R (see Figure 7.10):
r 2 + ( z-R ) 2 = R 2 (7.32)
Solving Eq. 7.32 for z gives:
z=R( 1 ± 1- ( r/R ) 2 ) (7.33)
For small x, 1-x 1-x/2 , so for rR we have:
z= r 2 2R (7.34) 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.

image: 47_home_ken_Mydocs_MSEcore_332_figures_hertzschematic.svg
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 ( δ h in Figure 7.10) is given by the following expression:
δ h = a 2 /R (7.35) 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:
P h = 4 E r a 3 3R (7.36)
We can use Eq. 7.35 to substitute for a and obtain a relationship between P h and δ h :
P h = 4 E r 3 R 1/2 δ h 3/2 (7.37)
The assumption here is that there is no adhesion between the indenter and the substrate, i.e., G = K I =0 . The fact that there is no stress concentration at the interface is consistent with the fact that the slope of the surface profile of the compliant material is continuous at r=a. This surface profile is plotted in Figure 7.10 and is given by the following expression:[6]
u z = δ π { ( 2- ( r/a ) 2 ) arcsin ( a/r ) +( r/a ) ( 1- ( a/r ) 2 ) 1/2 } (7.38)

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 ( δ 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 ( δ < δ h ) . These deviations from δ h and P h are related by the system compliance, which for this geometry is C 0 as given by Eq. 7.5:
δ - δ h P- P h = C 0 = 1 2 E r a (7.39)
Combination of Eqs. 7.5 and 7.39 gives the following relationship between δ , P and E r :
δ = a 2 3R + P 2 E r a (7.40)
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.
image: 48_home_ken_Mydocs_MSEcore_332_figures_jkrschematic.png
Figure 7.11: An example of the surface profile for adhesive contact for the case where δ 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):
G =- ( P t + P h ) 2 2 dC dA = ( P t + P h ) 2 8 π E r a 3 (7.41)
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 :
a 3 = 3R 4 E r ( P+3 π G R+ ( 6 π G RP+ ( 3 π G R ) 2 ) 1/2 ) (7.42)

7.5 Indentation with Berkovich Trips

Parabolic tips are often used in measurements of adhesion or of the elastic properties of materials. For Hardness measurements tips with sharp corners are more commonly used. One example is the Berkovich tip shown in Figure 7.12.
image: 49_home_ken_Mydocs_MSEcore_332_figures_berkovich_photo.jpg image: 50_home_ken_Mydocs_MSEcore_332_figures_Berkovich_schematic.jpg
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: H= P max A (7.43)
The projected area is related to the contact depth, δ c , by a relationship that depends on the shape of the indenter[9]. For a Berkovich tip the appropriate relationship is:
A=24.5 δ c 2 (7.44)
The procedure for determining the contact depth was developed by Oliver and Pharr, where the following expression is used to estimate the contact depth: δ c = δ max -0.75 P max S (7.45)
where δ 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 δ 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: S= 2 π E r A (7.46)
image: 51_home_ken_Mydocs_MSEcore_332_figures_nanoindentation.svg
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 .


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 ).
image: 52_home_ken_Mydocs_MSEcore_332_figures_tensile_tests_at_diff_temps.svg
Figure 8.1: Typical generic temperature behavior at different temperatures.
Here we are concerned with brittle behavior ( T 1 ) , 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.
image: 53_home_ken_Mydocs_MSEcore_332_figures_Fracture_modes_v2.svg
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.
image: 54_home_ken_Mydocs_MSEcore_332_figures_hole_force_lines.svg
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:
σ max = σ 0 ( 1+2 a c b c ) (8.1)
Note that we recover the behavior described above for a circular whole, where a c = b c and σ max / σ 0 =3. We can also write this in terms of the radius of curvature of the ellipse, ρ c , at the point of maximum stress:
ρ c = b c 2 a c (8.2)
Combination of Eqs. 8.1 and 8.2 gives:
σ max = σ 0 ( 1+2 a c / ρ c ) (8.3)
We are usually interested in very sharp cracks, where a c / ρ c 1 . In this case we can ignore the factor of 1 in Eq. 8.3 and we get the following proportionality:
σ max σ 0 a c (8.4)
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.
image: 55_home_ken_Mydocs_MSEcore_332_figures_elliptical_crack_schematic.svg
Figure 8.4: Elliptical crack with a crack tip radius of curvature, ρ 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:
σ = K 2 π d f( θ ) (8.5)
where K is the stress intensity factor, d is the distance from the crack tip and f( θ ) is some function of the angle θ that reduces to 1 for the direction directly in front of a crack ( θ =0 ). Different functional forms exist for f( θ ) for the different stress components σ xx , σ 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).
image: 56_home_ken_Mydocs_MSEcore_332_figures_crack_axes.svg
image: 57_home_ken_Mydocs_MSEcore_332_figures_crack_geometry_polar.svg
Figure 8.5: Cartesian (a) and polar (b) coordinate axes use d to define stresses in the vicinity of a crack tip.

Mode I loading

The stresses in the vicinity of a mode I crack are given by the following[23]:
( σ 11 σ 22 σ 12 )= K I 2 π d cos θ 2 ( 1- sin θ 2 sin 3 θ 2 1+ sin θ 2 sin 3 θ 2 cos 3 θ 2 sin θ 2 ) (8.6)
This compact notation is used to specify the three relevant values of f( θ ) . For example, for σ 11 we have the following:
σ 11 =( K I / 2 π d ) cos ( θ /2 ) ( 1- sin θ 2 sin 3 θ 2 ) (8.7)
These expressions assume that the crack tip is very sharp, with a very small radius of curvature, ρ c . If d is comparable to ρ c , these equations no longer apply. Consider for example, the presence of an internal crack of length a c and radius of curvature ρ c in a thin sheet of material, shown schematically in Figure 8.6. In this case the stress at the crack edge is σ max as given by Eq. 8.3. An assumption in the use of Eq. 8.6 is that the stresses are substantially less than σ 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, ρ c .
The mode I stress intensity factor for this geometry is given by the applied stress, σ 0 and the crack length a c :
K I = σ 0 π a c (8.8) For values of d that are substantially larger than ρ c but smaller than a c , we can determine the stresses from Eq. 8.6, with K I as given by Eq. 8.8.
image: 58_home_ken_Mydocs_MSEcore_332_figures_crack_geometry.svg
Figure 8.6: An internal crack in a homogeneous solid.

Mode II loading

For mode II loading the crack tip stress fields are given by the following set of expressions[23]:
( σ 11 σ 22 σ 12 )= K II 2 π d ( - sin θ 2 ( 2+ cos θ 2 cos 3 θ 2 ) sin θ 2 cos θ 2 cos 3 θ 2 cos θ 2 ( 1- sin θ 2 sin 3 θ 2 ) ) (8.9)
It is generally difficult to determine K II in a straightforward way, and finite element methods must often be used to determine it for a given loading condition and experimental geometry. Once K II is known, the crack tip stress fields can be obtained from Eq. 8.9.

Mode III loading

While mode III loading is often encountered in practical applications, it is generally avoided in experiments aimed at assessing the fracture behavior of materials, and is not considered further in this text.

8.4 Fracture condition

In the 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 σ 0 to the fracture stress, σ f , and setting K I to K IC in Eq. 8.8 gives:
K IC = σ f π a c (8.10) Rearranging gives:
σ f = K IC / π a c (8.11)
So the fracture stress decreases as the flaw size, a c , increases. This is why a material can appear to be fine, even though small cracks are present in the material. The cracks grow very slowly, but when the reach a critical size for which Eq. 8.11 is satisfied, the material fails catastrophically.
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, σ max f , at the point of fracture is: σ max f 2 σ f a c / ρ c (8.12)
where we have assumed that a c / ρ c 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 σ f . After rearranging we get:
K IC σ max f π 2 ρ c σ max f ρ c (8.13)
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:

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:
G = K I 2 + K II 2 + K III 2 E r (8.14)
Here E r is the reduced modulus that is slightly different for plane stress and plane strain conditions:
E r =E ( Planestressconditions ) E r = E 1- ν 2 ( Planestrainconditions ) (8.15)
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 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:
U f = 1 2 σ f e f = 1 2 σ f 2 E r (8.16)
The stress intensity factor, K I at the fracture point is proportional to the stress, σ 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:
G K I 2 E r (8.17)
This is consistent with Eq. 8.14, but we need to do a more detailed analysis to get the prefactor exactly right.
image: 59_home_ken_Mydocs_MSEcore_332_figures_G-K_relationship.svg
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.
image: 60_home_ken_Mydocs_MSEcore_332_figures_DCB_schematic.svg
Figure 8.8: Double cantilever beam geometry.
For the double cantilever beam geometry the compliance is given by the following expression:
C δ t P t = 8 a c 3 Ew t 3 (8.18)
The crack area, A c is obtained by the crack length, a c by the width of the sample:
A c =w a c (8.19)
So we have:
dC d A c = 1 w dC d a c (8.20)
We now combine Eqs. 8.18 and 8.20 to obtain the following for the energy release rate:
G = P t 2 2w dC d a c = 12 a c 2 P t 2 E w 2 t 3 (8.21)
At fixed load, G increases as the crack length increases - unstable geometry!
can use We Eq. 8.18 to substitute δ t for P t and write the compliance in the following way.
G = 3 δ t 2 t 3 E 16 a c 4 (8.22)
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= δ t in Figure 7.3):
σ zz σ avg =0.5 ( 1- ( r/a ) 2 ) -1/2 (8.23)
Here the average stress, σ avg is defined as follows:
σ avg P π a 2 (8.24)
Note that σ 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:
da-r (8.25)
Substituting d for r in Eq. 8.23 gives:
σ zz σ avg = 1 2 [ 2d/a- ( d/a ) 2 ] -1/2 (8.26)
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 σ zz that is valid near the contact edge: σ zz σ avg 1 2 3/2 ( d a ) -1/2 (8.27)
by comparing to Eq. 8.6 for K I we obtain the following: K I = 1 2 σ avg ( π a ) 1/2 (8.28)
Now we can use the following equation to obtain the following expression for G (assuming K II = K III =0 ):
G = K I 2 2 E r = π a σ avg 2 8 E r = P t 2 8 π E r a 3 (8.29)
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).
image: 61_home_ken_Mydocs_MSEcore_332_figures_flatpunchstress.svg
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.
image: 62_home_ken_Mydocs_MSEcore_332_figures_flat_punch_stress_distributions.svg
Figure 8.10: Dependence of the stress distribution under a flat punch for different values of the confinement ration, a/h .
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 σ 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:
K I = 2 π σ zz a c (8.30)
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):
G = K I 2 2 E r = 2 a c σ zz 2 π E r (8.31)

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:
G IC = K IC 2 E r (8.32)
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:
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.
Table 4: Typical fracture toughness values (plane strain) for different material.
E (GPa)
K IC (MPa m )
G IC J/m 2
High M polystyrene or PMMA
High Impact Polystyrene
Epoxy Resin
Rubber Toughened Epoxy
Glass Filled Epoxy Resin
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 / 2 π d , where d is the distance in front of the crack tip. If the maximum stress is equal to the yield stress, σ y , then this the material must be yielded for values of r that give a stress exceeding σ y . A propagating crack has K I = K IC , so if the size of the plastic zone is h p , we have:
σ y K IC 2 π h p (8.33)
Rearranging gives:
h p K I 2 2 π σ y 2 G IC σ y E 2 π σ y (8.34)
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 ) ? 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 1000 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.
image: 63_home_ken_Mydocs_MSEcore_332_figures_Fuel_Cell.svg
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 ).
image: 64_home_ken_Mydocs_MSEcore_332_figures_zirconia_structures.jpeg
Figure 8.12: Crystal structures of Zirconia (from ref. [#dmsc_machinable_nodate]).
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.
image: 65_home_ken_Mydocs_MSEcore_332_figures_phase_diagram_YO_sml.png
Figure 8.13: Zr O 2 / Y 2 O 3 Phase Diagram (from ref. [2]).
image: 66_home_ken_Mydocs_MSEcore_332_figures_zirconia_transformation_toughening.jpeg
Figure 8.14: Illustration of the transformation toughening mechanism of Zirconia.

8.8.2 Tempered Glass

Glass tempering is a process that leaves the external surfaces of the material in a state of compression. This compressive stress acts against the growth of any defects that exist in the material. The net forces on a material at rest must sum to zero, so if the external portion of the sample is in compression, the internal portions of the sample must be in tension. All of these tensile and compressive stresses store a lot of strain energy, so the sample will fracture in dramatic fashion if a crack manages to make into the compressive region of the sample. This effect is illustrated by “Prince Rupert's drops”, which are formed by rapidly quenching molten glass droplets by dropping them into water. For an entertaining and useful description of tempered glass , see the following video:
image: 67_home_ken_Mydocs_MSEcore_332_figures_tempered_glass.png
Figure 8.15: Comparison of tempered glass and float glass samples that have been fractured. (from http://www.bharatsafety.com/tempered_glass.html).
image: 68_home_ken_Mydocs_MSEcore_332_figures_laminated_glass.jpg
Figure 8.16: Laminated safety glass.

8.8.3 Fiber Composites

image: 69_home_ken_Mydocs_MSEcore_332_figures_fiber_composite_schematic.jpg
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)
image: 70_home_ken_Mydocs_MSEcore_332_figures_dreamliner_materials.gif
Figure 8.18: use of Composites in the Boeing 787 Dreamliner
image: 71_home_ken_Mydocs_MSEcore_332_figures_carbon_fiber_transverse_view.png image: 72_home_ken_Mydocs_MSEcore_332_figures_carbon_fiber_longitudinal_view.png
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).
image: 73_home_ken_Mydocs_MSEcore_332_figures_laminate_plies-700x372.PNG
Figure 8.20: Laminate plies in a carbon fiber composite.
image: 74_home_ken_Mydocs_MSEcore_332_figures_damage_in_multiply_plate.PNG image: 75_home_ken_Mydocs_MSEcore_332_figures_crack_of_resin-600x333.png
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

image: 76_home_ken_Mydocs_MSEcore_332_figures_delam_sizing_problem-700x527.PNG
Figure 8.22: Nondestructive evaluation of damage in a carbon fiber composite.

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:
image: 77_home_ken_Mydocs_MSEcore_332_figures_adhesive_transfer_schematic.png
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:
P s = exp [ - ( p max p cav ) M ] (9.1)
Here M is the Weibull modulus , which is a measure of the distribution of failure probabilities. We generally want M to be large, so that the distribution stresses is very narrow. For example, if M , P s =0 for p max > p cav and P s =1 for p max < p cav .
We can take natural logs of both sides a couple of times to convert Eq. 9.1 to the following:
ln [ ln ( 1/ P s ) ] =M ln p max -M ln p cav (9.2)
This means that we can obtain the Weibull modulus as the slope of a plot of ln [ ln ( 1/ P s ) ] 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 σ as the dependent variable, and with a characteristic tensile stress, σ avg for which the survival probability is 37%:
P s = exp [ - ( σ σ avg ) M ] (9.3)
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 ( -x ) 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:
P f ( σ σ avg ) M (9.4)

Toughening Mechanisms

Factors that increase the hardness or tensile strength ( σ ts ) 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, ρ c (recall that K IC ρ 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:
G C = K IC E r (10.1)
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 :
K I = σ 0 π a = K tip + K s (10.2)
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:
σ zz ( d ) = K tip 2 π d (10.3)
The factor K s describes a reduction in K tip compared to K I because of other internal stresses that are applied in some way to the sample in the vicinity of the crack tip. We'll consider two examples: transformation toughening and crack bridging in fiber reinforced composites. In transformation toughening , compressive stresses result from a phase transition in front of the crack tip to a lower density (higher volume) phase, thereby accommodating some of the strains in the stress field defined by K I , so that K tip < K I ( K s >0). The second example is a conceptually similar crack bridging example, where fibers bridge the crack, (much like craze fibrils in the polymer crazing example discussed above).
image: 78_home_ken_Mydocs_MSEcore_332_figures_toughness_vs_yield_strength.png
Figure 10.1: Schematic representation of toughness vs. tensile strength.

10.1 Hydraulic Fracturing (Fracking)

image: 79_home_ken_Mydocs_MSEcore_332_figures_fracking.jpg
Figure 10.2: The hydraulic fracturing (fracking) process.

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, σ 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 .
image: 80_home_ken_Mydocs_MSEcore_332_figures_tensile_testing.svg
Figure 11.1: Tensile testing sample and representative data for a ductile sample.
Table 5: Shear moduli and values of τ crss for different metals.
G (GPa)
τ crss (MPa)
Magnesium (basal slip)
Magnesium (prism slip)
Titanium (prism slip)
Beryllium (basal slip)
Beryllium (prism slip)

11.1 Critical Resolved Shear Stress

The simplest criterion for the yield of a material is that the resolved shear stress , τ RSS , exceeds a critical value referred to simply as the critical resolved shear stress , τ 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, σ , the angle φ between the slip plane normal and the tensile axis, and the angle λ between the tensile axis and the slip direction:
τ RSS = σ cos φ cos λ (11.1)
In other words, τ RSS = τ CRSS when σ = σ y , where σ y is the tensile yield strength. Substituting τ CRSS for τ RSS and σ y for σ in Eq. 11.1 and solving for σ y gives:
σ y = τ CRSS cos φ cos λ (11.2)
The factor cos φ cos λ is called the Schmid factor , and slip occurs along the slip system that the largest value of this quantity. The sum of φ and λ must be at least 90 , and it can be shown that cos φ cos λ is maximized when φ = λ =4 5 , in which case the yield stress is twice the critical resolved shear stress. We refer to this value of the yield stress as the minimum value, σ y min :
σ y min =2 τ CRSS (11.3)
Measured values of the critical resolved shear stress for different metals are shown in Table 5.
image: 81_home_ken_Mydocs_MSEcore_332_figures_resolved_shear_stress.svg
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:
M= 1 cos φ cos λ (11.4)
Here the average is taken over all possible grain orientations of the material. The tensile yield stress is given by the following expression:
σ y =M τ CRSS = M 2 σ y min (11.5)
This value of M is typically between 2.7 and 3, so we end up with a yield stress in a polycrystalline material that is between 35 and 50 percent larger than the minimum possible value given by Eq. 11.3.

11.2 Yield Surfaces

The full stress state of a material is defined by the 3 principal stresses, σ 1 p , σ 2 p and σ 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

The Tresca yield criterion is the simplest one in that we just assume that shear occurs whenever the maximum shear stress in the sample exceeds some critical value value τ crit . In mathematical terms, yield occurs under the following conditions:
| σ i p - σ j p | max 2 > τ crit (11.6)
where the 'max' subscript indicates that we take the principal stress difference ( σ 1 p - σ 2 p , σ 2 p - σ 3 p or σ 1 p - σ 3 p ) with the largest magnitude. The yield stress, σ y , is typically measured in a uniaxial tensile experiment, where σ 2 p = σ 3 p =0 , and plastic yielding of the material occurs when σ 1 p > σ y . In a uniaxial tensile experiment, the maximum shear stress is half the applied tensile stress, so τ crit = σ y /2 .

11.2.2 Von Mises Yield Criterion

A more complicated criterion is that the Von Mises stress , σ e , stress is larger than some critical value. The Von Mises stress is given as follows:
σ e = 2 2 ( σ 2 p - σ 1 p ) 2 + ( σ 3 p - σ 1 p ) 2 + ( σ 3 p - σ 2 p ) 2 (11.7)
Yielding occurs when σ e > σ y . The Tresca and Von Mises yield surfaces for a two dimensional stress state ( σ 3 p =0) are shown in Figure 11.3. Yielding does not occur inside the surface, but does occur outside the surface.
image: 82_home_ken_Mydocs_MSEcore_332_figures_tresca_von_mises.svg
Figure 11.3: Cross sections of the Tresca and Von Mises yield surfaces at the σ 3 p =0 plane.
Viewed down the hydrostatic line ( σ 1 p = σ 2 p = σ 3 p ), this is what the criteria look like:
image: 83_home_ken_Mydocs_MSEcore_332_figures_tresca_along_hydrostatic_line.svg

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 0 is stretched with a force, P . The engineering tensile stress, σ eng , is obtained by dividing the load by the undeformed cross section, and the true tensile stress, σ t is obtained by dividing the load by the actual cross sectional area of the deformed sample:
image: 84_home_ken_Mydocs_MSEcore_332_figures_tensile_test.svg
Figure 11.4: Uniaxial tensile test.
σ eng =P/ A 0 (11.8) σ true =P/A (11.9)
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:
A= A 0 0 (11.10)
The relationship between the true stress and the engineering stress is therefore as follows:
σ t = σ eng / 0 = σ eng λ (11.11)

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.
image: 85_home_ken_Mydocs_MSEcore_332_figures_considere_samples.svg
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:
dP d λ =0 (11.12)
Since the engineering stress is the load divided by the undeformed cross sectional area, which is a constant, we can write Eq. 11.12 as follows:
d σ eng d λ =0 (11.13)
We can rewrite this expression in terms of the true stress by recognizing that σ eng = σ t / λ (see Eq. 11.11), from which we obtain: λ d σ t d λ - σ t =0 (11.14)
which we rearrange to the following:
d σ t d λ = σ t λ (11.15)
This condition is met when a line drawn fro the origin of a plot of σ t vs. λ 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 λ 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.
d σ t d λ = σ t λ (11.16)
image: 86_home_ken_Mydocs_MSEcore_332_figures_unstable_neck.svg
image: 87_home_ken_Mydocs_MSEcore_332_figures_stable_neck.svg
Figure 11.6: Considére construction for a material that does not form a stable neck (a) and a material that does form a stable neck (b).

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:
σ t =K e t n (11.17)
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).
Strain Hardening Coefficient
Stainless Steel

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 λ (see Eq. 4.37), we have:
σ t =K ( ln ( λ ) ) n (11.18)
Differentiating with respect to λ gives:
d σ t d λ = nK ( ln ( λ ) ) n-1 λ = nK e t n-1 λ
and then using Eq. 11.15 to equate d σ λ /d λ to σ / λ gives the following:
nK e t n-1 λ = σ λ = K e t n λ
We can see that this equation is only true when the following condition holds:
n= ε t (11.19)
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

Deformation twinning is a mechanism of deformation in crystalline materials that doesn't involve dislocation motion. Instead, it involves the concerted, simultaneous motion of may atoms to form a twin boundary in the material. Twin boundaries were discussed in 316-1, and are low-energy, coherent grain boundaries corresponding to a change in stacking sequence in an FCC lattice. The easiest ones to visualize correspond to changes in the stacking sequence of the close-packed layers in face centered cubic (FCC) and hexagonal close packed (HCP) crystals. An example of a twin boundary in and FCC crystal is shown in Figure 11.7.
image: 88_home_ken_Mydocs_MSEcore_332_figures_twin_boundary_fcc.svg
Figure 11.7: Example of a series of twin boundaries in an FCC crystal (from ref.[10], adapted with permission from Macmillan Publishers Ltd.).
image: 89_home_ken_Mydocs_MSEcore_332_figures_deformation_twinning.png
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

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.
image: 90_home_ken_Mydocs_MSEcore_332_figures_frankread.svg
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, b .
  2. Dislocations move along slip planes that contain both 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 b .
  4. The energy per length of a dislocation, T s , is proportional to G b 2 , where b is the magnitude of 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:

12.2 Interactions Between Dislocations

Dislocations interact through the long range strain fields induced by the dislocations. Consider two screw dislocations with Burgers vectors b 1 and b 2 that separated by d . The stress at dislocation 2 due to the presence of dislocation 1 is given by the following expression: τ = G b 2 π d (12.1)
We use a vector notation for the stress here to remind us that the force associated with the shear stress is directed along b .
This stress induces a force on the dislocation, given by the following expression:
F s τ = G b 2 b 1 2 π d (12.2)
If b 1 and 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 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 , 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.
image: 91_home_ken_Mydocs_MSEcore_332_figures_edge_dislocation_interactions.svg
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:
F s τ = G b 2 b 1 2 π ( 1- ν ) d (12.3)
The expression is the same as the expression for the screw dislocations, with the extra factor of 1- ν .
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.
image: 92_home_ken_Mydocs_MSEcore_332_figures_Edge_dislocation_overlapping_stress_fields.svg
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.
image: 93_home_ken_Mydocs_MSEcore_332_figures_dislocation_interaction_map.svg
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 s ˆ pointing into the plane of the figure. The other dislocation is assumed to have the same value for s ˆ and 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 .
image: 94_home_ken_Mydocs_MSEcore_316-1_figures_Edge_dislocation_with_impurity.svg
Figure 12.5: Favored locations the segregation of substitutional impurities to the core of an edge dislocation.

12.3 Work Hardening

(Required supplementary reading:
Wikipedia article on work hardening[22]).
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 G/6 , implying a tensile strength of 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.
image: 95_home_ken_Mydocs_MSEcore_332_figures_strength_vs_dislocation_density.svg
Figure 12.6: Qualitative relationship between strength and dislocation density.

12.4 Grain Boundary Strengthening

(Required supplementary reading:
Wikipedia article on grain boundary strengthening[15]).
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, σ y and the grain size, d :
σ y = σ i + k y d -1/2 (12.4)
Here σ 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

(Supplementary reading:
Wikipedia article on precipitation hardening[18]).
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

(Required supplementary reading:
Wikipedia article on solid solution strengthening[19]).
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.
image: 96_home_ken_Mydocs_MSEcore_332_figures_melting_points.png
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.
image: 97_home_ken_Mydocs_MSEcore_332_figures_creep_load_and_stress_control.svg
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, ε ˙ 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.
image: 98_home_ken_Mydocs_MSEcore_332_figures_creep_rupture_data.png
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:
ε ˙ II =f( σ ) exp ( - Q c k B T ) (13.1)
For the specific case of power law creep, the stress dependence is given by a simple power law:
f( σ ) =A σ n (13.2)
where n is an empirically determined exponent, with typical values ranging between 1 and 7.
image: 99_home_ken_Mydocs_MSEcore_332_figures_creep_regimes.svg
Figure 13.4: Commonly observed creep regimes, illustrating the steady state creep rate, ε ss ˙ .
image: 100_home_ken_Mydocs_MSEcore_332_figures_creep_of_TiO2.svg
Figure 13.5: Creep behavior of Ti O 2 (from ref.[5]).
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.
image: 101_home_ken_Mydocs_MSEcore_332_figures_creep_and_diffusion_activation_energies.svg
Figure 13.6: Activation energies for creep (for T/ T m >0.5) 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, ρ 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 :
. de dt | glide = ρ d V g b (13.3)
Dislocation glide is a thermally activated process, and therefore be described by a model similar to the Eyring model discussed in Section 16.4:
V g sinh ( ( σ - σ crit ) v/2 k B T ) exp ( - Q glide / k B T ) (13.4)
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 σ > σ 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.
image: 102_home_ken_Mydocs_MSEcore_332_figures_n-h_creep.svg
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 :
C v 0 = 1 Ω exp ( - G v k B T ) (13.5)
where Ω 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 σ 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:
C v t = C v 0 exp ( σ Ω / k B T ) C v c = C v 0 exp ( - σ Ω / k B T ) (13.6)
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:
J v =- D v d C v dx (13.7)
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 :
d C v dx C v t - C v c d (13.8)
Combining Eqs.13.7 and13.8 gives:
J v - D v ( C v t - C v c d ) (13.9)
In order to relate the flux (units of vacancies/m 2 s) to the velocity of the grain edge (in m/s), we need to multiply J v by the atomic volume, Ω :
v= Ω J v Ω D v ( C v t - C v c d ) (13.10)
Now we get the strain rate by dividing this velocity by d:
de dt = Ω J v d Ω D v ( C v t - C v c d 2 ) = 2 D v d 2 sinh ( σ Ω k B T ) exp ( - G v k B T ) (13.11)
For real systems the stresses are small enough so that σ Ω k B T , so we can use the fact that sinh ( x ) x for small x to obtain:
de dt 2 D v d 2 σ Ω k B T exp ( - G v k B T ) (13.12)
Vacancy diffusion is a thermally activated process, so we can express the vacancy diffusion coefficient in the following form:
D v = D 0 v exp ( - Q D v k B T ) (13.13)
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:
Q vol = Q v + G v (13.14)
Combination of Eqs. 13.12-13.14 gives the following for the strain rate:
de dt 2 D 0 d 2 σ Ω k B T exp ( - Q vol k B T ) (13.15)
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:
de dt 2 D gb d gb d 3 σ Ω k B T exp ( - G v k B T ) (13.16)
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 :
D gb = D 0 gb exp ( - Q D gb k B T ) (13.17)
Q gb = Q v + G v (13.18)
de dt 2 D 0 gb d gb d 3 σ Ω k B T exp ( - Q gb k B T ) (13.19)
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.
image: 103_home_ken_Mydocs_MSEcore_332_figures_Edge_dislocation_climb.svg
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:
de dt =A σ 3 sinh ( B σ 1.5 / k B T ) exp ( - Q vol / k B T ) (13.20)

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:
image: 104_home_ken_Mydocs_MSEcore_332_figures_deformation_mechanism_map.svg
Figure 13.9: Deformation mechanism map.

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.
image: 105_home_ken_Mydocs_MSEcore_332_figures_amorhpous_polymers.svg
Figure 14.1: Structure of an amorphous polymer.
image: 106_home_ken_Mydocs_MSEcore_332_figures_semicrystalline_structure.gif
Figure 14.2: Structure at different length scales for a semicrystalline polymer.
image: 107_home_ken_Mydocs_MSEcore_332_figures_polymers_classification_scheme.svg
Figure 14.3: Polymers classification scheme.
image: 108_home_ken_Mydocs_MSEcore_332_figures_tensile_data_for_pvc_and_hdpe.svg
Figure 14.4: Load-Displacement and true stress-true strain curves for PVC (polyvinyl chloride) and HDPE (high density polyethylene).
image: 109_home_ken_Mydocs_MSEcore_332_figures_PMMA_compression_and_tension.svg
Figure 14.5: Stress-strain curves for PMMA (polymethylmethacrylate) under tensile and compressive loading conditions.
image: 110_home_ken_Mydocs_MSEcore_332_figures_PMMA_yield_surface.svg
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

image: 111_home_ken_Mydocs_MSEcore_332_figures_polymer_toughness_vs_M.svg
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. Deformation Mechanisms
Suppose we do a simple stress strain experiment on polystyrene. Polystyrene deforms by one of two different mechanisms:
  1. Shear bands due to strain softening (decrease in true stress after yield in shear).
  2. Crazing - requires net dilation of sample (fracture mechanism for PS and PMMA).
image: 112_home_ken_Mydocs_MSEcore_332_figures_shear_zone_schematic.svg
image: 113_home_ken_Mydocs_MSEcore_332_figures_craze_schematic.svg
Figure 14.8:
Crazes are load bearing - but they break down to form cracks - failure of specimen. Crazing
image: 114_home_ken_Mydocs_MSEcore_332_figures_craze_fibrils.svg
Figure 14.9: Conceptual drawing of fibrils at the interface between the crazed and uncrazed material (from ref.[8])
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 σ 1 + σ 2 + σ 3 >0 (crazes have voids between fibrils)
  4. Crazing occurs first for PMMA in uniaxial extension ( σ 2 =0 )
  5. G Ic is determined by energy required to form a craze ( 1000 J/ m 2 )
  6. Crazing requires strain hardening of fibrils - material must be entangled ( M> M c ), M c typically 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.
image: 115_home_ken_Mydocs_MSEcore_332_figures_meniscus_instability.svg
Figure 14.10: Meniscus instability mechanism for the formation of craze fibrils.
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)
image: 116_home_ken_Mydocs_MSEcore_332_figures_crazing_vs_shear_yielding.svg
Figure 14.11: Deformation map for the shear yielding and crazing for plane stress conditions ( σ 3 p =0) Dugdale model
1) Tensile stress throughout plastic zone is constant value, σ c
2) This stress acts to produce a crack opening displacement δ c
G IC = δ c σ c (14.1)
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.
image: 117_home_ken_Mydocs_MSEcore_332_figures_dugdale_stresses.svg
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.
image: 118_home_ken_Mydocs_MSEcore_332_figures_hips_morphology.svg
Figure 14.13: Morphology of high impact polystyrene.
image: 119_home_ken_Mydocs_MSEcore_332_figures_HIPS_multiple_crazes.png
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:
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.
image: 120_home_ken_Mydocs_MSEcore_332_figures_HIPS_tensile_schematic.svg
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)
image: 121_home_ken_Mydocs_MSEcore_332_figures_HIPS_fracture_schematic.svg
Figure 14.16: Schematic stress-strain curves for polystyrene (PS) and high impact polystyrene (HIPS) in the presence of a crack.
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.
image: 122_home_ken_Mydocs_MSEcore_332_figures_charpy_schematic.svg Figure 14.17: Charpy impact test.
For most materials the fracture toughness is rate dependent, but same general features for toughening materials often apply at both high and low fracture rates. For example, high impact polystyrene is much tougher than polystyrene at high strain rates, for the same general reasons outlined in Section 14.1.
image: 123_home_ken_Mydocs_MSEcore_332_figures_izod_geometry.svg
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, σ , enables us to define a time-dependent relaxation modulus , E( t ) :
E( t ) = σ ( t ) e 0 (15.1)
image: 124_home_ken_Mydocs_MSEcore_332_figures_uniaxial_viscoelastic.svg
Figure 15.1: Experimental geometry and time dependent strain for for the determination of the time-dependent relaxation modulus and the frequency-dependent dynamic moduli.
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, ω , and an amplitude, e 0 :
e( t ) = e 0 sin ( ω t ) (15.2)
Note that the strain rate, de dt , is also an oscillatory function, with the same angular frequency, but shifted with respect to the strain by 90 :
de dt = e 0 ω cos ( ω t ) = γ 0 ω sin ( ω t+ π /2 ) (15.3)
The resulting stress is also an oscillatory function with an angular frequency of ω , and is described by its amplitude and by the phase shift of the relative to the applied strain:
σ ( t ) = σ 0 sin ( ω t+ φ ) (15.4)
Now we can define a complex modulus with real and imaginary components as follows:
E * = E +i E =| E * | e i φ (15.5)
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 and E , respectively), or in terms of its magnitude, | E * | and phase, φ . The magnitude of the complex modulus is simply the stress amplitude normalized by the strain amplitude:
| E * | = σ 0 / e 0 (15.6)
The phase angle, φ , describes the lag between the stress and strain in the sample. Examples for φ =9 0 (the maximum value, characteristic of a liquid) and φ =4 5 are shown in Figure 15.2.
image: 125_home_ken_Mydocs_MSEcore_332_figures_sinewave.png
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:
e i φ = cos φ +i sin φ (15.7)
use this expression in Eq. 15.5, we see that the storage modulus , E gives the stress that is in phase with the strain (the solid-like part), and is given by the following expression:
E =| E * | cos φ (15.8)
Similarly, the loss modulus , E , gives the response of the material that is in phase with the strain rate (the liquid-like part):
E =| E * | sin φ (15.9)
We can combine Eqs. 15.8 and 15.9 to get the following expression for tan φ , commonly referred to simply as the loss tangent:
tan φ = E E (15.10)
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]:
T= K t θ +I d 2 θ d t 2 (15.11)
where θ is the rotation angle, T is the torque on the resonator, K t is the torsional stiffness and I is the moment of inertia of the system, which depends on the details of the inertial bar used in the experiment. For a cylindrical fiber K t is obtained from Eq. 5.10:
K t T θ = π G d 4 32 (15.12)
image: 126_home_ken_Mydocs_MSEcore_332_figures_torsional_resonator_no_numbers.svg
Figure 15.3: Torsional resonator.
The resonant angular frequency of the oscillator, ω n is given by the following expression:
ω n = K t I (15.13)
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:
θ = θ 0 cos ( ω n t ) (15.14)
This analysis assumes that the system is entirely elastic. What if the spring has some viscoelastic character to it? We begin by using Euler's formula (Eq. 15.7 ) to rewrite Eq. 15.14 in the following way:
θ = θ 0 Real ( exp ( i ω n t ) ) (15.15)
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:
To understand this last point, suppose write ω n * in the following way:
ω n * = ω n +i ω n (15.16)
Putting this back into Eq. 15.15 gives:
θ = θ 0 Real ( exp ( i ω n t ) ) exp ( - ω n t ) (15.17)
Using the Euler formula (Eq. 15.7) to expand exp ( i ω n t ) and taking the real part gives:
θ = θ 0 cos ( ω n t ) exp ( - ω n t ) (15.18)
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, σ 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 , η , as defined in the following expression:
σ = η de dt (15.19)
Schematic representations of springs and dashpots are shown in Figure 15.4.
image: 127_home_ken_Mydocs_MSEcore_332_figures_springs_and_dashpots.svg
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 μ 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:
image: 128_home_ken_Mydocs_MSEcore_332_figures_falling_ball.svg
Figure 15.5: Terminal velocity of a sphere descending in a viscous fluid.
F=g 4 3 π R 3 ( ρ s - ρ ) (15.20)
Here ρ s is the density of the solid sphere, ρ 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:
F=6 π V η R (15.21)
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:
V= 2 9 g R 2 η ( ρ s - ρ ) (15.22)
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 μ 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 η 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:
σ = η d e 1 dt (15.23)
σ =E e 2 (15.24)
image: 129_home_ken_Mydocs_MSEcore_332_figures_maxwell_model.svg
Figure 15.6: Maxwell model.
The total strain, e is e 1 + e 2 , so the time derivative of the total strain is:
de dt = d e 1 dt + d e 2 dt = σ η + 1 E d σ dt (15.25)
In a stress relaxation experiment the strain instantaneous increased to an initial value of e 0 and we follow the relaxation of the stress at this fixed strain. Because the dashpot cannot respond instantaneously, all of the initial strain is in the spring, i.e., e 0 = e 2 and the initial stress, σ 0 is equal to E e 0 . The solution to Eq. 15.25 for this initial condition and for de/dt=0 is:
σ ( t ) =E e 0 exp ( -t/ τ ) (15.26)
with τ = η /E . We divide by e 0 to obtain the relaxation modulus: E( t ) σ ( t ) / e 0 =E exp ( -t/ τ ) (15.27)
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:
E * ( ω ) = Ei ω τ 1+i ω τ (15.28)
We can obtain the real and imaginary parts of E * by multiplying the top and bottom of Eq. 15.28 by 1+i/( ω τ ) , remembering that i 2 =-1 . Taking the real and imaginary components of E * gives E and E :
E ( ω ) = E ( ω τ ) 2 1+ ( ω τ ) 2 (15.29)
E ( ω ) = E ω τ 1+ ( ω τ ) 2 (15.30)

15.3.2 Standard Linear Solid

For a single maxwell element the relaxation modulus decays to zero at long times ( t τ ). 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( t ) and E * ( w ) that were obtained from the Maxwell model:
E( t ) = E r + E 1 exp ( -t/ τ 1 ) (15.31)
E * ( ω ) = E r + E 1 i ω τ 1 1+i ω τ 1 (15.32)
image: 130_home_ken_Mydocs_MSEcore_332_figures_standard_linear_solid.svg 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( t ) and E * ( ω ) 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( t ) and E * ( w ) are given by the following expressions: E( t ) = E r +j E j exp ( -t/ τ j ) (15.33)
E * ( ω ) = E r +j E j i ω τ j 1+i ω τ j (15.34)
with τ j = η j / E j .
image: 131_home_ken_Mydocs_MSEcore_332_figures_generalized_maxwell.svg 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: σ = η de dt +Ee (15.35)
The Kelvin-Voigt model is the simplest model for the description of a creep experiment, where the stress jumps instantaneously from 0 to σ = σ 0 and we track the strain as a function of time. We have:
e( t ) = σ 0 E ( 1- exp ( t/ τ ) ) ; τ = η /E (15.36) The dynamic modulus for the Kelvin-Voigt element is given by the following expression:
E * =E+i ω η (15.37)
So that E is simply equal to E, and E is equal to ω η .
image: 132_home_ken_Mydocs_MSEcore_332_figures_Kelvin_Voight_model.svg Figure 15.9: Kelvin-Voigt Model.

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, σ and the time, t , since the load was applied.

16.1 Creep in the Linear Viscoelastic Regime

In the linear viscoelastic regime the strain is linearly dependent on the applied stress, σ , allowing us to define the creep compliance function , J( t ) , in the following way:
J( t ) e( σ ,t ) σ (16.1)
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 σ at t=0 can be represented as the sum of three contributions: e 1 , e 2 , and e 3 :
e( σ ,t ) = e 1 ( σ ) + e 2 ( σ ,t ) + e 3 ( σ ,t ) (16.2)
image: 133_home_ken_Mydocs_MSEcore_332_figures_linear_creep_model.svg
Figure 16.1: Linear creep model.
In this equation e 1 corresponds to an instantaneous elastic strain, e 2 is the recoverable viscoelastic strain and e 3 is the plastic strain, with the three different components illustrated in Figure 16.2, and given by the following expressions:
e 1 = σ E 1 e 2 = σ E 2 [ 1- exp ( -t/ τ 2 ) ] e 3 = σ η 3 t (16.3)
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.
image: 134_home_ken_Mydocs_MSEcore_332_figures_creep_model.svg
Figure 16.2: Strain contributions in a nonlinear creep model.

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:
e( σ ,t ) =f( σ ) J( t ) (16.4)
If this sort of separability holds, then it is possible to make predictions of creep based on limited experimental data. The procedure is illustrated in Figure 16.3 and involves the following steps:
  1. Measure e( t ) at different different stresses ( σ 1 and σ 2 ) in Figure 16.3. If the ratio e( σ 1 ,t ) /e( σ 2 ,t ) 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( t ) .
  2. To get the stress-dependent function, f( σ ) , it is sufficient to make a series of measurements at a single, experimentally convenient time.
image: 135_home_ken_Mydocs_MSEcore_332_figures_creep_data.svg 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( σ 1 ,t ) /e( σ 2 ,t ) 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:
e 1 = σ E 1 e 2 = C 1 σ n [ 1- exp ( - C 2 t ) ] e 3 = C 3 σ n t (16.5)
Not that the material behavior is specified by 5 constants, E, C 1 , C 2 , C 3 ,n , that we need to obtain by fitting to actual experimental data. For a linear response (n=1 ) we can make a connection to the spring and dashpot models described earlier. In this case the behavior of the material is represented by the model of linear viscoelastic elements shown in Figure 16.4, and the constants appearing in Eq. 16.5 correspond to the following linear viscoelastic elements from Figure 16.1:
E= E 1 C 1 =1/ η 2 C 2 = E 2 / η 2 C 3 =1/ η 3 (16.6)
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.
image: 136_home_ken_Mydocs_MSEcore_332_figures_nonlinear_creep_model.svg
Figure 16.4: Nonlinear creep model.

16.4 Eyring Model of Steady State Creep

The Eyring rate model of creep is a very general model aimed at understanding the effect of the stress on the flow properties of a material. It can be used to describe a very wide range of materials, and is based on the modification of the activation energy for material for flow by the applied stress.

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 σ v , where σ 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 σ /2 and to increase the activation barrier in the opposite direction by this same amount.
image: 137_home_ken_Mydocs_MSEcore_332_figures_Eyring.svg
Figure 16.5: Effect of an applied stress on a thermally activated creep process. 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:
f 1 = f 0 exp ( - Q * -v σ /2 k B T ) = exp ( - Q * k B T ) exp ( v σ 2 k B T ) (16.7)
f 2 = f 0 exp ( - Q * +v σ /2 k B T ) = exp ( - Q * k B T ) exp ( -v σ 2 k B T ) (16.8)
The net rate of strain is proportional to f 1 - f 2 , the net frequency of hops in the forward direction:
de dt =A( f 1 - f 2 ) =A f 0 exp ( - Q * k B T ) [ exp ( v σ 2 k B T ) - exp ( -v σ 2 k B T ) ] (16.9)
Where A is a dimensionless constant. Using the following definition of the hyperbolic sine function (sinh):
sinh ( x ) =( e x - e -x ) /2 (16.10)
we obtain the following expression for the strain rate:
de dt =2A f 0 exp ( - Q * k B T ) sinh ( v σ 2 k B T ) (16.11)
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 , sinh xx , and for large x, sinh ( x ) 0.5 exp ( x ) .
image: 138_home_ken_Mydocs_MSEcore_332_figures_sinhplot.svg
Figure 16.6: Behavior of the hyperbolic sine function. Low Stress Regime
In the low-stress regime we can use the approximation sinh ( x ) x to get the following expression, valid for v σ k B T .
de dt = A f 0 v σ k B T exp ( - Q * k B T ) (16.12)
In this regime the strain rate is linear in stress, which means that we can define a stress-independent viscosity from the following expression: de dt = σ η (16.13)
Comparing Eqs.16.12 and16.13 gives the following for the viscosity:
η = k B T A f 0 v exp ( Q * k B T ) (16.14)
So the Eyring theory reduces to and Arrhenius viscosity behavior in the linear, low-stress regime. High Stress Regime
In the high-stress regime, we use the fact that sinh ( x ) exp ( x ) 2 for large x to obtain the following expression for v σ k B T :
de dt =A f 0 exp ( - Q * k B T ) exp ( v σ 2 k B T ) (16.15)
Equivalently, we can write the following: de dt =A f 0 exp ( - Q * -v σ /2 k B T ) (16.16)
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 ( de/dt ) vs. σ , with the slope being equal to v/2 k B T . 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.
image: 139_home_ken_Mydocs_MSEcore_332_figures_Lec06creep-img19.svg 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 .
image: 140_home_ken_Mydocs_MSEcore_332_figures_eyring_creep_model.svg Figure 16.8: Nonlinear creep model, including an Eyring rate element for the steady-state creep component.

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.
image: 141_home_ken_Mydocs_MSEcore_332_figures_space_elevator_schematic.png
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):
F=m ω 2 r- G gr M e m/ r 2 (17.1)
where G gr is the gravitational constant and M e is the mass of the earth. The angular velocity of the earth is 2 π radians per day, or in more useful units:
ω = 2 π ( 24hr ) ( 3600s/hr ) =7.3x1 0 -5 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):
g 0 = G gr M e r 0 2 (17.2)
The net force on the mass at r can be written as follows:
F=m ω 2 r-m g 0 ( r 0 r ) 2 (17.3)
The net force is zero for r= r s (geosynchronous orbit):
r s = ( g 0 r 0 2 ω 2 ) 1/3 =4.1x1 0 7 m(22,000 miles ) (17.4)
This is an easier number to remember than the angular velocity of the earth, so we use this expression for r s to eliminate ω from Eq. 17.3, obtaining the following:
F=m g 0 r 0 2 ( r r s 3 - 1 r 2 ) (17.5)
image: 142_home_ken_Mydocs_MSEcore_332_figures_earth_schematic_for_elevator.svg
Figure 17.2: Radial forces acting on an orbiting mass.
Now consider a cable that extends from the earths surface ( r= r 0 ) to a distance r 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 ) we're in good shape. The mass increment for a cable of length dr is ρ Adr , where A is the cross sectional area of the cable and ρ 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:
F 0 = ρ A g 0 r 0 2 r 0 r ( r r s 3 - 1 r 2 ) r (17.6)
After integration we get:
F 0 = ρ A g 0 r 0 2 [ 1 2 r s 3 ( r 2 - r 0 2 ) +( 1 r - 1 r 0 ) ] (17.7)
image: 143_home_ken_Mydocs_MSEcore_332_figures_space_elevator_geometry.svg
Figure 17.3: A cable extending from the surface of the earth to a distance r away from the earth's center.
F 0 is plotted in Figure 17.4:
image: 144_home_ken_Mydocs_MSEcore_332_figures_space_elevator_cable.svg
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 :
F max = ρ A g 0 r 0 2 [ ( 1 r - 1 r s ) + 1 2 r s 3 ( r 2 - r s 2 ) ] (17.8)
We have the following numbers:
From these numbers we get F max ρ A = σ max ρ =4.8x1 0 7 N/ m 2 Kg/ m 3
Let's compare that to the best materials that are actually available. An Ashby plot of tensile strength ( σ 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 σ f / ρ . The line drawn on Figure 17.5 corresponds to σ f / ρ is 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.
image: 145_home_ken_Mydocs_MSEcore_332_figures_materials_availability_2010.svg 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 ) :
A s A 0 = exp ( 0.77 r 0 ρ g 0 σ ) (17.9)
If we assume σ / ρ =2.8x1 0 6 N/ m 2 Kg/ m 3 , so that the system is operating at the value of σ /f corresponding to the solid line in Figure 17.5 gives 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 10 cm at r= r s . We don't have much leeway in decreasing σ / ρ , however. If the best we can do is σ / ρ =1.0x1 0 6 N/ m 2 Kg/ m 3 , we get 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

18.2 Linear Elasticity

18.3 Fracture Mechanics

18.4 Yield Behavior

18.5 Viscoelasticity

18.6 Creep

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 φ ( x,y,z ) . In each element φ 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 φ are determined, in combination with the assumed spatial variation of φ within each element, the approximated field is completely determined within that element. Thus the function φ ( x,y,z ) 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:
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 σ can be related to nodal forces f 1 and f 2 by free-body diagrams, A σ + f 1 =0,-A σ + f 2 =0. (19.1) In turn, σ is related to elastic modulus E and axial strain ϵ , σ =E ϵ , ϵ = u 2 - u 1 L . (19.2) Thus, we obtain AE L ( u 1 - u 2 ) = f 1 (19.3) AE L ( u 2 - u 1 ) = f 2 (19.4) or, in matrix form,
[ k ] d = f (19.5) Here [ k ] is the element stiffness matrix: [ k ] =[ k -k -k k ] where k= AE L , (19.6) and d and f are the element nodal displacement and the element nodal force vectors, respectively, d ={ u 1 u 2 }, f ={ f 1 f 2 }. (19.7) 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.
image: 146_home_ken_Mydocs_MSEcore_332_figures_fig_bar_element.png
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 σ 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 ( 1 ) , A ( 1 ) , E ( 1 ) , L ( 2 ) , A ( 2 ) , E ( 2 ) , and L ( 3 ) , A ( 3 ) , E ( 3 ) , respectively. Since abrupt discontinuities in cross-Sectional area will lead to stress concentrations at the segment junctions, we will let A ( 1 ) = A ( 2 ) = A ( 3 ) . The structure is attached to a rigid wall at its left end, and a force 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 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 D and the nodal force vector F are: D ={ u 1 u 2 u 3 u 4 }, F ={ f 1 f 2 f 3 f 4 }. (19.8) The stiffness equation for the overall structure can be written as F =[ K ] D (19.9) where [ K ] is the 4 × 4 global stiffness matrix.
image: 147_home_ken_Mydocs_MSEcore_332_figures_fig_bar_structure.png
Figure 19.2: (a) A bar consisting of three segments, with stiffnesses k ( 1 ) = E ( 1 ) A ( 1 ) / L ( 1 ) , k ( 2 ) = E ( 2 ) A ( 2 ) / L ( 2 ) and k ( 3 ) = E ( 3 ) A ( 3 ) / L ( 3 ) , 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: k ( e ) d ( e ) = f ( e ) , (19.10) so that [ k ( 1 ) - k ( 1 ) - k ( 1 ) k ( 1 ) ]{ u 1 ( 1 ) u 2 ( 1 ) } = { f 1 ( 1 ) f 2 ( 1 ) }, [ k ( 2 ) - k ( 2 ) - k ( 2 ) k ( 2 ) ]{ u 2 ( 2 ) u 3 ( 2 ) } = { f 2 ( 2 ) f 3 ( 2 ) }, (19.11) [ k ( 3 ) - k ( 3 ) - k ( 3 ) k ( 3 ) ]{ u 3 ( 3 ) u 4 ( 3 ) } = { f 3 ( 3 ) f 4 ( 3 ) }. Note that each element quantity is denoted by a superscript ( e ) , 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 ( e ) = E ( e ) A ( e ) / L ( e ) .
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 ( 1 ) = u 2 ( 2 ) , and u 3 ( 2 ) = u 3 ( 3 ) . 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: [ k ( 1 ) - k ( 1 ) 0 0 - k ( 1 ) k ( 1 ) 0 0 0 0 0 0 0 0 0 0 ]{ u 1 u 2 u 3 u 4 } = { f 1 ( 1 ) f 2 ( 1 ) 0 0 }, [ 0 0 0 0 0 k ( 2 ) - k ( 2 ) 0 0 - k ( 2 ) k ( 2 ) 0 0 0 0 0 ]{ u 1 u 2 u 3 u 4 } = { 0 f 2 ( 2 ) f 3 ( 2 ) 0 }, (19.12) [ 0 0 0 0 0 0 0 0 0 0 k ( 3 ) - k ( 3 ) 0 0 - k ( 3 ) k ( 3 ) ]{ u 1 u 2 u 3 u 4 } = { 0 0 f 3 ( 3 ) f 4 ( 3 ) }. 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 ( 1 ) + f 2 ( 2 ) , and f 3 = f 3 ( 2 ) + f ( 3 ) , where the absence of a superscript denoting the element number indicates a global or structure nodal force. The global force vector F is therefore F = f ( 1 ) + f ( 2 ) + f ( 3 ) =( k ( 1 ) + k ( 2 ) + k ( 3 ) ) D = K D (19.13) 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 K =[ k ( 1 ) - k ( 1 ) 0 0 - k ( 1 ) k ( 1 ) + k ( 2 ) - k ( 2 ) 0 0 - k ( 2 ) k ( 2 ) + k ( 3 ) - k ( 3 ) 0 0 - k ( 3 ) k ( 3 ) ], (19.14) and the global stiffness equation is [ k ( 1 ) - k ( 1 ) 0 0 - k ( 1 ) k ( 1 ) + k ( 2 ) - k ( 2 ) 0 0 - k ( 2 ) k ( 2 ) + k ( 3 ) - k ( 3 ) 0 0 - k ( 3 ) k ( 3 ) ]{ u 1 u 2 u 3 u 4 }={ f 1 f 2 f 3 f 3 }. (19.15) 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 [ k ( 1 ) + k ( 2 ) - k ( 2 ) 0 - k ( 2 ) k ( 2 ) + k ( 3 ) - k ( 3 ) 0 - k ( 3 ) k ( 3 ) ]{ u 2 u 3 u 4 }={ 0 0 F } (19.16) 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 σ ( e ) = E ( e ) ϵ ( e ) . The element axial strain ϵ ( e ) = e ( e ) / L ( e ) , where e ( e ) is the elongation of the element. For example, the average axial stress of element 2 can be determined as σ ( 2 ) = E ( 2 ) L ( 2 ) ( u 3 - u 2 ) . (19.17)

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 [ k ( 1 ) - k ( 1 ) 0 0 - k ( 1 ) k ( 1 ) + k ( 2 ) - k ( 2 ) 0 0 - k ( 2 ) k ( 2 ) + k ( 3 ) 0 0 0 - k ( 3 ) k ( 3 ) ]{ u 1 u 2 u 3 u 4 }={ - k ( 1 ) u 2 0 0 F }. (19.18) Note that we have recovered thereaction force at the support: f 1 =- k ( 1 ) 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 [ u i v i ] T , each nodal force f i in the force vector becomes a vector [ f ix f iy ] 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
D = [ u 1 v 1 u 2 v 2 u 3 v 3 u 4 v 4 u 5 v 5 ] T F = [ f 1x f 1y f 2x f 2y f 3x f 3y f 4x f 4y f 5x f fy ] T (19.19)
while the stiffness matrix K that relates the displacements and forces in the relation K D = F is a 10 × 10 matrix.
image: 148_home_ken_Mydocs_MSEcore_332_figures_overhead_hoist.png
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:
σ x,x + τ xy,y + F x =0 τ xy,y + σ y,y + F y =0 (19.20)
or, in matrix form, T σ + F = 0 , (19.21) where =[ x 0 0 y y x ], σ ={ σ x σ y τ xy }, F ={ F x F y }. (19.22)
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.
image: 147_home_ken_Mydocs_MSEcore_332_figures_fig_bar_structure.png
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, σ = E ϵ . (19.23) In two-dimensions, the strain vector ϵ = [ ϵ x ϵ y γ xy ] T , and E =[ E 11 E 12 E 13 E 21 E 22 E 23 E 31 E 32 E 33 ]. (19.24) For example, in isotropic, plane stress conditions, E takes the form E = E 1- ν 2 [ 1 ν 0 ν 1 0 0 0 1- ν 2 ], (19.25) where E is the elastic modulus and ν is the Poisson's ratio of the material.
Note that the unknown field quantity are the displacements. The strain-displacement relation, ϵ = u , together with the constitutive equation (23), relates the stresses in the element to the displacements, where u = [ u v ] 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. Linear Interpolation
Linear interpolation between two points ( x 1 , u 1 ) and ( x 2 , u 2 ) . The linear interpolating function u ˆ is u ˆ ( x ) = a 0 + a 1 x, (19.26) so that
u 1 = a 0 + a 1 x 1 u 2 = a 0 + a 1 x x (19.27)
Solving for a 0 and a 1 and rearranging, we obtain u ˆ ( x ) = N 1 ( x ) u 1 + N 2 ( x ) u 2 = N u , (19.28) where N =[ N 1 N 2 ], u ={ u 1 u 2 }, (19.29) and N 1 ( x ) = x 2 -x x 2 - x 1 , N 2 ( x ) = x- x 1 x 2 - x 1 , (19.30) are the two linear shape functions.
Figure 19.5(a) shows a two-node element with linear interpolation of the unknown field u( x ) 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
image: 149_home_ken_Mydocs_MSEcore_332_figures_fig_shape.png
Figure 19.5: (a) Linear interpolation and shape functions. (b) Quadratic interpolation and shape functions. Quadratic Interpolation
Quadratic interpolation fits a quadratic function to the three points ( x 1 , u 1 ) , ( x 2 , u 2 ) , and ( x 3 , u 3 ) . The interpolation function is u ˆ ( x ) = N 1 ( x ) u 1 + N 2 ( x ) u 2 + N 3 ( x ) u 3 = N u , (19.31) where N =[ N 1 N 2 N 3 ], u ={ u 1 u 2 u 3 }, (19.32) and the three shape functions are N 1 = ( x 2 -x ) ( x 3 -x ) ( x 2 - x 1 ) ( x 3 - x 1 ) , N 2 = ( x 1 -x ) ( x 3 -x ) ( x 1 - x 2 ) ( x 3 - x 2 ) , N 3 = ( x 1 -x ) ( x 2 -x ) ( x 1 - x 3 ) ( x 2 - x 3 ) (19.33) Figure 19.5(b) shows a three-node element with quadratic interpolation of the unknown field u( x ) 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 . 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( x,y ) , and in the y -direction v( x,y ) .
image: 150_home_ken_Mydocs_MSEcore_332_figures_fea2_3.png
Figure 19.6: A linear triangle element.
The linear triangle uses linear interpolations for both u and v , we will call the interpolating functions u ˆ and v ˆ :
u ˆ = a 1 + a 2 x+ a 3 y v ˆ = a 4 + a 5 x+ a 6 y (19.34)
From the values of u and y at the 3 nodes, equation19.34 can be rearranged as
u ˆ = N 1 u 1 + N 2 u 2 + N 3 u 3 v ˆ = N 1 v 1 + N 2 v 2 + N 3 v 3 (19.35)
where the shape functions are
N 1 ( x,y ) =1- x x 2 + ( x 3 - x 2 ) y x 2 y 3 N 2 ( x,y ) =- x x 2 - x 3 y x 2 y 3 N 3 ( x,y ) = y y 3 (19.36)
Notice that the same shape functions are used for both u ˆ and v ˆ . 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( x,y ) and v( x,y ) are again interpolated in both x - and y -directions with the same linear variation, so we have:
u ˆ ( x,y ) = N 1 ( x,y ) u 1 + N 2 ( x,y ) u 2 + N 3 ( x,y ) u 3 + N 4 ( x,y ) u 4 v ˆ ( x,y ) = N 1 ( x,y ) v 1 + N 2 ( x,y ) v 2 + N 3 ( x,y ) v 3 + N 4 ( x,y ) v 4 (19.37)
where the shape functions are given by
N 1 ( x,y ) = ( a-x ) ( b-y ) 4ab , N( x,y ) = ( a+x ) ( b-y ) 4ab N 3 ( x,y ) = ( a+x ) ( b+y ) 4ab , N( x,y ) = ( a-x ) ( b+y ) 4ab (19.38)
image: 151_home_ken_Mydocs_MSEcore_332_figures_fig_bilinear.png
Figure 19.7: A bilinear quadrilateral element. 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.
image: 152_home_ken_Mydocs_MSEcore_332_figures_linear_quadratic_elements.png
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( x ) 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.
image: 153_home_ken_Mydocs_MSEcore_332_figures_fig_convergence.svg 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.
image: 154_home_ken_Mydocs_MSEcore_332_figures_cantilever.png