Open Educational Resources

Linear Elastic Materials: Examples and Problems

Example 1

The state of strain inside a continuum is represented by the small strain matrix:

    \[\varepsilon = \begin{pmatrix} 0.002 & 0.001 & 0.002 \\ 0.001 & -0.003 & 0.001 \\ 0.002 & 0.001 & 0.0015 \\ \end{pmatrix}\]

Assume that the material is steel with Young’s modulus of 210 GPa and Poisson’s ratio of 0.3. Assume that the material follows the von Mises yield criterion with \sigma_{max} = 750 MPa. Find the principal strains and the principal strains directions, the principal stresses and their directions, and comment on whether this state of strain is possible before yielding of the material.

Solution

The principal strains (eigenvalues) are

    \[\varepsilon_{1} = 0.004 \hspace{5mm} \varepsilon_{2} = -0.0033 \hspace{5mm} \varepsilon_{3} = -0.00026\]

and the corresponding principal directions (eigenvectors) of the strains are:

    \[ev_{\varepsilon 1} = \begin{pmatrix} 0.7324 \\ 0.1964 \\ 0.6519 \\ \end{pmatrix} \hspace{5mm} ev_{\varepsilon 2} = \begin{pmatrix} 0.1282 \\ -0.9801 \\ 0.1512 \\ \end{pmatrix} \hspace{5mm} ev_{\varepsilon 3} = \begin{pmatrix} -0.6687 \\ 0.02718 \\ 0.7431 \\ \end{pmatrix}\]

To apply the linear elastic constitutive relations, the strain matrix will be converted into a vector of strains, as follows:

    \[\varepsilon = \begin{pmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ 2\varepsilon_{12} \\ 2\varepsilon_{13} \\ 2\varepsilon_{23} \\ \end{pmatrix} = \begin{pmatrix} 0.002 \\ -0.003 \\ 0.0015 \\ 0.002 \\ 0.004 \\ 0.002 \\ \end{pmatrix}\]

Using Equation 7, the corresponding stress vector is:

    \[\sigma = \begin{pmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{12} \\ \sigma_{13} \\ \sigma_{23} \\ \end{pmatrix} = \begin{pmatrix} 383.66 \\ -424.04 \\ 302.89 \\ 161.54 \\ 323.08 \\ 161.54 \\ \end{pmatrix} MPa\]

The corresponding stress matrix is:

    \[\sigma = \begin{pmatrix} 383.66 & 161.54 & 323.08 \\ 161.54 & -424.04 & 161.54 \\ 323.08 & 161.54 & 302.89 \\ \end{pmatrix} MPa\]

The principal stresses (eigenvalues) of the stress matrix are:

    \[\sigma_1 = 714.53 MPa \hspace{5mm} \sigma_2 = -470.1 MPa \hspace{5mm} \sigma_3 = 18.064 MPa\]

the principal directions (the corresponding eigenvalues of the stress matrix) are:

    \[ev_{\sigma 1} = \begin{pmatrix} 0.7324 \\ 0.1964 \\ 0.6519 \\ \end{pmatrix} \hspace{5mm} ev_{\sigma 2} = \begin{pmatrix} 0.1282 \\ -0.9801 \\ 0.1512 \\ \end{pmatrix} \hspace{5mm} ev_{\sigma 3} = \begin{pmatrix} -0.6687 \\ 0.02718 \\ 0.7431 \\ \end{pmatrix}\]

Since the material is an isotropic linear elastic material, the eigenvectors of the stress and the strain matrix coincide. The maximum principal stress is 714.53 MPa; however, the von Mises stress is equal to 1031.19 MPa, which is higher than \sigma_{max}=750 MPa. This indicates that the material cannot sustain the stress matrix above, and in fact, the strain state above is impossible to occur before yielding or failure of the steel metal.

View Mathematica Code

str={{0.002,0.001,0.002},{0.001,-0.003,0.001},{0.002,0.001,0.0015}};
a=Eigensystem[str];
a//MatrixForm
str//MatrixForm
Ee=210000;
Nu=0.3;
G=Ee/2/(1+Nu);
Strainvector={str[[1,1]],str[[2,2]],str[[3,3]],2*str[[1,2]],2*str[[1,3]],2*str[[2,3]]};
Cc={{1/Ee,-Nu/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,1/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,-Nu/Ee,1/Ee,0,0,0},{0,0,0,1/G,0,0},
{0,0,0,0,1/G,0},{0,0,0,0,0,1/G}};
Dd=FullSimplify[Inverse[Cc]];
stressvector=Dd.Strainvector
Stressmatrix={{stressvector[[1]],stressvector[[4]],stressvector[[5]]},{stressvector[[4]],stressvector[[2]],stressvector[[6]]},
{stressvector[[5]],stressvector[[6]],stressvector[[3]]}};
Stressmatrix//MatrixForm
b=Eigensystem[Stressmatrix];
b//MatrixForm
VonMises[sigma_]:=Sqrt[1/2*((sigma[[1,1]]-sigma[[2,2]])^2+(sigma[[2,2]]-sigma[[3,3]])^2+(sigma[[3,3]]-
sigma[[1,1]])^2+6*(sigma[[1,2]]^2+sigma[[1,3]]^2+sigma[[2,3]]^2))];
VonMises[Stressmatrix]

View Python Code

import sympy as sp
import numpy as np
from sympy import *
from numpy.linalg import *
from scipy.linalg import *
sp.init_printing(use_latex = "mathjax")
def vonMises(M):
    return sp.sqrt(1/2*((M[0,0] - M[1,1])**2+(M[1,1] - M[2,2])**2+
                       (M[2,2] - M[0,0])**2+6*(M[0,1]**2+M[0,2]**2+M[1,2]**2)))
e = np.array([[0.002, 0.001, 0.002], [0.001, -0.003, 0.001], [0.002, 0.001, 0.0015]])
display("\u03b5 =", e)
eigen = eig(e) # eigen values and vectors
print("eigenvectors and eigenvalues\n",eigen)
# first index is eigen values
print("Eigenvalues (principal strains)= \n",eigen[0])
# Columns are normalized eigenvectors. Transpose is used to convert them into rows
print("Eigenvectors (principal directions) = \n",eigen[1].T)
Ee = 210000
Nu = 0.3
G = Ee/2/(1+Nu)
strainvector = Matrix([e[0,0], e[1,1], e[2,2], 2*e[0,1], 2*e[0,2], 2*e[1,2]])
display("strainvector:", strainvector)
CC = Matrix([[1/Ee, -Nu/Ee, -Nu/Ee, 0, 0, 0],
             [-Nu/Ee, 1/Ee, -Nu/Ee, 0, 0, 0], 
            [-Nu/Ee, -Nu/Ee, 1/Ee, 0, 0, 0],
            [0, 0, 0, 1/G, 0, 0], 
           [0, 0, 0, 0, 1/G, 0],
             [0, 0, 0, 0, 0, 1/G]])
display("CC: ",CC)
Dd = CC.inv()
display("D: ", Dd)
stressvector = Dd * strainvector
display("stressvector: ", stressvector)
stressmatrix = Matrix([[stressvector[0], stressvector[3], stressvector[4]], 
                       [stressvector[3], stressvector[1], stressvector[5]],
                       [stressvector[4], stressvector[5], stressvector[2]]])
display("stressmatrix: ", stressmatrix)
#convert the stress matrix to a numpy object for faster calculations of the eigenvalues and eigenvectors.
stressmatrix=np.array(stressmatrix).astype(np.float64)
eigen = eig(stressmatrix) # eigen values and vectors
print("eigenvectors and eigenvalues\n",eigen)
# first index is eigen values
print("Eigenvalues (principal stresses)= \n",eigen[0])
# Columns are normalized eigenvectors. Transpose is used to convert them into rows
print("Eigenvectors (principal directions) = \n",eigen[1].T)
vonMises(stressmatrix)
display("vonMises stress: ", vonMises(stressmatrix))

Example 2

In a confined compression test of a linear elastic isotropic material, the confining stress was -500MPa while the vertical stress was -100 MPa. The axial strains measured along the directions of the vectors e_1, e_2, and e_3 were -0.003, 0.001, and -0.003, respectively. Find Young’s modulus and Poisson’s ratio.

SOLUTION:

Using the constitutive relationships for linear elastic isotropic materials:

    \[\varepsilon_{11} = -\frac{500}{E} + \nu\frac{100}{E} + \nu\frac{500}{E} = -0.003 \Rightarrow -500 + 600\nu = -0.003E\]

    \[\varepsilon_{22} = \nu\frac{500}{E} - \frac{100}{E} + \nu\frac{500}{E} = 0.001 \Rightarrow -100 + 1000\nu = 0.001E\]

Solving the two equations:

    \[\nu = 0.222 \hspace{5mm} E = 122,222.00 MPa\]

View Python Code

from sympy import *
import sympy as sp
sp.init_printing(use_latex = "mathjax")
E, e1, e2, v, s1, s2, s3 = symbols("E \u03b5_{11} \u03b5_{22} \u03bd sigma_1 sigma_2 sigma_3")
Eq1 = Eq(e1, s1/E - v*s2/E -v*s3/E)
display(Eq1)
Eq1 = Eq1.subs({s1:-500, s2:-100, s3:-500, e1:-0.003})
display(Eq1)
Eq2 = Eq(e2, -v*s1/E + s2/E - v*s3/E)
display(Eq2)
Eq2 = Eq2.subs({s1:-500, s2:-100, s3:-500, e2:0.001})
display(Eq2)
Sol = solve((Eq1, Eq2), (E,v))
display("E = ", Sol[E])
display("\u03bd = ", Sol[v])

Example 3

Consider a transversely isotropic linear elastic piece of material with the isotropy plane being the plane of e_1 and e_2. The values of Young’s modulus and Poisson’s ratio in the plane of isotropy are equal to 100 MPa and 0.2, respectively. In a uniaxial tension test, when the material was stretched in the direction of e_1 by a stress of a value 1 MPa the strain measured in the direction of the vector e_3 was recorded to be –0.0005. In another uniaxial tension test, when the material was stretched in the direction of e_3 by a stress of a value 1 MPa, the strain measured in the direction of e_3 was 0.002. Find an estimate for the material constants E_{33}, \nu_{13}, and \nu_{31}.

SOLUTION:

The Young’s modulus in the direction of e_3 can be obtained from the second experiment:

    \[\varepsilon_{33} = \frac{\sigma_{33}}{E_{33}} \Rightarrow E_{33} = \frac{1}{0.002} = 500 MPa\]

From the first experiment, the value of Poisson’s ratio \nu_{13} can be obtained as follows:

    \[\varepsilon_{33} = -\nu_{13}\frac{\sigma_{11}}{E_{11}} \Rightarrow \nu_{13} = \frac{0.0005 \times 100}{1} = 0.05\]

The value of Poisson’s ratio \nu_{31} can be obtained from the symmetry of the coefficients matrix:

    \[\frac{\nu_{13}}{E_{11}} = \frac{\nu_{31}}{E_{33}} \Rightarrow = \nu_{31} = \frac{0.05 \times 500}{100} = 0.25\]

View Python Code

from sympy import *
import sympy as sp
sp.init_printing(use_latex = "mathjax")
E33, E11, e33, v13, v31, s11, s33 = symbols("E_{33} E_{11} \u03b5_{33} \u03bd_{13} \u03bd_{31} \u03c3_{11} \u03c3_{33}")
Eq1 = Eq(e33, s33/E33)
display(Eq1)
Eq1 = Eq1.subs({e33:0.002, s33:1})
display(Eq1)
sol = solve(Eq1, E33)
display(Eq(E33, sol[0]))
Eq2 = Eq(e33, -v13*s11/E11)
display(Eq2)
Eq2 = Eq2.subs({e33:-0.0005, s11:1, E11:100})
display(Eq2)
sol = solve(Eq2, v13)
display(Eq(v13, sol[0]))
Eq3 = Eq(v13/E11, v31/E33)
display(Eq3)
Eq3 = Eq3.subs({v13:0.05, E33:500, E11:100})
display(Eq3)
sol = solve(Eq3, v31)
display(Eq(v31, sol[0]))

Example 4

Consider an orthotropic material whose major axes of orthotropy lie in the directions of the vectors:

    \[e_1' = \begin{pmatrix} \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \\ 0 \\ \end{pmatrix} \hspace{5mm} e_2' = \begin{pmatrix} -\frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \\ 0 \\ \end{pmatrix} \hspace{5mm} e_3' = \begin{pmatrix} 0 \\ 0 \\ 1 \\ \end{pmatrix}\]

The material properties in the coordinate systems defined by the orthotropy axes are: E_{11}=100MPa, E_{22}=200MPa, E_{33}=350MPa, \nu_{31}=0.2, \nu_{32}=0.1, \nu_{21}=0.05, G_{12}=50MPa, G_{13}=70MPa, and G_{23}=90MPa. Consider the coordinate system defined by the basis set B=\{e_1,e_2,e_3\} where

    \[e_1 = \begin{pmatrix} 1 \\ 0 \\ 0 \\ \end{pmatrix} \hspace{5mm} e_2 = \begin{pmatrix} 0 \\ 1 \\ 0 \\ \end{pmatrix} \hspace{5mm} e_3 = \begin{pmatrix} 0 \\ 0 \\ 1 \\ \end{pmatrix}\]

If the material is in a state of uniaxial stress such that the only nonzero stress component is \sigma_{11}=1MPa, find the engineering strain matrix in the coordinate systems defined by the basis sets B and B'=\{e'_1,e'_2, e'_3\}.

Solution

Since the material is orthotropic, the coefficients matrix C depends on the coordinate system involved. The values of the material coefficients are given in the coordinate system defined by B'. Therefore, it is best to apply the constitutive relationship in the coordinate system of B'.

The stress matrix in the coordinate system of B is given by:

    \[\sigma = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ \end{pmatrix}\]

The rotation matrix required for the coordinate transformation is:

    \[Q = \begin{pmatrix} \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & 0 \\ -\frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & 0 \\ 0 & 0 & 1 \\ \end{pmatrix}\]

The stress matrix in the coordinate system defined by B' is given by:

    \[\sigma ' = Q\sigma Q^T = \begin{pmatrix} \frac{1}{2} & \frac{1}{2} & 0 \\ -\frac{1}{2} & \frac{1}{2} & 0 \\ 0 & 0 & 0 \\ \end{pmatrix}\]

The relationship between the vector representations of the stress and the strain in the coordinate system of B' can be used to find the strains in the coordinate system of B' (See Equation 2) as follows:

    \[\begin{pmatrix} \varepsilon_{11}' \\ \varepsilon_{22}' \\ \varepsilon_{33}' \\ 2\varepsilon_{12}' \\ 2\varepsilon_{13}' \\ 2\varepsilon_{23}' \\ \end{pmatrix} = \begin{pmatrix} \frac{1}{100} & -\frac{0.05}{200} & -\frac{0.2}{250} & 0 & 0 & 0 \\ -\frac{0.05}{200} & \frac{1}{200} & -\frac{0.1}{350} & 0 & 0 & 0 \\ -\frac{0.2}{350} & -\frac{0.1}{350} & \frac{1}{350} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{50} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{70} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{90} \\ \end{pmatrix} \begin{pmatrix} \sigma_{11}' \\ \sigma_{22}' \\ \sigma_{33}' \\ \sigma_{12}' \\ \sigma_{13}' \\ \sigma_{23}' \\ \end{pmatrix} = \begin{pmatrix} 0.00488 \\ 0.00238 \\ -0.00043 \\ -0.01 \\ 0 \\ 0 \\ \end{pmatrix}\]

Therefore, the strain matrix in the coordinate system of B' is given by:

    \[\varepsilon ' = \begin{pmatrix} 0.00488 & -0.005 & 0 \\ -0.005 & 0.00238 & 0 \\ 0 & 0 & -0.00043 \\ \end{pmatrix}\]

On the other hand, the strain matrix \varepsilon in the coordinate system of B is given by:

    \[\varepsilon ' = Q\varepsilon Q^T \Rightarrow \varepsilon = Q^T \varepsilon ' Q = \begin{pmatrix} 0.00863 & 0.00125 & 0 \\ 0.00125 & -0.001375& 0 \\ 0 & 0 & -0.00043 \\ \end{pmatrix}\]

It is very important to notice that, while a uniaxial state of stress was applied in the coordinate system defined by B, shear strain was developed in the same coordinate system. This is because the uniaxial state of stress was applied in a coordinate system different from the orthotropy axes coordinate system!

View Mathematica Code

E11=100;E22=200;E33=350;G12=50;G13=70;G23=90;
nu31=0.2;nu32=0.1;nu21=0.05;
nu13=nu31/E33*E11;nu12=nu21/E22*E11;nu23=nu32/E33*E22;
Ss={{1/E11,-nu21/E22,-nu31/E33,0,0,0} ,{-nu12/E11,1/E22,-nu32/E33,0,0,0} ,
{-nu13/E11,-nu23/E22,1/E22,0,0,0},{0,0,0,1/G12,0,0},{0,0,0,0,1/G13,0}, {0,0,0,0,0,1/G23}};
s={{1,0,0},{0,0,0},{0,0,0}};
Q=RotationMatrix[-45Degree,{0,0,1}]
sdash=Q.s.Transpose[Q];
sdash//MatrixForm
mat=sdash
sigmavector={mat[[1,1]],mat[[2,2]],mat[[3,3]],mat[[1,2]],mat[[1,3]],mat[[2,3]]}
sigmavector//MatrixForm
strvector=Ss.sigmavector;
strvector//MatrixForm
strdash={{strvector[[1]],strvector[[4]]/2,strvector[[5]]/2},{strvector[[4]]/2,strvector[[2]],strvector[[6]]/2},
{strvector[[5]]/2,strvector[[6]]/2,strvector[[3]]}} ;
strdash//MatrixForm
str=Transpose[Q].strdash.Q;
str//MatrixForm

View Python Code

from sympy import *
import sympy as sp
sp.init_printing(use_latex = "mathjax")
E11, E22, E33, G12, G13, G23, nu31, nu32, nu21 = 100, 200, 350, 50, 70, 90, 0.2, 0.1, 0.05
nu13 = nu31/E33*E11
nu12 = nu21/E22*E11
nu23 = nu32/E33*E22
Ss = Matrix([[1/E11, -nu21/E22, -nu31/E33, 0, 0, 0],
             [-nu12/E11, 1/E22, -nu32/E33, 0, 0, 0],
             [-nu13/E11, -nu23/E22, 1/E22, 0, 0, 0],
             [0, 0, 0, 1/G12, 0, 0],
             [0, 0, 0, 0, 1/G13, 0],
             [0, 0, 0, 0, 0, 1/G23]])
s = Matrix([[1,0,0], [0,0,0], [0,0,0]])
display("stress matrix in B: ", s)
#this rotates clockwise, so use - of the angle
Q = rot_axis3(pi/4)
display("rotation matrix required: ", Q)
sdash = Q*s*Transpose(Q)
display("stress matrix in B': ", sdash)
sigmavector = Matrix([[sdash[0,0]], [sdash[1,1]], [sdash[2,2]], [sdash[0,1]], [sdash[0,2]], [sdash[1,2]]])
strvector = Ss*sigmavector
strdash = Matrix([[strvector[0], strvector[3]/2, strvector[4]/2], 
                  [strvector[3]/2, strvector[1], strvector[4]/2], 
                  [strvector[4]/2, strvector[5]/2, strvector[2]]])
display("strain matrix in  B': ", strdash)
strain = Transpose(Q) * strdash * Q
display("strain matrix in B: ", strain)

Example 5

The in-plane strain components \varepsilon_{11}, \varepsilon_{22}, and \gamma_{12} of a specimen are 0.001, -0.002, and 0.003, respectively. If Young’s modulus is 10,000MPa and Poisson’s ratio is 0.2, find the stress and strain matrices if the material is in the state of plane strain and if the material is in the state of plane stress.

SOLUTION:

For the state of plane strain, the strain matrix has the following form:

    \[\varepsilon = \begin{pmatrix} 0.001 & 0.0015 & 0 \\ 0.0015 & -0.002 & 0 \\ 0 & 0 & 0 \\ \end{pmatrix}\]

Equation 11 can be used to find the stress components \sigma_{11}, \sigma_{22}, and \sigma_{12} as follows:

    \[\begin{pmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \\ \end{pmatrix} = \begin{pmatrix} 11111.1 & 2777.78 & 0 \\ 2777.78 & 11111.1 & 0 \\ 0 & 0 & 4166.67 \\ \end{pmatrix} \begin{pmatrix} 0.001 \\ -0.002 \\ 0.003 \\ \end{pmatrix} = \begin{pmatrix} 5.56 \\ -19.44 \\ 12.5 \\ \end{pmatrix} MPa\]

In a plane strain condition, the movement in the direction of e_3 is restrained, and thus, a restraining stress component \sigma_{33} develops according to the relationship:

    \[\sigma_{33} = \nu(\sigma_{11} + \sigma_{22}) = 0.2(5.56-19.4) = -2.78 MPa\]

Thus, the stress matrix has the following form:

    \[\sigma = \begin{pmatrix} 5.56 & 12.5 & 0 \\ 12.5 & -19.4 & 0 \\ 0 & 0 & -2.78 \\ \end{pmatrix}\]

View Mathematica Code for the State of Plane Strain:

Ee=10000;Nu=0.2;
Cc=Ee/(1+Nu)*{{(1-Nu)/(1-2Nu),Nu/(1-2Nu),0} ,{Nu/(1-2Nu),(1-Nu)/(1-2Nu),0},{0,0,1/2}};
Cc//MatrixForm
strvector={0.001,-0.002,0.003};
stressvector=Cc.strvector;
stressvector//MatrixForm
s=Table[0,{i,1,3},{j,1,3}];
s[[1,1]]=stressvector[[1]];
s[[2,2]]=stressvector[[2]];
s[[1,2]]=stressvector[[3]];
s[[2,1]]=s[[1,2]];
s[[3,3]]=Nu (s[[1,1]]+s[[2,2]]);
s//MatrixForm

View Python Code

from sympy import *
import sympy as sp
sp.init_printing(use_latex = "mathjax")
Ee, Nu = 10000, 0.2
Cc = Matrix([[(1 - Nu)/(1 - 2*Nu), Nu/(1 - 2*Nu), 0], 
             [Nu/(1 - 2*Nu), (1 - Nu)/(1 - 2*Nu), 0],
             [0, 0, 1/2]])
scalar = Ee/(1+Nu) 
Cc = scalar * Cc
strvector  = Matrix([0.001, -0.002, 0.003])
stressvector = Cc * strvector
s = Matrix([[stressvector[0], stressvector[2], 0], 
            [stressvector[2], stressvector[1], 0],
            [0, 0, 0]])
s[2,2] = Nu*(s[0,0] + s[1,1])
display("stress matrix: ", s)

In a plane stress condition, a nonzero strain component \varepsilon_{33} can develop, and thus, the strain matrix has the form:

    \[\varepsilon = \begin{pmatrix} 0.001 & 0.0015 & 0 \\ 0.0015 & -0.002 & 0 \\ 0 & 0 & \varepsilon_{33} \\ \end{pmatrix}\]

The relationship between the stresses and the strains in a plane stress condition (Equation 10) can be utilized to find the values of the non-zero stress components, as follows:

    \[\begin{pmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \\ \end{pmatrix} = \begin{pmatrix} 10416.7 & 2083.33 & 0 \\ 2083.33 & 10416.7 & 0 \\ 0 & 0 & 4166.67 \\ \end{pmatrix} \begin{pmatrix} 0.001 \\ -0.002 \\ 0.003 \\ \end{pmatrix} = \begin{pmatrix} 6.25 \\ -18.75 \\ 12.5 \\ \end{pmatrix} MPa\]

Thus, the stress matrix has the following form:

    \[\sigma = \begin{pmatrix} 6.25 & 12.5 & 0 \\ 12.5 & -18.75 & 0 \\ 0 &; 0 & 0 \\ \end{pmatrix}\]

In a plane stress condition, \varepsilon_{33} develops according to the relationship:

    \[\varepsilon_{33} = -\frac{\nu}{E}(\sigma_{11} + \sigma_{22}) = -\frac{0.2}{10000}(6.25-18.75) = 0.00025\]

View Mathematica Code for the State of Plane Strain:

Ee=10000;Nu=0.2;
Cc=Ee/(1+Nu)*{{1/(1-Nu),Nu/(1-Nu),0},{Nu/(1-Nu),1/(1-Nu),0},{0,0,1/2}};
Cc//MatrixForm
strvector={0.001,-0.002,0.003};
stressvector=Cc.strvector;
stressvector//MatrixForm
s=Table[0,{i,1,3},{j,1,3}];
s[[1,1]]=stressvector[[1]];
s[[2,2]]=stressvector[[2]];
s[[1,2]]=stressvector[[3]];
s[[2,1]]=s[[1,2]];
s[[3,3]]=0;
s//MatrixForm
str=Table[0,{i,1,3} ,{j,1,3}];
str[[1,1]]=strvector[[1]];
str[[2,2]]=strvector[[2]];
str[[1,2]]=str[[2,1]]=strvector[[3]]/2;
str[[3,3]]=-Nu/Ee*(s[[1,1]]+s[[2,2]]);
str//MatrixForm

View Python Code

from sympy import *
import sympy as sp
sp.init_printing(use_latex = "mathjax")
e33 = symbols("\u03b5_{33}")
Ee, Nu = 10000, 0.2
Cc = Matrix([[1/(1 - Nu), Nu/(1 - Nu), 0], 
             [Nu/(1 - Nu), 1/(1 - Nu), 0],
             [0, 0, 1/2]])
scaler = Ee/(1+Nu) 
Cc = scaler * Cc
strvector  = Matrix([0.001, -0.002, 0.003])
stressvector = Cc * strvector
s = Matrix([[stressvector[0], stressvector[2], 0], 
            [stressvector[2], stressvector[1], 0],
            [0, 0, 0]])
display("stress matrix: ", s)
strain = Matrix([[strvector[0], strvector[2]/2, 0], 
              [strvector[2]/2, strvector[1], 0],
              [0, 0, -Nu/Ee * (s[0,0] + s[1,1])]])
display("strain matrix: ", strain)
display(Eq(e33, strain[2,2]))

Example 6

The position function in dimensions of m of a plate adheres to the following form:

    \[x_1 = X_1 + 0.001X_2 + 0.0002X_1^2\]

    \[x_2 = 1.001X_2 - 0.0002X_2^2\]

    \[x_3 = X_3\]

If the plate has dimensions 2m in the direction of e_1, and 1m in the direction of e_2, with the bottom left corner coinciding with the origin of the coordinate system, then:

  • Show that this is a state of plane strain.
  • Find the engineering small strain matrix, and the Green strain matrix and show that the engineering small strain is adequate to describe the motion.
  • Draw the vector plot of the displacement function.
  • If Young’s modulus and Poisson’s ratio are equal to 210 GPa and 0.3, respectively, then find the equilibrium body forces to which this plate is subjected.
  • Draw the contour plot of the von Mises stress and the normal out-of-plane stress component.
SOLUTION:

Unlike the previous examples, this example shows a nonlinear displacement function resulting in a strain field that varies according to position. I.e., every material point will have a different local strain matrix. The displacement function for this plate has the following form:

    \[u = x - X = \begin{pmatrix} 0.0002X_1^2 + 0.001X_2 \\ 0.001X_2 - 0.0002X_2^2 \\ 0 \\ \end{pmatrix}\]

The following figures shows the vector plot of the displacement function as drawn by Mathematica:

The gradient of the displacement tensor is:

    \[\nabla u = \begin{pmatrix} 0.0004X_1 & 0.001 & 0 \\ 0 & 0.001-0.0004X_2 & 0 \\ 0 & 0 & 0 \\ \end{pmatrix}\]

The engineering and the Green strain matrices are:

    \[\varepsilon_{eng} = \frac{1}{2}(\nabla u + \nabla u^T) = \begin{pmatrix} 0.0004X_1 & 0.0005 & 0 \\ 0.0005& 0.001-0.0004X_2 & 0 \\ 0 & 0 & 0 \\ \end{pmatrix}\]

    \[\varepsilon_{Green} = \frac{1}{2}(\nabla u + \nabla u^T + \nabla u^T \nabla u)\]

    \[=\begin{pmatrix} 0.0004X_1+8 \times 10^{-8}X_1^2 & 0.0005+2 \times 10^{-7}X_1 & 0 \\ 0.0005 + 2 \times 10^{-7}X_1 & 0.001-0.0004X_2 + 8 \times 10^{-8}X_2^2 & 0 \\ 0 & 0 & 0 \\ \end{pmatrix}\]

The difference between the engineering and the Green strain matrices is very small and can be ignored. Thus, the engineering strain can be used to adequately describe the motion (i.e., the displacement is not associated with large rotations). In addition, the displacement component u_3=0 and the other displacement components are independent of the coordinate X_3. Therefore, the strain components \varepsilon_{33}=\varepsilon_{13}=\varepsilon_{23}=0, i.e., a state of plane strain.

In order to calculate the stresses, the vector representation of the strain will be used:

    \[\varepsilon = \begin{pmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ 2\varepsilon_{12} \\ 2\varepsilon_{13} \\ 2\varepsilon_{23} \\ \end{pmatrix} = \begin{pmatrix} 0.0004X_1 \\ 0.001 - 0.0004X_2 \\ 0 \\ 0.001 \\ 0 \\ 0 \\ \end{pmatrix}\]

Using the constitutive equation for the linear elastic isotropic material (Equation 7), the vector representation of the stress is:

    \[\sigma = \begin{pmatrix} 121.154+113.077X_1-48.4615X_2 \\ 282.692 + 48.4615X_1-113.077X_2 \\ 121.154+48.4615X_1+48.4615X_2 \\ 80.7692 \\ 0 \\ 0 \\ \end{pmatrix} MPa\]

Since the displacement is small, the coordinates x_i can be used instead of the coordinates X_i to find the equilibrium body forces:

    \[pb_1 = -\sum_{j=1}^3 \frac{\partial \sigma_{j1}}{\partial x_j} = -\frac{\partial \sigma_{11}}{\partial x_1} - \frac{\partial \sigma_{21}}{\partial x_2} - \frac{\partial \sigma_{31}}{\partial x_3} = -113.077 MN/m^3\]

    \[pb_2 = -\sum_{j=1}^3 \frac{\partial \sigma_{j2}}{\partial x_j} = -\frac{\partial \sigma_{12}}{\partial x_1} - \frac{\partial \sigma_{22}}{\partial x_2} - \frac{\partial \sigma_{32}}{\partial x_3} = 113.077 MN/m^3\]

    \[pb_3 = -\sum_{j=1}^3 \frac{\partial \sigma_{j3}}{\partial x_j} = -\frac{\partial \sigma_{13}}{\partial x_1} - \frac{\partial \sigma_{23}}{\partial x_2} - \frac{\partial \sigma_{33}}{\partial x_3} = 0 MN/m^3\]

The von Mises Stress:

    \[\sigma_{vM} = \sqrt{45665.7+4175.15(-2.5+X_1)X_1-20875.7X_2+4175.15X_1X_2+4175.15X_2^2}\]

The contour plot of the von Mises stress shows that the areas with the maximum von Mises stresses are the right and left bottom corners of the plate. Compare with the vector plot of the displacement (shown above) and note that the maximum displacements are not necessarily correlated with the maximum strain or stress. Since the plate is under a state of plane strain, the restraining stress component \sigma_{33} is non-zero. The contour plot of \sigma_{33} is shown below as well (Units are in MPa).

View Mathematica Code

Clear[Ee,Nu]
x1=X1+0.001*X2+0.0002*X1^2;
x2=1.001*X2-0.0002*X2^2;
x3=X3;
X={X1,X2,X3};
x={x1,x2,x3};
u=x-X;
u//MatrixForm
uplane=Table[u[[i]],{i,1,2}]
VectorPlot[uplane,{X1,0,2},{X2,0,1},AspectRatio->Automatic,PlotLabel->”u”,VectorStyle->Black]
gradu=Table[D[u[[i]],X[[j]]],{i,1,3},{j,1,3}];
gradu//MatrixForm
strGreen=FullSimplify[1/2*(gradu+Transpose[gradu]+Transpose[gradu].gradu)];
Expand[strGreen]//MatrixForm
str=FullSimplify[1/2*(gradu+Transpose[gradu])];
str//MatrixForm
Strainvector={str[[1,1]],str[[2,2]],str[[3,3]],2*str[[1,2]],2*str[[1,3]],2*str[[2,3]]};
Strainvector//MatrixForm
Ee=210000;
Nu=0.3;
G=Ee/2/(1+Nu);
Cc={{1/Ee,-Nu/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,1/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,-Nu/Ee,1/Ee,0,0,0}
,{0,0,0,1/G,0,0},{0,0,0,0,1/G,0},{0,0,0,0,0,1/G}};
Dd=FullSimplify[Inverse[Cc]];
stressvector=FullSimplify[Chop[Dd.Strainvector]];
stressvector//MatrixForm
Stressmatrix={{stressvector[[1]],stressvector[[4]],stressvector[[5]]},{stressvector[[4]],stressvector[[2]],stressvector[[6]]},
{stressvector[[5]],stressvector[[6]],stressvector[[3]]}};
Stressmatrix=FullSimplify[Chop[Stressmatrix]];
Stressmatrix//MatrixForm
X={X1,X2,X3};
s=Stressmatrix;
-Sum[D[s[[i,1]],X[[i]]],{i,1,3}]
-Sum[D[s[[i,2]],X[[i]]],{i,1,3}]
-Sum[D[s[[i,3]],X[[i]]],{i,1,3}]
VonMises[sigma_]:=Sqrt[1/2*((sigma[[1,1]]-sigma[[2,2]])^2+(sigma[[2,2]]-sigma[[3,3]])^2+(sigma[[3,3]]-
sigma[[1,1]])^2+6*(sigma[[1,2]]^2+sigma[[1,3]]^2+sigma[[2,3]]^2))];
a=FullSimplify[Chop[VonMises[Stressmatrix]]]
ContourPlot[a,{X1,0,2},{X2,0,1},AspectRatio->Automatic,ContourLabels->All,PlotLabel->”vonMises”,ColorFunction->”GrayTones”]
ContourPlot[Stressmatrix[[3,3]],{X1,0,2},{X2,0,1},AspectRatio->Automatic,ContourLabels->All,PlotLabel->”s33″,ColorFunction->”GrayTones”]

View Python Code

from sympy import *
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
sp.init_printing(use_latex = "mathjax")
def vonMises(M):
    return sp.sqrt(1/2*((M[0,0] - M[1,1])**2+(M[1,1] - M[2,2])**2+
                       (M[2,2] - M[0,0])**2+6*(M[0,1]**2+M[0,2]**2+M[1,2]**2)))

def plotVector(f, limits, title):
    fx, fy = [lambdify((X1,X2),fi) for fi in f]
    x1, xn, y1, yn = limits
    dx, dy = 10/100*(xn-x1),10/100*(yn-y1)
    xrange = np.arange(x1,xn,dx)
    yrange = np.arange(y1,yn,dy)
    X, Y = np.meshgrid(xrange, yrange)
    dxplot, dyplot = fx(X,Y), fy(X,Y)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.quiver(X, Y, dxplot, dyplot)
    plt.title(title)

def plotContour(f, limits, title):
    x1, xn, y1, yn = limits
    dx, dy = 10/100*(xn-x1),10/100*(yn - y1)
    xrange = np.arange(x1,xn,dx)
    yrange = np.arange(y1,yn,dy)
    X, Y = np.meshgrid(xrange, yrange)
    lx, ly = len(xrange), len(yrange)
    F = lambdify((X1,X2),f)
    # if F is constant, then generates array
    # if F generates array, doesn't change anything
    Z = F(X,Y)*np.ones(lx*ly).reshape(lx, ly)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    cp = ax.contourf(X,Y,Z)
    fig.colorbar(cp)
    plt.title(title)

x, x1, x2, x3, X1, X2, X3, u, pb1, pb2, pb3= symbols("x x_1 x_2 x_3 X_1 X_2 X_3 u pb_1 pb_2 pb_3")
x1 = X1 + 0.001*X2 + 0.0002*X1**2
x2 = 1.001*X2 - 0.0002*X2**2
x3 = X3
display("displacement function: ")
X = Matrix([[X1], [X2], [X3]])
x = Matrix([[x1], [x2], [x3]])
u = x - X
display(u)
limits = [0,2,0,1]
uplane = u[:2]
display(uplane)
plotVector(uplane, limits, "u")
gradu = Matrix([[diff(ui,Xj) for Xj in X] for ui in u])
display("displacement tensor gradient", gradu)
strGreen = 1/2 * (gradu + Transpose(gradu) + Transpose(gradu) * gradu)
strain = 1/2 * (gradu + Transpose(gradu))
display("engineering strain matrix: ", strain)
display("Green strain matrix: ", strGreen)
strainvector = Matrix([[strain[0,0]], [strain[1,1]], [strain[2,2]], [2 * strain[0,1]], [2 * strain[0,2]                       ], [2 * strain[1, 2]]])
display("small strain vector: ", strainvector)
Ee = 210000
Nu = 0.3
G = Ee/2/(1+Nu)
Cc = Matrix([[1/Ee, -Nu/Ee, -Nu/Ee, 0, 0, 0],
             [-Nu/Ee, 1/Ee, -Nu/Ee, 0, 0, 0], 
             [-Nu/Ee, -Nu/Ee, 1/Ee, 0, 0, 0],
             [0, 0, 0, 1/G, 0, 0], 
             [0, 0, 0, 0, 1/G, 0],
             [0, 0, 0, 0, 0, 1/G]])
Dd = Inverse(Cc)
stressvector = Dd * strainvector
display("stressvector: ", stressvector)     
stressmatrix = Matrix([[stressvector[0], stressvector[3], stressvector[4]],
                       [stressvector[3], stressvector[1], stressvector[5]],
                       [stressvector[4], stressvector[5], stressvector[2]]])
display("stressmatrix: ", stressmatrix)
s = stressmatrix
sum1 = -sum(diff(s[i,0], X[i]) for i in range(3)) 
display(Eq(pb1, sum1))
sum2 = -sum(diff(s[i,1], X[i]) for i in range(3)) 
display(Eq(pb2, sum2))
sum3 = sum(diff(s[i,2], X[i]) for i in range(3)) 
display(Eq(pb3, sum3))
a = vonMises(stressmatrix)
display("von Mises stress: ", a)
plotContour(a,limits,"vonMises")
plotContour(s[2,2],limits, "s33")

Example 7

The in plane displacements of a plate of dimensions of 5m and 1m in the directions of e_1 and e_2, respectively, with the bottom left corner situated at the origin are according to the relationship:

    \[u_1 = 0.02X_1 + 0.02X_1X_2\]

    \[u_2 = 0.0015X_2X_1^2 + 0.02X_2\]

If the plate is in the state of plane stress, draw the contour plot of the change of thickness of the plate, assuming that Young’s modulus is 1000 MPa, Poisson’s ratio is 0.4, and the thickness of the plate is 1mm.

SOLUTION:

Since the plate is in a state of plane stress, \sigma_{33}= 0 while a strain compoenent \varepsilon_{33} develops according to the equation:

    \[\varepsilon_{33} = -\frac{\nu}{E}(\sigma_{11} + \sigma_{22})\]

The in plane strain components can be obtained from the in plane gradient of the displacement function, as follows:

    \[\nabla u = \begin{pmatrix} 0.02+0.02X_2 & 0.02X_1 \\ 0.0003X_1X_2 & 0.02+0.0015X_1^2 \end{pmatrix}\]

    \[\begin{pmatrix} \varepsilon_{11} & \varepsilon_{12} \\ \varepsilon_{21} & \varepsilon_{22} \\ \end{pmatrix} = \frac{1}{2}(\nabla u + \nabla u ^T) = \begin{pmatrix} 0.02+0.02X_2 & (0.01+0.0015X_2)X_1 \\ (0.01+0.0015X_2)X_1 & 0.02+0.0015X_1^2 \end{pmatrix}\]

The stresses can be obtained using the plane stress constitutive relation (Equation 10):

    \[\begin{pmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \\ \end{pmatrix} = \begin{pmatrix} 33.33+0.714X_1^2+23.81X_2 \\ 33.33+1.786X_1^2+9.524X_2 \\ X_1(7.143+1.07X_2) \\ \end{pmatrix} MPa\]

Thus, the strain component \varepsilon_{33} is:

    \[\varepsilon_{33} = \frac{\nu}{E}(\sigma_{11} + \sigma_{22}) = -0.027-0.001X_1^2-0.013X_2\]

The change in thickness in mm can be obtained using the formula:

    \[\Delta = \varepsilon_{33}t\]

The required contour plot is shown below:

View Mathematica Code

X={X1,X2,X3};
u={0.02X1+0.02X1*X2,0.0015X2*X1^2+0.02X2};
VectorPlot[u,{X1,0,5},{X2,0,1},AspectRatio->Automatic]
gradu=Table[D[u[[i]],X[[j]]],{i,1,2},{j,1,2}];
gradu//MatrixForm
esmall=1/2*(gradu+Transpose[gradu]);
esmall//MatrixForm
strvector={esmall[[1,1]],esmall[[2,2]],2*esmall[[1,2]]};
strvector//MatrixForm
Ee=1000;Nu=0.4;Cc=Ee/(1+Nu)*{{1/(1-Nu),Nu/(1-Nu),0},{Nu/(1-Nu),1/(1-Nu),0},{0,0,1/2}};
stressvector=Cc.strvector;
stressvector//MatrixForm
strain=Table[0,{i,1,3},{j,1,3}];
strain[[1;;2,1;;2]]=esmall;
strain[[3,3]]=-Nu/Ee*(stressvector[[1]]+stressvector[[2]]);
strain//MatrixForm
thickness=1;
delta=FullSimplify[strain[[3,3]]*thickness]
ContourPlot[delta,{X1,0,5},{X2,0,1},ContourLabels->All,AspectRatio->Automatic,PlotLabel->”delta Thickness mm”]

View Python Code

from sympy import *
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
sp.init_printing(use_latex = "mathjax")
E, e33, v, s11, s22, X1, X2, X3= symbols("E \u03b5_{33} \u03bd \u03c3_{11} \u03c3_{22} X_1 X_2 X_3")

def plotVector(f, limits, title):
    fx, fy = [lambdify((X1,X2),fi) for fi in f]
    x1, xn, y1, yn = limits
    dx, dy = 10/100*(xn-x1),10/100*(yn-y1)
    xrange = np.arange(x1,xn,dx)
    yrange = np.arange(y1,yn,dy)
    X, Y = np.meshgrid(xrange, yrange)
    dxplot, dyplot = fx(X,Y), fy(X,Y)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.quiver(X, Y, dxplot, dyplot)
    plt.title(title)

def plotContour(f, limits, title):
    x1, xn, y1, yn = limits
    dx, dy = 10/100*(xn-x1),10/100*(yn - y1)
    xrange = np.arange(x1,xn,dx)
    yrange = np.arange(y1,yn,dy)
    X, Y = np.meshgrid(xrange, yrange)
    lx, ly = len(xrange), len(yrange)
    F = lambdify((X1,X2),f)
    # if F is constant, then generates array
    # if F generates array, doesn't change anything
    Z = F(X,Y)*np.ones(lx*ly).reshape(lx, ly)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    cp = ax.contourf(X,Y,Z)
    fig.colorbar(cp)
    plt.title(title)

X = Matrix([[X1], [X2], [X3]])
u = Matrix([[0.02*X1 + 0.02*X1*X2], [0.0015*X2*X1**2 + 0.02*X2]])
limits = [0,5,0,1]
plotVector(u,limits, "u")
gradu = Matrix([[diff(ui,Xj) for Xj in X[:2]] for ui in u])
display("plane gradient for displacement function: ", gradu)
esmall = 1/2 * (gradu + Transpose(gradu))
strvector = Matrix([esmall[0,0], esmall[1,1], 2*esmall[0,1]])
Ee = 1000
Nu = 0.4
Cc = Matrix([[1/(1-Nu), Nu/(1-Nu), 0],
             [Nu/(1-Nu), 1/(1-Nu), 0],
             [0,0,1/2]])
Cc = Ee/(1+Nu) * Cc
stressvector = Cc * strvector
display("stressvector: ", stressvector)
strain = Matrix([[esmall[0,0], esmall[0,1], 0],
                 [esmall[1,0], esmall[1,1], 0],
                 [0, 0, -Nu/Ee * (stressvector[0] + stressvector[1])]])
thickness = 1
delta = strain[2,2] * thickness
display("strain component: ")
display(Eq(e33, v/E*(s11+s22)))
display(Eq(e33, delta))
plotContour(delta, limits, "delta thickness mm")

Example 8

A cylindrical pressure vessel is made out of a steel material. The material’s Young modulus and Poisson’s ratio are 210 GPa and 0.3 respectively. The pressure vessel is used to hold a gas with an internal pressure P of 3 MPa. The vessel’s average diameter is 2 m, and its thickness is 10mm. Find the von Mises stress, the change in diameter, and the change in thickness of the cylinder in the following two conditions:

  • Condition 1: Capped ends.
  • Condition 2: Restrained ends.
SOLUTION:

As shown in the sketch, the orthonormal basis vectors are chosen in this example such that e_1 is aligned with the longitudinal axis of the pipe and e_2 is vertical. The analyzed part of the pressure wall is taken to be perpendicular to e_3.

For condition 1: The normal stress component in the longitudinal direction is:

    \[\sigma_{11} = \frac{P \pi r^2}{2 \pi rt} = \frac{Pr}{2t} = 150MPa\]

The normal stress component in the direction of the vector e_2 is:

    \[\sigma_{22} = \frac{Pr}{t} = 300MPa\]

All the remaining components of the stress matrix are zero. Thus, the stress matrix has the following form:

    \[\sigma = \begin{pmatrix} 150 & 0 & 0 \\ 0 & 300 & 0 \\ 0 & 0 & 0 \\ \end{pmatrix} MPa\]

The von Mises Stress is:

    \[\sigma_{vM} = 150\sqrt{3} = 259.81 MPa\]

The strain-stress relationship (Equation 6) can be used to find the vector representation of the strain:

    \[\begin{pmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ 2\varepsilon_{12} \\ 2\varepsilon_{13} \\ 2\varepsilon_{23} \\ \end{pmatrix} = \begin{pmatrix} 0.000286 \\ 0.001214 \\ -0.00064 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix}\]

The strain component \varepsilon_{22} gives the relative change of diameter as follows:

    \[\varepsilon_{22} = \frac{2\pi (r+\Delta r) - 2\pi r}{2 \pi r} = \frac{\Delta r}{r} = \frac{\Delta D}{D} = 0.001214\]

Therefore the change in diameter is:

    \[\Delta D = 0.001214D = 0.001214 \times 2 = 0.0024m = 2.4mm\]

The component \varepsilon_{33} provides the relative change in the thickness of the pressure vessel. Therefore:

    \[\Delta t = \varepsilon_{33}t = -0.00064(10) = -0.0064mm\]

View Mathematica Code for Condition 1:

P=3;r=1000;t=10;Ee=210000;Nu=0.3;
s=Table[0,{i,1,3},{j,1,3}];
s[[1,1]]=P*r/2/t
s[[2,2]]=P*r/t
VonMises[sigma_]:=Sqrt[1/2*((sigma[[1,1]]-sigma[[2,2]])^2+(sigma[[2,2]]-sigma[[3,3]])^2+(sigma[[3,3]]-
sigma[[1,1]])^2+6*(sigma[[1,2]]^2+sigma[[1,3]]^2+sigma[[2,3]]^2))];
a=VonMises[s]
G=Ee/2/(1+Nu);
Ss={{1/Ee,-Nu/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,1/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,-Nu/Ee,1/Ee,0,0,0},{0,0,0,1/G,0,0},
{0,0,0,0,1/G,0},{0,0,0,0,0,1/G}};
stressvector={s[[1,1]],s[[2,2]],s[[3,3]],s[[1,2]],s[[1,3]],s[[2,3]]}
strainvector=Ss.stressvector;
strainvector//MatrixForm
deltadiameter=strainvector[[2]]*2*r
deltathickness=strainvector[[3]]*t

View Python Code

from sympy import *
import sympy as sp
sp.init_printing(use_latex = "mathjax")
def vonMises(M):
    return sp.sqrt(1/2*((M[0,0] - M[1,1])**2+(M[1,1] - M[2,2])**2+
                       (M[2,2] - M[0,0])**2+6*(M[0,1]**2+M[0,2]**2+M[1,2]**2)))

P, r, t, Ee, Nu =  3, 1000, 10, 210000, 0.3
s = Matrix([[P*r/2/t, 0, 0], 
            [0, P*r/t, 0],
            [0, 0, 0]])
display("stress matrix: ", s)
a = vonMises(s)
display("von Mises Stress: ", a)
G = Ee/2/(1+Nu)
Ss = Matrix([[1/Ee ,-Nu/Ee, -Nu/Ee, 0, 0, 0], 
             [-Nu/Ee, 1/Ee, -Nu/Ee, 0, 0, 0],
             [-Nu/Ee, -Nu/Ee, 1/Ee, 0, 0, 0],
             [0, 0, 0, 1/G, 0, 0],
             [0, 0, 0, 0, 1/G, 0],
             [0, 0, 0, 0, 0, 1/G]])
stressvector = Matrix([[s[0,0]], [s[1,1]], [s[2,2]], [s[0,1]], [s[0,2]], [s[1,2]]])
strainvector = Ss * stressvector
display("strainvector: ", strainvector)
deltadiameter = strainvector[1] * 2 * r
display("change in diameter: ", deltadiameter)
deltathickness = strainvector[2] * t
display("change in thickness: ", deltathickness)

For condition 2, the normal stress in the longitudinal direction (\sigma_{11}) is not known since the vessel is restrained in the longitudinal direction. The circumferential stress (normal stress in the direction of the vector e_2) is given by:

    \[\sigma_{22} = \frac{Pr}{t} = 300 MPa\]

Thus, the stress matrix has the form:

    \[\sigma = \begin{pmatrix} \sigma_{11} & 0 & 0 \\ 0 & 300 & 0 \\ 0 & 0 & 0 \\ \end{pmatrix} MPa\]

Since for condition 2, the pipe is restrained from expanding or contracting, therefore, \varepsilon_{11}=0. Therefore:

    \[0 = \varepsilon_{11} = \frac{\sigma_{11}}{E} - \frac{\nu}{E}(\sigma_{22} + \sigma_{33}) \Rightarrow \sigma_{11} = 0.3 \times 300 = 90 MPa\]

The von Mises Stress is thus given by:

    \[\sigma_{vM} = 266.65 MPa\]

The strain-stress relationship (Equation 6) can be used to find the vector representation of the strain:

    \[\begin{pmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ 2\varepsilon_{12} \\ 2\varepsilon_{13} \\ 2\varepsilon_{23} \\ \end{pmatrix} = \begin{pmatrix} 0 \\ 0.0013 \\ -0.00056 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix}\]

The strain component \varepsilon_{22} gives the relative change of diameter as follows:

    \[\varepsilon_{22} = \frac{2 \pi (r+\Delta r) - 2\pi r}{2 \pi r} = \frac{\Delta r}{r} = \frac{\Delta D}{D} = 0.0013\]

Therefore the change in diameter is:

    \[\Delta D = 0.0013D = 0.0013 \times 2 = 0.0026m = 2.6m\]

The component \varepsilon_{33} provides the relative change in the thickness of the pressure vessel. Therefore:

    \[\Delta t = \varepsilon_{33}t = -0.00056(10) = -0.0056mm\]

View Mathematica Code for Condition 2:

P=3;r=1000;t=10;Ee=210000;Nu=0.3;
s=Table[0,{i,1,3},{j,1,3}];
s[[2,2]]=P*r/t
s[[1,1]]=Nu*s[[2,2]]
VonMises[sigma_]:=Sqrt[1/2*((sigma[[1,1]]-sigma[[2,2]])^2+(sigma[[2,2]]-
sigma[[3,3]])^2+(sigma[[3,3]]-sigma[[1,1]])^2+6*(sigma[[1,2]]^2+sigma[[1,3]]^2+sigma[[2,3]]^2))];
a=VonMises[s]
G=Ee/2/(1+Nu);
Ss={{1/Ee,-Nu/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,1/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,-Nu/Ee,1/Ee,0,0,0},{0,0,0,1/G,0,0},
{0,0,0,0,1/G,0},{0,0,0,0,0,1/G}};
stressvector={s[[1,1]],s[[2,2]],s[[3,3]],s[[1,2]],s[[1,3]],s[[2,3]]}
strainvector=Ss.stressvector;
strainvector//MatrixForm
deltadiameter=strainvector[[2]]*2*r
deltathickness=strainvector[[3]]*t

View Python Code

from sympy import *
import sympy as sp
sp.init_printing(use_latex = "mathjax")
def vonMises(M):
    return sp.sqrt(1/2*((M[0,0] - M[1,1])**2+(M[1,1] - M[2,2])**2+
                       (M[2,2] - M[0,0])**2+6*(M[0,1]**2+M[0,2]**2+M[1,2]**2)))

P, r, t, Ee, Nu =  3, 1000, 10, 210000, 0.3
s = Matrix([[0, 0, 0], 
            [0, P*r/t, 0],
            [0, 0, 0]])
s[0,0] = Nu * s[1,1]
display("stress matrix: ", s)
a = vonMises(s)
display("von Mises Stress: ", a)
G = Ee/2/(1+Nu)
Ss = Matrix([[1/Ee ,-Nu/Ee, -Nu/Ee, 0, 0, 0], 
             [-Nu/Ee, 1/Ee, -Nu/Ee, 0, 0, 0],
             [-Nu/Ee, -Nu/Ee, 1/Ee, 0, 0, 0],
             [0, 0, 0, 1/G, 0, 0],
             [0, 0, 0, 0, 1/G, 0],
             [0, 0, 0, 0, 0, 1/G]])
stressvector = Matrix([[s[0,0]], [s[1,1]], [s[2,2]], [s[0,1]], [s[0,2]], [s[1,2]]])
strainvector = Ss * stressvector
display("strainvector: ", strainvector)
deltadiameter = strainvector[1] * 2 * r
display("change in diameter: ", deltadiameter)
deltathickness = strainvector[2] * t
display("change in thickness: ", deltathickness)

Problems

  1. What does a negative Poisson’s ratio mean? Conduct a literature search and write a small paragraph supporting or opposing that materials with negative Poisson’s ratio exist (cite any references used).

  2. What does a Poisson’s ratio>0.5 mean? Conduct a literature search and write a small paragraph supporting or opposing that materials with Poisson’s ratio>0.5 exist (cite any references used).

  3. Assuming a linear elastic material such that the relationship between the stress and the strain is described by the following matrix:

    \begin{pmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ 2\varepsilon_{12} \\ 2\varepsilon_{13} \\ 2\varepsilon_{23} \\ \end{pmatrix} = \begin{pmatrix} \frac{1}{E} & -\frac{\nu}{E} & -\frac{\nu}{E} & 0 & 0 & 0 \\ -\frac{\nu}{E} & \frac{1}{E} & -\frac{\nu}{E} & 0 & 0 & 0 \\ -\frac{\nu}{E} & -\frac{\nu}{E} & \frac{1}{E} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{G} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{G} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{G} \\ \end{pmatrix} \begin{pmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{12} \\ \sigma_{13} \\ \sigma_{23} \\ \end{pmatrix}

    Assuming that the material is isotropic (i.e., the relationship above applies in any coordinate system) and starting from a state of stress described by the matrix:

    \begin{pmatrix} \sigma_{11} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ \end{pmatrix}

    Show that:

    G = \frac{E}{2(1+\nu)}

    (Hint: use a procedure similar to the one shown below the definition of the shear modulus).

  4. Show the following relationships for linear elastic isotropic materials:

    S_{ij} = 2\mu \varepsilon_{ij}'

    \sum_{k=1}^3 \frac{\sigma_{kk}}{3} = (\frac{2}{3} \mu + \lambda) \sum_{m=1}^3 \varepsilon_{mm}

    \varepsilon_{ij} = \frac{\sigma_{ij}}{2\mu} - \frac{\lambda}{2 \mu(2\mu+3\lambda)} \sum_{k=1}^3 \sigma_{kk} \delta_{ij}

    where: S_{ij} are the components of the deviatoric stress tensor
    \varepsilon'_{ij} are the components of the deviatoric strain tensor

    \varepsilon_{ij}' = \varepsilon_{ij} - \sum_{k=1}^3 \frac{\varepsilon_{kk}}{3} \delta_{ij}

    \lambda and \mu are Lamé’s constants.

  5. A state of stress of a linear isotropic material is represented by a positive definite symmetric stress matrix. Find the volumetric strain in terms of the principal stresses, Young’s modulus and Poisson’s ratio. Comment on whether the volume increases, decreases, stays constant or whether the information provided is not enough to answer this question.

  6. In a confined compression test of a linear elastic isotropic material, the testing apparatus was such that the confining pressure was equal to 50% of the applied vertical compressive stress P. Find the relationship between P and the volumetric strain in terms of Young’s modulus and Poisson’s ratio. Also, find the maximum shear stress, and the von Mises stress in terms of P.

  7. In a confined compression test of a linear elastic isotropic material, the testing apparatus was such that the specimen was restrained from expanding (or contracting) laterally. Find the relationship between the vertical compressive stress and the vertical strain in terms of Young’s modulus and Poisson’s ratio.

  8. Find the value of Poisson’s ratio at which both the shear modulus and the bulk modulus are equal to each other in value.

  9. A cylindrical sample is subjected to a confined compression test in which the horizontal stress is kept constant at –200 MPa, while the magnitude of the compressive vertical stress is increased to reach the value of –600 MPa. Find the volumetric strain if the shear and bulk moduli are 22.72 GPa and 20.833 GPa, respectively. (Answer: –1.6%)

  10. A linear elastic isotropic material with Young’s modulus of 210GPa and Poisson’s ratio of 0.3 is in a state of strain described by the infinitesimal strain matrix:

    \varepsilon = \begin{pmatrix} 0.03 & 0.002 & 0 \\ 0.002 & -0.015 & 0 \\ 0 & 0 & 0 \\ \end{pmatrix}

    • Find the principal strains and their directions.

    • Find the principal stresses and their directions.

    • Are the principal stresses and the principal strains aligned?

  11. Consider a linear elastic orthotropic material in a state of strain described by the infinitesimal strain matrix as in the previous problem. Assume the following material properties: E_{11}=200GPa, E_{22}=100GPa, E_{33}=50GPa, \nu_{12}=\nu_{13}=\nu_{23}=0.3, G_{12}=G_{13}=G_{23}=50GPa.

    • Find the principal stresses and their directions.

    • Are the principal stresses and the principal strains aligned?

  12. Using sound mathematical arguments show that the principal stresses and the principal strains are always aligned for linear elastic isotropic materials. (Hint: use the fact that Equations 6 and 7 are valid in any coordinate system)

  13. The in-plane displacement function of the two dimensional plate shown below is given by:

    u_1 = (-0.001X_1^2+0.005X_1)X_2

    u_2 = -0.0016x_1^3+0.0006x_1^2

    Assuming that the plate is of a linear elastic isotropic material with Young’s modulus E=210 GPa and Poisson’s ratio =0.2, draw the contour plot of the stress components \sigma_{11}, \sigma_{22}, \sigma_{12}, \sigma_{33}, \sigma_{vM}, and the hydrostatic stress for each of the following two cases:

    • the plate is in a state of plane stress.

    • the plate is in a state of plane strain.

  14. A cylindrical pressure vessel with an average diameter of 1m and a wall thickiness of 10mm is made of a material with Young’s modulus of 2.1 GPa and Poisson’s ratio 0.3. A constant pressure P of 1.3MPa is maintained inside the vessel

    • Assuming that the length of the pressure vessel is 5m, find the increase in length. (Answer: 31mm)

    • Find the increase in the diameter of the pressure vessel. (Answer: 26.3mm)

    • Find the change in thickness of the pressure vessel. (Answer: –0.14mm)

    • Find the volumetric strain throughout the pressure vessel wall. (Answer: 1.9%)

  15. An elastic cylindrical pressure vessel with average diameter d, thickness t, length L, Young’s modulus E, Poisson’s ratio \nu is subjected to an internal pressure P. Assume that the thickness is very small compared to the diameter. Find the circumferential stress, the longitudinal stress, the circumferential strain, and the longitudinal strain in the following two conditions:

    • Condition 1: The pressure vessel has two caps on its ends.

    • Condition 2: The pressure vessel is infinitely long and, thus, cannot expand or contract in the longitudinal direction.

    16. Leave a Reply

    Your email address will not be published.