Open Educational Resources

Plasticity: Examples and Problems

Examples

Example 1: Isotropic Hardening

The “true” stress versus “true” strain of a metallic material follows a bilinear relationship with true strains of 0.002 and 0.25 corresponding to stress levels of 400 and 550 MPa respectively. If Young’s modulus and Poisson’s ratio are equal to 200GPa and 0.3 respectively. Using von Mises incremental plasticity with isotropic hardening find the components of the elastic strain, the plastic strain and the total strain, the magnitude of the equivalent plastic strain parameter as a function of \sigma_{11}. Evaluate the total strains at \sigma_{11}=550MPa. Plot the relationship between \sigma_{11} and the total longitudinal strain component \varepsilon_{11} in each of the following conditions.

  • Condition 1: If a specimen is under a uniaxial state of stress.

  • Condition 2: If a specimen is under a stress \sigma_{11} and \varepsilon_{22} is kept at 0 in a plane stress condition.

  • Condition 3: If a specimen is under a stress \sigma_{11} and \varepsilon_{22} is kept at 0 in a plane strain condition.

Solution

First, the relationship between \overline{\varepsilon^p} and \sigma_Y is obtained from the information given. At \sigma_{11}=400 MPa, the elastic strain \varepsilon_{11}^e=400/200000=0.002 while the plastic strain \varepsilon_{11}^p=0. At \sigma_{11}=550 MPa, the elastic strain \varepsilon_{11}^e=550/200000=0.00275 while the plastic strain \varepsilon_{11}^p=0.24725. Thus, the relationship between \sigma_Y and the plastic strain parameter \overline{\varepsilon^p} is as follows:

    \[\sigma_Y = 400 + 150\frac{\overline{\varepsilon^p}}{0.24725}\]

The slope of this relationship is constant and can be calculated as follows:

    \[h = \frac{d\sigma_Y}{d \overline{\varepsilon^p}} = 606.673\]

Thus, the relationship between \sigma_Y and the plastic strain parameter \overline{\varepsilon^p} can be written as:

    \[\sigma_Y = 400 + h\overline{\varepsilon^p}\]

Condition 1: Uniaxial State of Stress

In a uniaxial state of stress condition, the hydrostatic stress, the von Mises stress, the Cauchy stress the deviatoric stress matrices have the following forms:

    \[p = \frac{\sigma_{11}}{3} \hspace{5mm} \sigma_{vM} = \sigma_{11} \hspace{5mm} \sigma = \begin{pmatrix} \sigma_{11} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ \end{pmatrix} \hspace{5mm} S = \begin{pmatrix} \frac{2 \sigma_{11}}{3} & 0 & 0 \\ 0 & -\frac{\sigma_{11}}{3} & 0 \\ 0 & 0 & -\frac{\sigma_{11}}{3} \\ \end{pmatrix}\]

The elastic strain components can be obtained using the elastic constitutive relationships as follows:

    \[\varepsilon_{11}^e = \frac{\sigma_{11}}{E} \hspace{5mm} \varepsilon_{22}^e = \varepsilon_{33}^3 = -\nu \frac{\sigma_{11}}{E} \hspace{5mm} \varepsilon_{12}^e = \varepsilon_{13}^e = \varepsilon_{23}^e = 0\]

The plastic strain components can be obtained as follows:

When \sigma_{vM}=\sigma_{11}\leq 400, \forall i,j: \varepsilon_{ij}^p=0 and the total strain components are equal to the elastic strain components.

When \sigma_{vM}=\sigma_{11}> 400, the increment in the plastic strain parameter \Delta\overline{\varepsilon^p} can be obtained from the consistency condition:

    \[d\phi = 0 \approx \Delta \phi = \frac{\partial \sigma_{vM}}{\partial \sigma_{ij}} \Delta \sigma_{ij} - \frac{d\sigma_Y}{d \overline{\varepsilon^p}} \Delta \overline{\varepsilon^p} = \frac{3S_{11}}{2\sigma_{vM}} \Delta \sigma_{11} - h \Delta \overline{\varepsilon^p} = \Delta \sigma_{11} - h \Delta \overline{\varepsilon^p}\]

Thus, an increment \Delta\sigma_{11} will produce an increment in the plastic strain parameter \Delta \overline{\varepsilon^p} that is equal to:

    \[\Delta \overline{\varepsilon^p} = \frac{\Delta \sigma_{11}}{h}\]

The magnitude of the plastic strain parameter can be obtained by summing (or integrating) the increments as follows:

    \[\overline{\varepsilon^p} = \int_{400}^{\sigma_{11}} \frac{1}{h} d\sigma_{11} = \frac{\sigma_{11} - 400}{h}\]

The increments in the plastic strain components can be obtained as follows:

    \[\Delta \varepsilon_{ij}^p = \Delta \overline{\varepsilon^p} \frac{\partial \phi}{\partial \sigma_{ij}} = \Delta \overline{\varepsilon^p} \frac{3S_{ij}}{2\sigma_{vM}}\]

Thus:

    \[\Delta \varepsilon_{11}^p = \Delta \overline{\varepsilon^p} \hspace{5mm} \Delta \varepsilon_{22}^p = \Delta \varepsilon_{33}^p = -\frac{\Delta \overline{\varepsilon^p}}{2} \hspace{5mm} \Delta \varepsilon_{12}^p = \Delta \varepsilon_{13}^p = \Delta \varepsilon_{23}^p = 0\]

The plastic strain components are equal to the sum of the respective increments and since h is constant for this particular material and there is no unloading:

    \[\varepsilon_{11}^p = \int_{400}^{\sigma_{11}} \frac{1}{h} d\sigma_{11} = \frac{\sigma_{11} - 400}{h} \hspace{5mm} \varepsilon_{22}^p = \varepsilon_{33}^p = -\frac{\varepsilon_{11}^p}{2}\]

The total strains at a stress level \sigma_{11}>400 are obtained as:

    \[\varepsilon_{11} = \varepsilon_{11}^e +\varepsilon_{11}^p = \frac{\sigma_{11}}{E} + \frac{\sigma_{11} - 400}{h}\]

    \[\varepsilon_{22} = \varepsilon_{33} = \varepsilon_{22}^e + \varepsilon_{22}^p = -\nu \frac{\sigma_{11}}{E} - \frac{\sigma_{11}-400}{2h}\]

    \[\varepsilon_{12} = \varepsilon_{13} = \varepsilon_{23} = 0\]

Notice that the relationship between \sigma_{11} and the total \varepsilon_{11} strain follows that given by the material curve:

    \[\varepsilon_{11} = \begin{cases} \frac{\sigma_{11}}{E} & \sigma_{11} \le 400 \\ \frac{\sigma_{11}}{E} + \frac{\sigma_{11}-400}{h} & \sigma_{11} > 400 \\ \end{cases} \Rightarrow \sigma_{11} = \begin{cases} E \varepsilon_{11} & \varepsilon_{11} \le 0.002 \\ \frac{E(400+h\varepsilon_{11})}{E+h} & \varepsilon_{11} > 0.002 \\ \end{cases}\]

Condition 2: Plane Stress Condition with \varepsilon_{22} = 0

The Cauchy stress and the deviatoric stress matrices have the following forms:

    \[\sigma = \begin{pmatrix} \sigma_{11} & 0 & 0 \\ 0 & \sigma_{22} & 0 \\ 0 & 0 & 0 \\ \end{pmatrix} \hspace{5mm} S = \begin{pmatrix} \frac{2\sigma_{11} - \sigma_{22}}{3} & 0 & 0 \\ 0 & \frac{2\sigma_{22} - \sigma_{11}}{3} & 0 \\ 0 & 0 & \frac{-\sigma_{11} - \sigma_{22}}{3} \\ \end{pmatrix}\]

In this case, the onset of initial yielding occurs when \sigma_{vM}=400 MPa, i.e., when:

    \[\sigma_{11}^2 - \sigma_{11}\sigma_{22} + \sigma_{22}^2 = (400)^2\]

Before developing any plastic strains, the total strain component \varepsilon_{22} is equal to the elastic strain component \varepsilon_{22}^e. The condition of loading ensures that \varepsilon_{22}=0 and thus:

    \[\varepsilon_{22}^e = \frac{\sigma_{22}}{E} - \frac{\nu \sigma_{11}}{E} = 0 \Rightarrow \sigma_{22} = \nu \sigma_{11}\]

Therefore, at the onset of initial yielding:

    \[\sigma_{11} = \frac{400}{\sqrt{1-\nu + \nu^2}} \approx 450MPa \hspace{5mm} \sigma_{22} = \frac{\nu 400}{\sqrt{1-\nu+ \nu^2}} \approx 135MPa\]

Therefore, when \sigma_{11}\leq 450MPa, the plastic strain components are zero while the elastic strain components are given by:

    \[\varepsilon_{11}^e = \frac{\sigma_{11}}{E} (1-\nu^2) \hspace{5mm} \varepsilon_{22}^e = 0 \hspace{5mm} \varepsilon_{33}^e = -\frac{\nu \sigma_{11}}{E} (1+\nu)\]

Once yielding occurs, the relationship between the stress components obtained above no longer follows. We seek to find expressions of the evolution of the elastic and plastic strain components as a function of the changes in the stress component. To do that, we first note that there are two nonzero stress components in the stress matrix, these are \sigma_{11} and \sigma_{22}. We also note that the condition of the total strain \varepsilon_{22}=0 can be used to relate the rate of changes of the stress components together. However, the total strain component \varepsilon_{22} has elastic and plastic components. The plastic component is a function of the plastic strain parameter \overline{\varepsilon^p}. So, the consistency condition is also used to relate the rates of the nonzero stress components with the rate of change of the plastic strain parameter. In other words, two equations are used to relate the rates of change of \sigma_{22} and of \overline{\varepsilon^p} to the rate of change of \sigma_{11}. The following shows the details of the procedure here.

First, the consistency condition can be used to relate the rate of change of the equivalent plastic strain parameter to the rate of change of the stress components:

    \[\dot{\phi} = \dot{\sigma}_{vM} - \dot{\sigma}_{Y} = 0 \Rightarrow \frac{\partial \sigma_{vM}}{\partial \sigma_{ij}} \dot{\sigma}_{ij} - \frac{d\sigma_Y}{d\overline{\varepsilon^p}} \dot{\overline{\varepsilon^p}} = 0 \Rightarrow \frac{3S_{11}}{2\sigma_{vM}} \dot{\sigma}_{11} + \frac{3S_{22}}{2\sigma_{vM}} \dot{\sigma}_{22} - h\overline{\varepsilon^p} = 0\]

Rearranging yields the relationship between the rate of change of the equivalent plastic strain parameter and the rate of change of the stress components:

    \[\dot{\overline{\varepsilon^p}} = \frac{(2\sigma_{11} - \sigma_{22}) \dot{\sigma}_{11} + (2\sigma_{22} - \sigma_{11}) \dot{\sigma}_{22}} {2h\sigma_{vM}}\]

The boundary conditions dictate that \dot{\varepsilon}_{22}=\dot{\varepsilon}_{22}^e+\dot{\varepsilon}_{22}^p=0, which leads to the following relationship between the rates of the equivalent plastic strain and the stress components:

    \[\frac{\dot{\sigma}_{22}}{E} - \frac{\nu \dot{\sigma}_{11}}{E} + \dot{\overline{\varepsilon^p}} \frac{\partial \phi}{\partial \sigma_{22}} = 0\]

Substituting for \dot{\overline{\varepsilon^p}} in the above equation yields \dot{\sigma}_{22} as an explicit function of \dot{\sigma}_{11}, \sigma_{11} and \sigma_{22}:

    \[\dot{\sigma}_{22} = \dot{\sigma}_{11} \frac{E(\sigma_{11} - 2\sigma_{22})(2\sigma_{11}-\sigma_{22}) + 4h\nu \sigma_{vM}^2}{E(\sigma_{11} - 2\sigma_{22})^2 + 4h\sigma_{vM}^2}\]

Substituting the value of \dot{\sigma}_{22} back into the expression of \dot{\overline{\varepsilon^p}} yields the following:

    \[\dot{\overline{\varepsilon^p}} = 2\dot{\sigma}_{11} \frac{((2 - \nu)\sigma_{11} + \sigma_{22}(2\nu - 1))\sigma_{vM}} {E(\sigma_{11} - 2\sigma_{22})^2 + 4h\sigma_{vM}^2}\]

The rate of change of the individual plastic strain components can then be calculated as a function of \dot{\overline{\varepsilon^p}} as follows:

    \[\dot{\varepsilon}_{ij}^p = \dot{\overline{\varepsilon^p}} \frac{\partial \phi}{\partial \sigma_{ij}} = \dot{\overline{\varepsilon^p}} \frac{3S_{ij}}{2\sigma_{vM}}\]

Substituting for \dot{\overline{\varepsilon^p}}, the following expressions are obtained:

    \[\dot{\varepsilon}_{11}^p = \dot{\sigma_{11}}(2\sigma_{11}-\sigma_{22}) \frac{\sigma_{11}(2-\nu) + \sigma_{22}(2\nu - 1)}{E(\sigma_{11}-2\sigma_{22})^2 + 4h\sigma_{vM}^2}\]

    \[\dot{\varepsilon}_{22}^p = \dot{\sigma_{11}}(2\sigma_{22}-\sigma_{11}) \frac{\sigma_{11}(2-\nu) + \sigma_{22}(2\nu - 1)}{E(\sigma_{11}-2\sigma_{22})^2 + 4h\sigma_{vM}^2}\]

    \[\dot{\varepsilon}_{33}^p = \dot{\sigma_{11}}(\sigma_{11}+\sigma_{22}) \frac{\sigma_{11}(\nu -2) + \sigma_{22}(1 - 2\nu)}{E(\sigma_{11}-2\sigma_{22})^2 + 4h\sigma_{vM}^2}\]

The problem reduces to finding a solution to the ordinary differential equations of the following form:

    \[\dot{\sigma}_{22} = f(\sigma_{11}, \sigma_{22}, \dot{\sigma}_{11}) \hspace{5mm} \dot{\overline{\varepsilon^p}} = g(\sigma_{11}, \sigma_{22}, \dot{\sigma}_{11}) \hspace{5mm} \dot{\varepsilon}_{ij}^p = h_{ij}(\sigma_{11}, \sigma_{22}, \dot{\sigma}_{11})\]

The initial conditions for these equations are given by \sigma_{11}=450 MPa, \sigma_{22}=135 MPa, \overline{\varepsilon^p}=\varepsilon_{ij}^p=0. The rate of change of \sigma_{11}, namely \sigma_{11} can be arbitrarily chosen. In order to solve the above ordinary differential equations, an independent time parameter t is introduced such that \sigma_{11}=t, therefore, \dot{\sigma}_{11}=1. The integration of the above equations can be performed numerically using a spreadsheet or using the NDSolve[] function in Mathematica with an arbitrary upper limit for \sigma_{11} or t, which in this example was taken as 650MPa. Note that if a spreadsheet is to be used, then the backward Euler integration method would produce more accurate results than the forward Euler integration method. The following figures show the plots of \sigma_{22} and the various nonzero plastic strain components as a function of \sigma_{11}. The plot of \sigma_{22} versus \sigma_{11} reveals three different regions during the loading process. In the initial elastic region, the stresses are related using Poisson’s ratio. In the second region, as plastic strain develops, the stress \sigma_{22} increases rapidly up to around half \sigma_{11}. In the third region, as \sigma_{11} increases, the stress \sigma_{22} is almost half that of \sigma_{11}. The plastic strain components reveal that as the specimen is pulled beyond its initial elasticity limit, positive \varepsilon_{11}^p strain is developed. In contract, negative \varepsilon_{22}^p and \varepsilon_{33}^p are developed. This indicates that if the stresses are removed, the specimen will be longer in the direction of the basis vector e_1 while shorter in the other two directions. The sum of the plastic strain components \varepsilon_{11}^p+\varepsilon_{22}^p+\varepsilon_{33}^p is always equal to zero.

The following figures show the plots of the various nonzero elastic strain components as a function of \sigma_{11}. Notice that the elastic strain \varepsilon_{22}^e is equal to zero when \sigma_{11} is less than the initial elastic limit. Beyond that, \varepsilon_{22}^e increases rapidly similar to the plastic strain component \varepsilon_{22}^p but with a negative sign to maintain a zero value for the total strain \varepsilon_{22}.

The following figures show the plots of the various nonzero total strain components as a function of \sigma_{11}:

View Mathematica Code


Clear[ds11,ds22,Svm,s11,s22,s33]
(*Values of material constants*)
epmax=25/100-550/200000;
hslope=N[150/epmax];
rule={Ee->200000,nu->0.3,h->hslope};
(*Stress, hydrostatic stress, deviatoric stress and von Mises stress*)
s={{s11,0,0},{0,s22,0},{0,0,s33}};
hydr=(s11+s22+s33)/3;
s33=0;
Sdev=FullSimplify[s-hydr*IdentityMatrix[3]];
Sdev//MatrixForm
Sigmavm=Sqrt[1/2*((s11-s22)^2+(s11-s33)^2+(s22-s33)^2)];
s=FullSimplify[(Sigmavm)^2]
(*The stresses at the elastic limit utilizing ee22=0*)
b=Solve[{Sqrt[s]==sys,s22==nu*s11},{s11,s22}]
(*The strains at the elastic limit*)
ee11=(s11/Ee-nu*s22/Ee)/.s22->nu*s11
ee22=(-nu*s11/Ee+s22/Ee)/.s22->nu*s11
ee33=(-nu*s11/Ee-nu*s22/Ee)/.s22->nu*s11
Sdev=Sdev/.{s11->s11[t],s22->s22[t]}
Sigmavm=Sigmavm/.{s11->s11[t],s22->s22[t]}
(*Consistency condition*)
cep'[t]=(3/2*Sdev[[1,1]]*s11'[t]/h/Svm+3/2*Sdev[[2,2]]*s22'[t]/h/Svm)
(*ep22 Plastic strain Rate *)
cep22'[t]=FullSimplify[cep'[t]*3/2*Sdev[[2,2]]/Svm]
(*ee22 elastic strain rate*)
Clear[ee22]
ee22'[t]=s22'[t]/Ee-nu*s11'[t]/Ee
(*Finding the relationship between s11'[t] and s22'[t]*)
a=FullSimplify[Solve[ee22'[t]+cep22'[t]==0,{s22'[t]}]]
(*Making everything function of s11'[t] alone*)
eq1=ep'[t]==FullSimplify[cep'[t]/.a[[1]]]
eq2=ep11'[t]==FullSimplify[(cep'[t]/.a[[1]])*3/2*Sdev[[1,1]]/Svm]
eq3=ep22'[t]==FullSimplify[(cep'[t]/.a[[1]])*3/2*Sdev[[2,2]]/Svm]
eq4=ep33'[t]==FullSimplify[(cep'[t]/.a[[1]])*3/2*Sdev[[3,3]]/Svm]
(*Finding s11 and s22 as a function of t using NDsolve, the procedure  assumes s11=t*)
s22equation=s22'[t]/.a[[1]]
s22equation=s22equation/.rule;
Svm=Sigmavm/.rule;
s11initial=400/Sqrt[1-nu+nu^2]/.rule;
ff=NDSolve[{(s22'[t]==s22equation),s11'[t]==1,s11[s11initial]==s11initial,
(s22[s11initial]==(nu/.rule)*s11initial)},{s22,s11},{t,s11initial,650}]
s22function=Piecewise[{{nu*t/.rule,0<=t<=s11initial},{s22[t]/.ff[[1]],t>s11initial}}]
(*Finding plastic strain components as a function of t using NDsolve, the procedure assumes s11=t*)
rule2={s11[t]->t,s22[t]->s22function,s22'[t]->D[s22function,t],s11'[t]->1};
crule=Flatten[{rule,rule2}]
eq1=eq1/.crule
eq2=eq2/.crule;
eq3=eq3/.crule;
eq4=eq4/.crule;
a=NDSolve[{eq1,ep[s11initial]==0},ep,{t,s11initial,650}]
b=NDSolve[{eq2,ep11[s11initial]==0},ep11,{t,s11initial,650}]
c=NDSolve[{eq3,ep22[s11initial]==0},ep22,{t,s11initial,650}]
d=NDSolve[{eq4,ep33[s11initial]==0},ep33,{t,s11initial,650}]
(*Plastic Strain functions of s11=t*)
epfunction=Piecewise[{{0,0<=t<=s11initial},{ep[t]/.a[[1]],t>s11initial}}]
ep11function=Piecewise[{{0,0<=t<=s11initial},{ep11[t]/.b[[1]],t>s11initial}}]
ep22function=Piecewise[{{0,0<=t<=s11initial},{ep22[t]/.c[[1]],t>s11initial}}]
ep33function=Piecewise[{{0,0<=t<=s11initial},{ep33[t]/.d[[1]],t>s11initial}}]
(*Elastic Strain functions of s11=t and s22*)
ee11function=(t/Ee-s22function/Ee*nu)/.rule;
ee22function=(-nu*t/Ee+s22function/Ee)/.rule;
ee33function=(-nu*t/Ee-nu*s22function/Ee)/.rule;
(*Producing the plots*)
Plot[s22function,{t,0,660},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","s22(MPa)"},PlotStyle->Black]
Plot[ep11function,{t,0,660},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","ep11(MPa)"},PlotStyle->Black]
Plot[ep22function,{t,0,660},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","ep22(MPa)"},PlotStyle->Black]
Plot[ep33function,{t,0,660},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","ep33(MPa)"},PlotStyle->Black]
Plot[epfunction,{t,0,660},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","ep(MPa)"},PlotStyle->Black]
Plot[ee11function,{t,0,660},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","ee11(MPa)"},PlotStyle->Black]
Plot[ee22function,{t,0,660},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","ee22(MPa)"},PlotStyle->Black]
Plot[ee33function,{t,0,660},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","ee33(MPa)"},PlotStyle->Black]
Plot[ee11function+ep11function,{t,0,660},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","e11(MPa)"},PlotStyle->Black]
Plot[Chop[ee22function+ep22function,10^-7],{t,0,660},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","e22(MPa)"},PlotStyle->Black]
Plot[ee33function+ep33function,{t,0,660},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","e33(MPa)"},PlotStyle->Black]

View Python Code

from sympy import *
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.interpolate import InterpolatedUnivariateSpline
sp.init_printing(use_latex = "mathjax")
s11, s22, nu, E, h, Svm, t = symbols("\u03c3_{11} \u03c3_{22} nu E h \u03c3_{vM} t")
epmax = 25/100-550/200000
hslope = 150/epmax
display(epmax)
display(hslope)
rule = {E:200000, nu:0.3, h:hslope}
display("Stress, hydrostatic stress, deviatoric stress and von Mises stress")
s33 = 0
s = Matrix([[s11,0,0],[0,s22,0],[0,0,s33]])
hydr = (s11+s22+s33)/3
display("Cauchy stress: ", s)
Sdev = s-hydr*eye(3)
display("Deviatoric stress: ", Sdev)
Sigmavn = sqrt(1/2*((s11-s22)**2+(s11-s33)**2+(s22-s33)**2))
s = simplify(Sigmavn**2)
display(s)
display("The stresses at the elastic limit utilizing ee22=0")
b = solve((sqrt(s) - 400, s22 - s11*nu), (s11,s22))
display(b)
display("The strains at the elastic limit")
ee11=(s11/E-nu*s22/E).subs(s22,nu*s11)
ee22=(-nu*s11/E+s22/E).subs(s22,nu*s11)
ee33=(-nu*s11/E-nu*s22/E).subs(s22,nu*s11)
display(ee11,ee22,ee33)
s11t = Function('\u03c3_{11}')
s22t = Function('\u03c3_{22}')
Sdev = Sdev.subs({s11:s11t(t), s22:s22t(t)})
Sigmavn = Sigmavn.subs({s11:s11t(t), s22:s22t(t)})
display("Consistency condition")
display(Sdev)
ceptp = (3/2*Sdev[0,0]*s11t(t).diff(t)/h/Svm+3/2*Sdev[1,1]*s22t(t).diff(t)/h/Svm)
display(ceptp)
display("ep22 Plastic strain Rate")
cep22tp = simplify(ceptp*3/2*Sdev[1,1]/Svm)
display(cep22tp)
display("ee22 Elastic strain Rate")
ee22tp = s22t(t).diff(t)/E-nu*s11t(t).diff(t)/E
display(ee22tp)
display("Finding the relationship between s11t and s22t")
a = solve(ee22tp + cep22tp, s22t(t).diff(t))
display(a)
display("Making everything a function of s11t alone")
eq1 = eptp = simplify(ceptp.subs(s22t(t).diff(t), a[0]))
eq2 = ep11tp = simplify(ceptp.subs(s22t(t).diff(t), a[0])*3/2*Sdev[0,0]/Svm)
eq3 = ep22tp = simplify(ceptp.subs(s22t(t).diff(t), a[0])*3/2*Sdev[1,1]/Svm)
eq4 = ep33tp = simplify(ceptp.subs(s22t(t).diff(t), a[0])*3/2*Sdev[2,2]/Svm)
display(eq1, eq2, eq3, eq4)
display("Finding s11 and s22 as a function of t using NDsolve, the procedure  assumes s11=t")
s22equation = a[0].subs(rule)
Svmsubs = Sigmavn.subs(rule)
s11initial = (400/sqrt(1-nu+nu**2)).subs(rule)
display(s22equation, Svmsubs, s11initial)
s22equation = s22equation.subs(Svm,Svmsubs)
s22equ = s22equation.subs({s11t(t):s11, s22t(t):s22, s11t(t).diff(t):1})
def rhs1(y,t):
    f = sp.lambdify([s11, s22], s22equ)
    f = f(y[0],y[1])
    ds22tdt = f
    dss11tdt = 1
    return [dss11tdt,ds22tdt]
s11initial = float(s11initial)
t1 = np.linspace(s11initial,650,20)
s22t0 = nu.subs(rule)*s11initial
s11t0 = s11initial
yzero = [s11t0,s22t0]
ff = odeint(rhs1, yzero, t1)
# display(ff)
# plot results
# aa is for plot1
lambda1 = lambda t : t*nu.subs(rule)
xL = np.linspace(0,s11initial,2)
fig, ax = plt.subplots(6,2, figsize = (20,30))
plt.setp(ax[0,0], xlabel = "s11", ylabel = "s22")
ax[0,0].plot(ff[:,0],ff[:,1],'b-')
ax[0,0].plot(xL,lambda1(xL),'b-')
def integralsolve(eq, ff):
    eq = eq.subs(Svm,Svmsubs)
    eptp = eq.subs({s11t(t):s11, s22t(t):s22, s11t(t).diff(t):1, nu:0.3, h:hslope, E:200000})
    f = sp.lambdify([s11,s22], eptp)
    sol = f(ff[:,0],ff[:,1])
    f2 = InterpolatedUnivariateSpline(ff[:,0],sol, k=1)     
    return [f2.integral(s11initial,ff[i,0]) for i in range(20)]
epfunction = integralsolve(eq1, ff)
plt.setp(ax[4,0], xlabel = "s11", ylabel = "e11")
ax[4,0].plot(ff[:,0],epfunction,'b-')
ax[4,0].plot(xL,[0,0], 'b-')
ep11function = integralsolve(eq2, ff)
plt.setp(ax[0,1], xlabel = "s11", ylabel = "ep11")
ax[0,1].plot(ff[:,0],ep11function,'b-')
ax[0,1].plot(xL,[0,0], 'b-')
ep22function = integralsolve(eq3, ff)
plt.setp(ax[0,1], xlabel = "s11", ylabel = "ep22")
ax[1,0].plot(ff[:,0],ep22function,'b-')
ax[1,0].plot(xL,[0,0], 'b-')
ep33function = integralsolve(eq4, ff)
plt.setp(ax[0,1], xlabel = "s11", ylabel = "ep33")
ax[1,1].plot(ff[:,0],ep33function,'b-')
ax[1,1].plot(xL,[0,0], 'b-')
# Now eefunctions
ee11function = (s11/E - s22/E*nu).subs(rule)
ee11_1 = sp.lambdify([s11, s22], ee11function)
sol_ee11 = ee11_1(ff[:,0],ff[:,1])
lambda2 = lambdify(t,(t/E-0.3*t/E*nu).subs(rule))
plt.setp(ax[2,0], xlabel = "s11", ylabel = "ee11")
ax[2,0].plot(ff[:,0],sol_ee11,'b-')
ax[2,0].plot(xL,lambda2(xL),'b-')
ee22function = (-nu*s11/E + s22/E).subs(rule)
ee22_1 = sp.lambdify([s11, s22], ee22function)
sol_ee22 = ee22_1(ff[:,0],ff[:,1])
#lambda3 = lambdify(t,((-nu*t)/E + (0.3*t)/E).subs(rule))
plt.setp(ax[2,1], xlabel = "s11", ylabel = "ee22")
ax[2,1].plot(ff[:,0],sol_ee22,'b-')
ax[2,1].plot(xL,[0,0],'b-')
ee33function = (-nu*s11/E - nu*s22/E).subs(rule)
ee33_1 = sp.lambdify([s11, s22], ee33function)
sol_ee33 = ee33_1(ff[:,0],ff[:,1])
lambda4 = lambdify(t,(-nu*t/E - nu*0.3*t/E).subs(rule))
plt.setp(ax[3,0], xlabel = "s11", ylabel = "ee33")
ax[3,0].plot(ff[:,0],sol_ee33,'b-')
ax[3,0].plot(xL,lambda4(xL),'b-')
plt.setp(ax[3,1], xlabel = "s11", ylabel = "e11")
ax[3,1].plot(ff[:,0],sol_ee11 + ep11function,'b-')
ax[3,1].plot(xL,lambda2(xL),'b-')
plt.setp(ax[4,1], xlabel = "s11", ylabel = "e22")
ax[4,1].plot(ff[:,0],sol_ee22 + ep22function,'b-')
ax[4,1].plot(xL,[0,0],'b-')
ax[4,1].set_ylim([-1,1])
plt.setp(ax[5,1], xlabel = "s11", ylabel = "e33")
ax[5,1].plot(ff[:,0],sol_ee33 + ep33function,'b-')
ax[5,1].plot(xL,lambda4(xL),'b-')
for i in range(6):
    for j in range(2):
        ax[i,j].grid(True,which='both')
        ax[i,j].axvline(x=0,color='k',alpha=0.5)
        ax[i,j].axhline(y=0,color = 'k',alpha=0.5)
plt.show()

Condition 3: \varepsilon_{22} = \varepsilon_{33}=0

In this loading condition, the Cauchy stress and the deviatoric stress matrices have the following forms:

    \[\sigma = \begin{pmatrix} \sigma_{11} & 0 & 0 \\ 0 & \sigma_{22} & 0 \\ 0 & 0 & \sigma_{33} \\ \end{pmatrix} \hspace{5mm} S = \begin{pmatrix} \frac{2\sigma_{11} - \sigma_{22} - \sigma_{33}}{3} & 0 & 0 \\ 0 & \frac{2\sigma_{22} - \sigma_{11} - \sigma_{33}}{3} & 0 \\ 0 & 0 & \frac{2\sigma_{33} - \sigma_{11} - \sigma_{22}}{3}\\ \end{pmatrix}\]

The onset of initial yielding occurs when \sigma_{vM}=400 MPa. Before developing any plastic strains, the total strain components \varepsilon_{22}and \varepsilon_{33} are equal to the elastic strain component \varepsilon_{22}^e and \varepsilon_{33}^e, respectively. The conditions of loading ensure that \varepsilon_{22}=\varepsilon_{33}=0 and thus:

    \[\varepsilon_{22}^e = \frac{\sigma_{22}}{E} - \frac{\nu \sigma_{11}}{E} - \frac{\nu \sigma_{33}}{E} = 0 \hspace{5mm} \varepsilon_{33}^e = \frac{\sigma_{33}}{E} - \frac{\nu \sigma_{11}}{E} - \frac{\nu \sigma_{22}}{E} = 0\]

Which lead to the following:

    \[\sigma_{22} = \sigma_{33} = \frac{\nu \sigma_{11}}{1-\nu}\]

Therefore, at the onset of initial yielding:

    \[\sigma_{11} = \frac{(1-\nu)400}{1-2\nu} = 700 MPa \hspace{5mm} \sigma_{22} = \sigma_{33} = \frac{\nu 700}{1-\nu} = 300 MPa\]

Therefore, when \sigma_{11}\leq 700 MPa, the plastic strain components are zero while the elastic strain components are given by:

    \[\varepsilon_{11}^e = \frac{\sigma_{11}}{E(1-\nu)}(1-\nu - 2\nu^2) \hspace{5mm} \varepsilon_{22}^e = \varepsilon_{33}^e = 0\]

Once yielding occurs, the relationship between the stress components obtained above no longer follows. We utilize the three equations: the consistency condition, \varepsilon_{22}=0, and \varepsilon_{33}=0 to find relationships between the four variables: \dot{\sigma}_{22}, \dot{\sigma}_{33}, \dot{\overline{\varepsilon^p}}, and \dot{\sigma}_{11}. First, the consistency condition can be used to relate the rate of change of the equivalent plastic strain parameter to the rates of change of the stress components:

    \[\dot{\phi} = \dot{\sigma}_{vM} - \dot{\sigma}_Y = 0 \Rightarrow \frac{\partial \sigma_{vM}}{\partial \sigma_{ij}} \dot{\sigma}_{ij} - \frac{d\sigma_Y}{d\overline{\varepsilon^p}} \dot{\overline{\varepsilon^p}} = 0 \Rightarrow \frac{3S_{11}}{2\sigma_{vM}} \dot{\sigma}_{11} + \frac{3S_{22}}{2\sigma_{vM}} \dot{\sigma}_{22} + \frac{3S_{33}}{2\sigma_{vM}} \dot{\sigma}_{33} - h \dot{\overline{\varepsilon^p}} = 0\]

Rearranging:

    \[\dot{\overline{\varepsilon^p}} = \frac{(2\sigma_{11} - \sigma_{22} -\sigma_{33}) \dot{\sigma}_{11} + (2\sigma_{22} - \sigma_{11} - \sigma_{33}) \dot{\sigma}_{22} + (2\sigma_{33} -\sigma_{22} -\sigma_{11})\dot{\sigma}_{33}}{2h\sigma_{vM}}\]

The boundary conditions dictate that \dot{\varepsilon}_{22}=\dot{\varepsilon}_{33}=\dot{\varepsilon}_{22}^e+\dot{\varepsilon}_{22}^p=\dot{\varepsilon}_{33}^e+\dot{\varepsilon}_{33}^p=0, which lead to the following relationships:

    \[\frac{\dot{\sigma}_{22}}{E} - \frac{\nu \dot{\sigma}_{11}}{E} - \frac{\nu \dot{\sigma}_{33}}{E} + \dot{\overline{\varepsilon^p}} \frac{\partial \phi}{\partial \sigma_{22}} = 0 \hspace{5mm} \frac{\dot{\sigma}_{33}}{E} - \frac{\nu \dot{\sigma}_{11}}{E} - \frac{\nu \dot{\sigma}_{22}}{E} + \dot{\overline{\varepsilon^p}} \frac{\partial \phi}{\partial \sigma_{33}} = 0\]

Utilizing the above relationships, the explicit relationships of \dot{\sigma}_{22} and \dot{\sigma}_{33} as functions of \dot{\sigma}_{11} can be obtained with the following forms:

    \[\dot{\sigma_{22}} = f_1 (\sigma_{11}, \sigma_{22}, \sigma_{33}, \dot{\sigma_{11}}, E, \nu, h) \hspace{5mm} \dot{\sigma_{33}} = f_2 (\sigma_{11}, \sigma_{22}, \sigma_{33}, \dot{\sigma_{11}}, E, \nu, h)\]

An ordinary differential equation for \dot{\overline{\varepsilon^p}} can also be obtained by substituting the obtained expressions f_1 and f_2 into the expression for \dot{\overline{\varepsilon^p}} and thus a third expression f_3 is obtained:

    \[\dot{\overline{\varepsilon^p}} = f_3 (\sigma_{11}, \sigma_{22}, \sigma_{33}, \dot{\sigma_{11}}, E, \nu, h)\]

The plastic strain components can then be calculated by integrating the following expressions for every i and j:

    \[\dot{\varepsilon}_{ij}^p = \dot{\overline{\varepsilon^p}} \frac{\partial \phi}{\partial \sigma_{ij}} = \dot{\overline{\varepsilon^p}} \frac{3S_{ij}}{2\sigma_{vM}}\]

The three ordinary differential equations for \dot{\sigma}_{22}, \dot{\sigma}_{33}, and \dot{\overline{\varepsilon^p}} can be integrated numerically by assuming that \sigma_{11} is equal to an arbitrary time parameter t. The initial conditions for the stresses are the elastic limits, while the initial conditions for the plastic strains are zero. The integration can be performed using a spreadsheet (using backward Euler numerical integration) or using the NDSolve[] function in Mathematica with an arbitrary upper limit for \sigma_{11}, which in this example was taken as 1050MPa.

The following figures show the plots of \sigma_{22} and the various nonzero plastic strain components as a function of \sigma_{11}. The plot of \sigma_{22} versus \sigma_{11} reveals two different regions during the loading process. In the initial elastic region, the stresses are related using Poisson’s ratio. In the second region after the development of plastic strains, the ratio of \sigma_{22} to \sigma_{11} is higher. The plastic strain components reveal that as the specimen is pulled beyond its initial elasticity limit, positive \varepsilon_{11}^p strain component is developed while negative \varepsilon_{22}^p and \varepsilon_{33}^p are encountered. Notice that the boundary and loading conditions ensure that \varepsilon_{22}^p=\varepsilon_{33}^p. The signs of the plastic strain components indicate that if the stresses are removed, the specimen will be longer in the direction of the basis vector e_1 while shorter in the other two directions. The sum of the axial plastic strain components is always equal to zero.

The following figures show the plots of the various nonzero elastic strain components as a function of \sigma_{11}. Notice that the elastic strain components \varepsilon_{22}^e=\varepsilon_{33}^e are equal to zero when \sigma_{11} is less than the initial elastic limit. Beyond that, \varepsilon_{22}^e increases rapidly similar to the plastic strain component \varepsilon_{22}^p but with a negative sign to maintain a zero value for the total strain \varepsilon_{22}.

The following figures show the plots of the only nonzero total strain component \varepsilon_{11} as a function of \sigma_{11}:

It is important to note that this condition predicts that the stresses can reach very high values before the actual failure of the material. As shown in the plot of \overline{\varepsilon^p} versus \sigma_{11} given below, the value of the plastic strain parameter is only equal to 0.0014 when \sigma_{11} is equal to 1050MPa. This is because the von Mises yield criterion in this condition only predicts development of plastic strains when there is a difference between the normal stress components \sigma_{11}, \sigma_{22}, and \sigma_{33}. However, the plane strain condition \varepsilon_{33}=0 along with the condition \varepsilon_{22}=0 leads to a very small increase of \sigma_{vM} with the increase in \sigma_{11}.

View Mathematica Code

Clear[Svm,s11,s22,s33]
(*stress, hydrostatic stress, deviatoric stress and parameters*)
s={{s11,0,0},{0,s22,0},{0,0,s33}};
hydr=(s11+s22+s33)/3;
Sdev=FullSimplify[s-hydr*IdentityMatrix[3]];
Sdev//MatrixForm
epmax=25/100-550/200000;
hslope=N[150/epmax];
rule={Ee->200000,nu->0.3,h->hslope};
(*Relationship between stress components in the elastic regime*)
a=Solve[{0==s22/Ee-nu*s33/Ee-nu*s11/Ee,0==s33/Ee-nu*s11/Ee-nu*s22/Ee},{s22,s33}]
a=a/.rule;
Sigmavm=Sqrt[1/2*((s11-s22)^2+(s11-s33)^2+(s22-s33)^2)];
(*s11,s22 and s33 at the elastic limit*)
b=FullSimplify[Solve[{(Sigmavm/.a[[1]])==400},s11],Assumptions->{0.5>nu>0,s11>0}]
s11initial=s11/.b[[2]]
s22initial=(s22/.a[[1]])/.b[[2]]
s33initial=(s33/.a[[1]])/.b[[2]]
(*Elastic strains right before plastifying*)
ee11=FullSimplify[(s11/Ee-nu*s22/Ee-nu*s33/Ee)/.a[[1]]]/.rule
ee22=FullSimplify[(-nu*s11/Ee+s22/Ee-nu*s33/Ee)/.a[[1]]]/.rule
ee33=FullSimplify[(-nu*s11/Ee-nu*s22/Ee+s33/Ee)/.a[[1]]]/.rule
(*Manipulation of rates*)
Clear[ee22,ee11,ee33]
Sdev=Sdev/.{s11->s11[t],s22->s22[t],s33->s33[t]}
Sigmavm=Sigmavm/.{s11->s11[t],s22->s22[t],s33->s33[t]}
ee22'[t]=s22'[t]/Ee-nu*s11'[t]/Ee-nu*s33'[t]/Ee;
ee33'[t]=s33'[t]/Ee-nu*s11'[t]/Ee-nu*s33'[t]/Ee;
cep'[t]=(3/2*Sdev[[1,1]]*s11'[t]/h/Svm+3/2*Sdev[[2,2]]*s22'[t]/h/Svm+3/2*Sdev[[3,3]]*s33'[t]/h/Svm)
cep22'[t]=FullSimplify[cep'[t]*3/2*Sdev[[2,2]]/Svm]
cep33'[t]=FullSimplify[cep'[t]*3/2*Sdev[[3,3]]/Svm]
(*Solving for s22'[t] and s33'[t] in terms of s11'[t] using the zero total strain conditions*)
a=FullSimplify[Solve[{ee22'[t]+cep22'[t]==0,ee33'[t]+cep33'[t]==0},{s22'[t],s33'[t]}]]
(*Solving the ordinary differential equations for s11'[t],s22'[t] and s33'[t]*)
eq1=(((s22'[t]/.a[[1]])/.rule)/.Svm->Sigmavm)
eq2=(((s33'[t]/.a[[1]])/.rule)/.Svm->Sigmavm)
ff=NDSolve[{s22'[t]==eq1,s33'[t]==eq2,s11'[t]==1,s11[s11initial]==s11initial,s22[s11initial]==s22initial, s33[s11initial]==s33initial},{s11,s22,s33},{t,s11initial,1050}]
s22function=Piecewise[{{nu*t/(1-nu)/.rule,0<=t<=s11initial},{s22[t]/.ff[[1]],t>s11initial}}]
s33function=Piecewise[{{nu*t/(1-nu)/.rule,0<=t<=s11initial},{s33[t]/.ff[[1]],t>s11initial}}]
Plot[s22function,{t,0,1050},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","s22(MPa)"},PlotStyle->Black]
Plot[s33function,{t,0,1050},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","s33(MPa)"},PlotStyle->Black]
(*Solving the ordinary differential equations for ep'[t] and ep_ij'[t]*) cep'[t]=(cep'[t]/.Svm->Sigmavm)/.rule
eq1=ep'[t]==FullSimplify[cep'[t]]
eq2=ep11'[t]==FullSimplify[(cep'[t])*3/2*Sdev[[1,1]]/Sigmavm]
eq3=ep22'[t]==FullSimplify[(cep'[t])*3/2*Sdev[[2,2]]/Sigmavm]
eq4=ep33'[t]==FullSimplify[(cep'[t])*3/2*Sdev[[3,3]]/Sigmavm]
rule2={s11[t]->t,s22[t]->s22function,s33[t]->s33function,s22'[t]->D[s22function,t],s33'[t]->D[s33function,t],s11'[t]->1}
eq1=eq1/.rule2;eq2=eq2/.rule2;eq3=eq3/.rule2;eq4=eq4/.rule2;
a=NDSolve[{eq1,ep[s11initial]==0},ep,{t,s11initial,1050}]
b=NDSolve[{eq2,ep11[s11initial]==0},ep11,{t,s11initial,1050}]
c=NDSolve[{eq3,ep22[s11initial]==0},ep22,{t,s11initial,1050}]
d=NDSolve[{eq4,ep33[s11initial]==0},ep33,{t,s11initial,1050}]
(*Plastic strains as functions of t and their plots*)
ep11function=Piecewise[{{0,0<=t<=s11initial},{ep11[t]/.b[[1]],t>s11initial}}]
ep22function=Piecewise[{{0,0<=t<=s11initial},{ep22[t]/.c[[1]],t>s11initial}}]
ep33function=Piecewise[{{0,0<=t<=s11initial},{ep33[t]/.d[[1]],t>s11initial}}]
epfunction=Piecewise[{{0,0<=t<=s11initial},{ep[t]/.a[[1]],t>s11initial}}]
ee11function=(t/Ee-s22function/Ee*nu-s33function/Ee*nu)/.rule;
ee22function=(-nu*t/Ee+s22function/Ee-s33function/Ee*nu)/.rule;
ee33function=(-nu*t/Ee-nu*s22function/Ee+s33function/Ee)/.rule;
Plot[ep11function,{t,0,1050},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","ep11(MPa)"},PlotStyle->Black]
Plot[ep22function,{t,0,1050},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","ep22(MPa)"},PlotStyle->Black]
Plot[ep33function,{t,0,1050},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","ep33(MPa)"},PlotStyle->Black]
Plot[epfunction,{t,0,1050},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","ep(MPa)"},PlotStyle->Black]
Plot[ee11function,{t,0,1050},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","ee11(MPa)"},PlotStyle->Black]
Plot[ee22function,{t,0,1050},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","ee22(MPa)"},PlotStyle->Black]
Plot[ee33function,{t,0,1050},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","ee33(MPa)"},PlotStyle->Black]
Plot[ee11function+ep11function,{t,0,1050},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","e11(MPa)"},PlotStyle->Black]
Plot[Chop[ee22function+ep22function,10^-7],{t,0,1050},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","e22(MPa)"},PlotStyle->Black]
Plot[Chop[ee33function+ep33function],{t,0,1050},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"s11(MPa)","e33(MPa)"},PlotStyle->Black]

View Python Code

from sympy import *
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.interpolate import InterpolatedUnivariateSpline
sp.init_printing(use_latex = "mathjax")
s11,s22,s33,nu,E,h,Svm,t = symbols("\u03c3_{11} \u03c3_{22} \u03c3_{33} nu E h \u03c3_{vM} t")
s = Matrix([[s11,0,0],[0,s22,0],[0,0,s33]])
hydr = (s11 + s22 + s33)/3;
Sdev = simplify(s-hydr*eye(3))
display(Sdev)
epmax = 25/100-550/200000
hslope = 150/epmax
display(epmax)
display(hslope)
rule = {E:200000, nu:0.3, h:hslope}
display("Relationship between stress components in the elastic regime")
a = solve((s22/E-nu*s33/E-nu*s11/E,s33/E-nu*s11/E-nu*s22/E),s22,s33)
for key in a:
    a[key] = a[key].subs(rule)
display(a)
Sigmavm = sqrt(1/2*((s11-s22)**2+(s11-s33)**2+(s22-s33)**2))
Sigmavmsub = simplify(Sigmavm.subs(a))
display("s11, s22 and s33 at the elastic limit")
b = solve((Sigmavmsub-400),sqrt(s11**2), dict = true)
display(b[0])
display("Since s_11>0:")
b = {s11:700}
display(b)
s11initial = s11.subs(b)
s22initial = (s22.subs(a)).subs(b)
s33initial = (s33.subs(a)).subs(b)
display(s11initial,s22initial,s33initial)
display("Elastic strains right before plastifying")
ee11 = simplify((s11/E-nu*s22/E-nu*s33/E).subs(a)).subs(rule)
ee22 = simplify((-nu*s11/E+s22/E-nu*s33/E).subs(a)).subs(rule)
ee33 = simplify((-nu*s11/E-nu*s22/E+s33/E).subs(a)).subs(rule) 
display(ee11,ee22,ee33)
display("Manipulation of rates")
s11t,s22t,s33t = Function('\u03c3_{11}')(t),Function('\u03c3_{22}')(t),Function('\u03c3_{33}')(t)
Sdev = Sdev.subs({s11:s11t,s22:s22t,s33:s33t})
Sigmavm = Sigmavm.subs({s11:s11t,s22:s22t,s33:s33t})
display(Sdev,Sigmavm)
ee22p = s22t.diff(t)/E-nu*s11t.diff(t)/E-nu*s33t.diff(t)/E
ee33p = s33t.diff(t)/E-nu*s11t.diff(t)/E-nu*s33t.diff(t)/E
cepp = 3/2*Sdev[0,0]*s11t.diff(t)/h/Svm+3/2*Sdev[1,1]*s22t.diff(t)/h/Svm+3/2*Sdev[2,2]*s33t.diff(t)/h/Svm
cep22p = simplify(cepp*3/2*Sdev[1,1]/Svm)
cep33p = simplify(cepp*3/2*Sdev[2,2]/Svm)
display(cepp,cep22p,cep33p)
display("Solving for s22'[t] and s33'[t] in terms of s11'[t] using the zero total strain conditions")
a = solve((ee22p+cep22p,ee33p+cep33p),s22t.diff(t),s33t.diff(t))
display(a)
display("Solving the ordinary differential equations for s11'[t],s22'[t] and s33'[t]")
eq1 = ((s22t.diff(t).subs(a)).subs(rule)).subs(Svm,Sigmavm)
eq2 = ((s33t.diff(t).subs(a)).subs(rule)).subs(Svm,Sigmavm)
equ1 = eq1.subs({s11t:s11,s22t:s22,s33t:s22,s11t.diff(t):1})
display(eq1,eq2)
equ2 = eq2.subs({s11t:s11,s22t:s33,s33t:s33,s11t.diff(t):1})
def rhs1(y,t):
    f = sp.lambdify([s11,s22],equ1)
    f = f(y[0],y[1])
    return [1,f]
def rhs2(y,t):
    f = sp.lambdify([s11,s33],equ2)
    f = f(y[0],y[1])
    return [1,f]
s11initial = float(s11initial)
t1 = np.linspace(s11initial,1050,20)
yzero1 = [s11initial,s22initial]
yzero2 = [s11initial,s33initial]
ff1 = odeint(rhs1,yzero1,t1)
ff2 = odeint(rhs2,yzero2,t1)
display("Solving the ordinary differential equations for ep'[t] and ep_ij'[t]")
cepp = (cepp.subs(rule)).subs(Svm,Sigmavm)
display(cepp)
eq1 = simplify(cepp)
eq2 = simplify(cepp*3/2*Sdev[0,0]/Sigmavm)
eq3 = simplify(cepp*3/2*Sdev[1,1]/Sigmavm)
eq4 = simplify(cepp*3/2*Sdev[2,2]/Sigmavm)
display(eq1,eq2,eq3,eq4)
rule2 = {s11t:s11,s22t:s22,s33t:s33,s11t.diff(t):1}
def integralsolve(eq, ff1,ff2):
    eptp = simplify((((eq.subs(a)).subs(Svm,Sigmavm)).subs(rule)).subs(rule2))
    f = sp.lambdify([s11,s22,s33],eptp)
    sol = f(ff1[:,0],ff1[:,1],ff2[:,1])
    f2 = InterpolatedUnivariateSpline(ff1[:,0],sol,k=1)     
    return [f2.integral(s11initial,ff1[i,0]) for i in range(20)]
epfunction = integralsolve(eq1,ff1,ff2)
ep11function = integralsolve(eq2,ff1,ff2)
ep22function = integralsolve(eq3,ff1,ff2)
ep33function = integralsolve(eq4,ff1,ff2)
#eeFunctions
ee11function = sp.lambdify([s11,s22,s33],(s11/E-s22/E*nu-s33/E*nu).subs(rule))
sol_ee11 = ee11function(ff1[:,0],ff1[:,1],ff2[:,1])
ee22function = sp.lambdify([s11,s22,s33],(-nu*s11/E+s22/E-s33/E*nu).subs(rule))
sol_ee22 = ee22function(ff1[:,0],ff1[:,1],ff2[:,1])
ee33function = sp.lambdify([s11,s22,s33],(-nu*s11/E-nu*s22/E+s33/E).subs(rule))
sol_ee33 = ee33function(ff1[:,0],ff1[:,1],ff2[:,1])
lambda1 = lambdify(t,(t*nu/(1-nu)).subs(rule))
lambda2 = lambdify(t,(t/E-nu*t/(1-nu)/E*nu-nu*t/(1-nu)/E*nu).subs(rule))
xL = np.linspace(0,s11initial,2)
fig,ax = plt.subplots(6,2, figsize = (20,30))
plt.setp(ax[0,0], xlabel = "s11",ylabel="s22")
ax[0,0].plot(ff1[:,0],ff1[:,1],'b-')
ax[0,0].plot(xL,lambda1(xL),'b-')
plt.setp(ax[0,1], xlabel = "s11",ylabel="s33")
ax[0,1].plot(ff2[:,0],ff2[:,1],'b-')
ax[0,1].plot(xL,lambda1(xL),'b-')
plt.setp(ax[1,0], xlabel = "s11",ylabel="ep")
ax[1,0].plot(ff1[:,0],epfunction,'b-')
ax[1,0].plot(xL,[0,0], 'b-')
plt.setp(ax[1,1], xlabel = "s11",ylabel="ep11")
ax[1,1].plot(ff1[:,0],ep11function,'b-')
ax[1,1].plot(xL,[0,0], 'b-')
plt.setp(ax[2,0], xlabel = "s11",ylabel="ep22")
ax[2,0].plot(ff1[:,0],ep22function,'b-')
ax[2,0].plot(xL,[0,0], 'b-')
plt.setp(ax[2,1], xlabel = "s11",ylabel="ep33")
ax[2,1].plot(ff1[:,0],ep33function,'b-')
ax[2,1].plot(xL,[0,0], 'b-')
plt.setp(ax[3,0], xlabel = "s11",ylabel="ee11")
ax[3,0].plot(ff1[:,0],sol_ee11,'b-')
ax[3,0].plot(xL,lambda2(xL),'b-')
plt.setp(ax[3,1], xlabel = "s11",ylabel="ee22")
ax[3,1].plot(ff1[:,0],sol_ee22,'b-')
ax[3,1].plot(xL,[0,0],'b-')
plt.setp(ax[4,0], xlabel = "s11",ylabel="ee33")
ax[4,0].plot(ff1[:,0],sol_ee33,'b-')
ax[4,0].plot(xL,[0,0],'b-')
plt.setp(ax[4,1], xlabel = "s11",ylabel="e11")
ax[4,1].plot(ff1[:,0],sol_ee11+ep11function,'b-')
ax[4,1].plot(xL,lambda2(xL),'b-')
plt.setp(ax[5,0], xlabel = "s11",ylabel="e22")
ax[5,0].plot(ff1[:,0],sol_ee22+ep22function,'b-')
ax[5,0].plot(xL,[0,0],'b-')
ax[5,0].set_ylim([-1,1])
plt.setp(ax[5,1], xlabel = "s11",ylabel="e33")
ax[5,1].plot(ff1[:,0],sol_ee33+ep33function,'b-')
ax[5,1].plot(xL,[0,0],'b-')
ax[5,1].set_ylim([-1,1])
for i in range(6):
    for j in range(2):
        ax[i,j].grid(True,which='both')
        ax[i,j].axvline(x=0,color='k',alpha=0.5)
        ax[i,j].axhline(y=0,color = 'k',alpha=0.5)
plt.show()

Example 2: Nonlinear Kinematic Hardening

The “true” stress versus “true strain” of a metal that follows the nonlinear kinematic hardening model with a constant size of the yield surface is given in the following table.

True Strain (%)00.090.182.0454.9059.49518.23523.975
Stress(MPa)0180360395450510570590

Assume that E=200 GPa and \nu=0.3. Find the values of C and \gamma such that the Armstrong-Frederick evolution law for the backstress can be used to describe the metal’s behaviour. Then, find the different nonzero elastic and plastic strain components and the components of the backstress tensor associated with the following loading steps:

  • Step 1: A specimen of this material is incrementally loaded in shear such that \sigma_{12}=300 MPa.

  • Step 2: The stress \sigma_{12} is removed

  • Step 3: \sigma_{11} is increased to 500 MPa.

Solution

The true plastic strain versus true stress can be computed by subtracting the elastic strain component and thus, the relationship between \overline{\varepsilon^p} and \sigma_Y is according to the following table:

\overline{\varepsilon^p}0000.01850.04680.09240.17950.2368
Stress(MPa)0180360395450510570590

The material starts to yield when \sigma_{11}=360 MPa and since the material is assumed to have a constant size for the yield surface, then \sigma_{Y\alpha}=360 MPa. Integration of the Armstrong-Frederick material model yields the following relationship for uniaxial loading:

    \[\alpha_{11} = \frac{C}{\gamma} (1 - e^{-\gamma \overline{\varepsilon^p}})\]

Where \alpha_{11}=\sigma_{11}-\sigma_{Y\alpha}. The values of C and \gamma that would produce the best fit between the given data and the model can be found to be equal to 2352.3 and 8.97417 respectively using Mathematica.

View Mathematica Code

Clear[Cc,gamma]
ss={{0,0},{0.0185,395-360},{0.0468,450-360},{0.0924,510-360},{0.1795,570-360},{0.2368,590-360}};
a=FindFit[ss,Cc/gamma*(1-Exp[-gamma*ep]),{Cc,gamma},ep]

Cc=Cc/.a[[1]]

gamma=gamma/.a[[2]]

ss2=Table[{ss[[i,1]],Cc/gamma*(1-Exp[-gamma*ss[[i,1]]])},{i,1,6}];
ListPlot[{ss,ss2},Joined->True,AxesOrigin->{0,0}]

View Python Code

import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np
def func(x,a,b):
    return a/b*(1-np.exp(-b*x)) 
xdata = np.array([0,0.0185,0.0468,0.0924,0.1795,0.2368])
ydata = np.array([0,35,90,150,210,230])
popt,pcov = curve_fit(func,xdata,ydata)
display("a value =",popt[0],"b value =",popt[1])
xfundata = np.array([i/20*0.25 for i in range(20)])
yfundata = func(xfundata,*popt)
plt.plot(xdata,ydata,'o',label='Data')
plt.plot(xfundata,yfundata,'r-',label='Best Fit: a=%5.1f, b=%5.5f'%tuple(popt))
plt.xlabel('x')
plt.ylabel('y')
plt.legend()

Step 1

At the beginning of the analysis, the values of the components \alpha_{ij} are equal to zero. As the material is first loaded in shear, the first initial yield is obtained when \sigma_{vM}=360 MPa, therefore:

    \[\sqrt{3\sigma_{12}^2} = 360 MPa \Rightarrow \sigma_{12} = 207.85 MPa\]

At this instant, the only nonzero elastic strain components are \varepsilon_{12}^e=\varepsilon_{21}^e=\frac{\sigma_{12}}{2G}. Beyond the elastic limit, an increment \Delta\sigma_{12} is accompanied by plastic strains and a transition of the yield surface according to the evolution of the backstress. The consistency condition can be used to find the value of \Delta\overline{\varepsilon^p}:

    \[\Delta f = 0 = \sum_{i,j=1}^3 \frac{\partial f}{\partial \sigma_{ij}} \Delta \sigma_{ij} + \sum_{i,j=1}^3 \frac{\partial f}{\partial \alpha_{ij}} \Delta \alpha_{ij} - \frac{d\sigma_{Y \alpha}}{d\overline{\varepsilon^p}} \Delta \overline{\varepsilon^p} = 0\]

Since the size of the yield surface is given in this problem as a constant, then the value of \frac{d\sigma_{Y\alpha}}{d\overline{\varepsilon^p}} is equal to 0. The only nonzero stress increment components in the above equations are \Delta \sigma_{12} and \Delta \sigma_{21}. Utilizing the Armstrong-Frederick evolution law to substitute for the components \Delta \alpha_{ij} results in the following expression:

    \[\frac{\partial f}{\partial \sigma_{12}} \Delta \sigma_{12} + \frac{\partial f}{\partial \sigma_{21}}\Delta \sigma_{21} + \sum_{i,j=1}^3 \Big(\frac{\partial f}{\partial \alpha_{ij}} \Big(\frac{C(\sigma_{ij} - \alpha_{ij})}{\sigma_{Y \alpha}} - \gamma \alpha_{ij}\Big)\Big) \Delta \overline{\varepsilon^p} = 0\]

The only nonzero components of the stress and the backstress tensor components would be \sigma_{12}, \sigma_{21}, \alpha_{12}, and \alpha_{21}. Therefore:

    \[\frac{\partial f}{\partial \sigma_{12}} \Delta \sigma_{12} + \frac{\partial f}{\partial \sigma_{21}} \Delta \sigma_{21} + \Big(\frac{\partial f}{\partial \alpha_{21}}\Big(\frac{C(\sigma_{12} - \alpha_{12})}{\sigma_{Y \alpha}} - \gamma \alpha_{12}\Big) + \frac{\partial f}{\partial \alpha_{21}} \Big(\frac{C(\sigma_{21} - \alpha_{21})}{\sigma_{Y \alpha}} - \gamma \alpha_{21} \Big) \Big) \Delta \overline{\varepsilon^p} = 0\]

The above relationship can be integrated numerically to find the value of \overline{\varepsilon^p} as \sigma_{12} increases from the initial yield at \sigma_{12}=207.85 MPa up to \sigma_{12}=300 MPa. The following table shows the numerical integration process with an increment \Delta \sigma_{12}=20 MPa.

\sigma_{12}=S_{12}\alpha_{12}\frac{\partial f}{\partial \sigma_{12}} = -\frac{\partial f}{\partial \alpha_{12}}\Delta \sigma_{12}\Delta \overline{\varepsilon^p}\Delta \alpha_{12}
207.8500.866200.014720
227.85200.866200.017020
247.85400.866200.02020
267.85600.866200.024420
287.85800.86612.150.019012.15
30092.15

The value of the plastic strain parameter at the end of the loading process is the sum of the individual increments:

    \[\overline{\varepsilon^p} = 0.0951\]

Also, the value of the non-zero backstress tensor components \alpha_{12}=\alpha_{21} is equal to the sum of the individual increments:

    \[\alpha_{12} = \alpha_{21} = 92.15 MPa\]

The nonzero elastic strain components are:

    \[\varepsilon_{12}^e = \varepsilon_{21}^e = \frac{300}{2G} = 0.00195\]

The nonzero plastic strain components can be calculated by numerical integration of the following formula:

    \[\Delta \varepsilon_{12}^p = \Delta \varepsilon_{21}^p = \Delta \overline{\varepsilon^p} \frac{\partial f}{\partial \sigma_{12}}\]

The following table shows the values of the plastic strain components at the end of this load step:

\frac{\partial f}{\partial \sigma_{12}}\Delta \overline{\varepsilon^p}\Delta \varepsilon_{12}^p = \Delta \varepsilon_{21}^p\sigma_{12} (MPa)\varepsilon_{12}^p = \varepsilon_{21}^p
0.8660.01470.0128207.850
0.8660.01700.0147227.850.0128
0.8660.02000.0173247.850.0274
0.8660.02440.0211267.850.0448
0.8660.01900.0164287.850.0659
3000.0824
Step 2

The unloading in step 2 occurs in the elastic regime since as \sigma_{12}is reduced to zero, the yield function is negative:

    \[\sigma_{vM}|_{\sigma - \alpha} = \sqrt{3(\sigma_{12} - 92.15)(\sigma_{12}-92.15)} = \sqrt{3} (\sigma_{12} - 92.15) < \sigma_{Y\alpha} = 360\]

At the end of step 2, the elastic strains are zero while no change occurs in the plastic strain components, the plastic strain parameter nor in the components of the backstress tensor:

    \[\varepsilon_{12}^p = \varepsilon_{21}^p = 0.0824 \hspace{5mm} \alpha_{12} = \alpha_{21} = 92.15 MPa \hspace{5mm} \overline{\varepsilon^p} = 0.0951\]

Step 3

Loading in step 3 will occur elastically until the yield function reaches the value of zero. In this instance, the value of \sigma_{11} can be obtained by solving the following equation:

    \[\sigma_{vM}|_{\sigma - \alpha} = \sqrt{3(-92.15)^2 + \frac{3}{2} \Big(\frac{2}{3} \sigma_{11} \Big)^2 + 3\Big(\frac{-1}{3} \sigma_{11} \Big)^2} = 360 \Rightarrow \sigma_{11} = 322.684 MPa\]

As \sigma_{11} increases beyond this limit, plastic strain develops and the consistency condition can be used to find the increment in the plastic strain parameter as follows:

    \[\Delta f = 0 = \sum_{i,j=1}^3 \frac{\partial f}{\partial \sigma_{ij}} \Delta \sigma_{ij} + \sum_{i,j=1}^e \frac{\partial f}{\partial \alpha_{ij}} \Delta \alpha_{ij} - \frac{d\sigma_{Y \alpha}}{d \overline{\varepsilon^p}} \Delta \overline{\varepsilon^p} = 0\]

In the above equation, \frac{d\sigma_{Y\alpha}}{d\overline{\varepsilon^p}}=0and the only nonzero stress increment component is \Delta \sigma_{11}. Utilizing the Armstrong-Frederick evolution law to substitute for the components \Delta\alpha_{ij} results in the following expression:

    \[\frac{\partial f}{\partial \sigma_{11}} \Delta \sigma_{11} + \sum_{i,j=1}^3 \Big(\frac{\partial f}{\partial \alpha_{ij}} \Big(\frac{C(\sigma_{ij} - \alpha_{ij})}{\sigma_{Y \alpha}} - \gamma \alpha_{ij} \Big) \Big) \Delta \overline{\varepsilon^p} = 0\]

The only nonzero components of the stress and the backstress tensor components would be \sigma_{11}, \alpha_{11}, \alpha_{12}, and \alpha_{21}. Therefore:

    \[\frac{\partial f}{\partial \sigma_{11}} \Delta \sigma_{11} + \Big(\frac{\partial f}{\partial \alpha_{11}} \Big(\frac{C(\sigma_{11} - \alpha_{11})}{\sigma_{Y \alpha}} -\gamma \alpha_{11}\Big) + \frac{\partial f}{\partial \alpha_{12}} \Big(\frac{C(\sigma_{12} -\alpha_{12})}{\sigma_{Y \alpha}} - \gamma \alpha_{12} \Big)\]

    \[+ \frac{\partial f}{\partial \alpha_{21}} \Big(\frac{C(\sigma_{21}-\alpha_{21})}{\sigma_{Y \alpha}} - \gamma \alpha_{21} \Big)\Big) \Delta \overline{\varepsilon^p} = 0\]

The above relationship can be integrated numerically to find the value of \overline{\varepsilon^p} as \sigma_{11} increases from the initial yield at \sigma_{11}=322.684 MPa up to \sigma_{11}=500 MPa. The following table shows the numerical integration process with an increment \Delta \sigma_{11}=20 MPa.

\sigma_{11}MPa\alpha_{11}\alpha_{12}=\alpha_{21} (MPa)\Delta \sigma_{11} (MPa)\frac{\partial f}{\partial \sigma_{11}} = -\frac{\partial f}{\partial \alpha_{11}}\frac{\partial f}{\partial \alpha_{12}}\Delta \overline{\varepsilon^p}\Delta \alpha_{11} (MPa)\Delta \alpha_{21} (MPa)
322.68092.15200.89630.38400.006012.65-8.58
342.6812.6583.58200.91680.34820.006713.50-8.56
362.6826.1575.02200.93480.32160.007314.33-8.48
382.6840.4766.53200.950.60.27720.008115.14-8.34
402.6855.6158.19200.96410.24250.009015.92-8.12
422.6871.5350.07200.97540.20860.010116.65-7.82
442.6888.1842.25200.98470.17600.011417.34-7.45
462.68105.5234.80200.99210.14500.012917.96-6.99
482.68123.4727.8117.3160.99780.11590.012916.02-5.58
500.00139.4922.23

The increments in the plastic strain components can be computed as follows:

    \[\Delta \varepsilon_{12}^p = \Delta \varepsilon_{21}^p = \Delta \overline{\varepsilon^p} \frac{\partial f}{\partial \sigma_{12}} = -\Delta \overline{\varepsilon^p} \frac{\partial f}{\partial \alpha_{12}} \hspace{5mm} \Delta \varepsilon_{11}^p = \Delta \overline{\varepsilon^p} \frac{\partial f}{\partial \sigma_{11}}\]

The following table shows the calculations of the increments in the plastic strain components:

\Delta \overline{\varepsilon^p}(MPa)\Delta \varepsilon_{12}^p\Delta \varepsilon_{11}^p (MPa)\sigma_{11} (MPa)\overline{\varepsilon^p}\varepsilon_{11}^p\varepsilon_{22}^p = \varepsilon_{33}^p\varepsilon_{12}^p = \varepsilon_{21}^p
0.0060-0.00230.0054322.680.0951000.0824
 0.0067-0.00230.0061342.680.10110.0054-0.00270.0801
 0.0073-0.00230.0068362.680.10770.0114-0.00570.0778
 0.0081-0.00220.0077382.680.11500.0183-0.00910.0755
 0.0090-0.00220.0087402.680.12310.0259-0.01300.0733
 0.0101-0.00210.0098422.680.13210.0346-0.01730.0711
 0.0114-0.00200.0112442.680.14220.0444-0.02220.0690
 0.0129-0.00190.0128462.680.15350.0556-0.02780.0670
 0.0129-0.00150.0129482.680.16650.0685-0.03420.0651
 500.000.17940.0814-0.04070.0636

Thus, at the end of step 3, the following are the values of the non-zero backstress components, elastic strain components, and plastic strain components:

    \[\sigma_{11} = 139.49 \hspace{5mm} \alpha_{12} = \alpha_{21} = 22.23 \hspace{5mm} \varepsilon_{11}^e = \frac{500}{E} = 0.0025 \hspace{5mm} \varepsilon_{22}^e = \varepsilon_{33}^e = -0.00075\]

    \[\varepsilon_{11}^p = 0.0814 \hspace{5mm} \varepsilon_{22}^p = \varepsilon_{33}^p = -0.0407 \hspace{5mm} \varepsilon_{12}^e = \varepsilon_{21}^p = 0.0636\]

This example shows the behaviour of the kinematic hardening model under multiple load steps. The application of a shear stress in an initial load step developed plastic shear strains and a backstress components \varepsilon_{12}. When the shear stress was removed and a uniaxial stress \sigma_{11} was applied, the material started to yield at a stress value of \sigma_{11}<360 MPa. In addition, contrary to the isotropic hardening model, this model predicts a decrease in the plastic shear strain components with the increase in \sigma_{11}.

Problems

  1. The engineering strain-engineering stress curve for a steel material follows the numbers in the table shown. Assuming that E=210,000MPa and Poisson’s ratio \nu=0.3.

    Strain (%)00.0850.172.065.029.952027.07
    Stress(MPa)0178.4356.4382428.5463.8472.5456.2

    • Plot the true strain – true stress curve.

    • Plot the relationship between the stress and the true plastic strain.

    • Find the yield stress given the two criteria: 0.2% offset and the 0.5% strain.

    • Find a Ramberg-Osgood approximation n,k for the relationship between \overline{\varepsilon^p} and \sigma_Y by minimizing the square root of the error in the observations compared to the model. Hints:

      • You can use one of the Mathematica functions NMinimize, Fit or FindFit to find n,k such that the experimental values would have the best fit to the equation \sigma_Y=E\left(\frac{\overline{\varepsilon^p}}{k}\right)^\frac{1}{n}

      • The Ramberg-Osgood approximation assumes that the point (0,0) belongs to the curve of the yield stress vs. plastic strain parameter. When curve fitting, you shouldn’t use a point with a yield stress corresponding to zero plastic strain.

    • Use the incremental Von-Mises isotropic hardening theory to find the elastic strains and the plastic strain parameter \overline{\varepsilon^p} for the following three plane stress states. Assume that the stresses increase linearly together starting from zero (i.e., no unloading).

      • State 1: \sigma_{11}=-134 MPa, \sigma_{22}=241 MPa, \sigma_{12}=53 MPa

      • State 2: \sigma_{11}=60 MPa, \sigma_{22}=370 MPa, \sigma_{12}=100 MPa

      • State 3: \sigma_{11}=100 MPa, \sigma_{22}=400 MPa, \sigma_{12}=120 MPa

    • Describe three different history paths that lead to the stress state \sigma_1=430 MPa and \sigma_2=\sigma_3=0

  2. In a three dimensional state of stress, during plastic deformation the increment in the plastic strain energy that is dissipated (irrecoverable plastic strain) in the material per unit volume is equal to dW_p=\sum_{i,j=1}^3\sigma_{ij} d\varepsilon_{ij}^p. Show that when the associated flow rule of von-Mises plasticity is used to calculate the plastic strain components, then the increment in the plastic strain energy is equal to dW_p=\sigma_{vM} d\overline{\varepsilon^p}, i.e., the product of the equivalent stress and the increment in the equivalent plastic strain.

  3. Let \phi=\sigma_{vM}-\sigma_Y be the von Mises yield function, where, \sigma_{vM} is the von Mises stress and \sigma_Y is the uniaxial yield stress. If \sigma,S\in\mathbb{M}^3 are the Cauchy and the deviatoric stress tensors respectively, show the following equality:

    \frac{\partial \phi}{\partial \sigma_{ij}} = \frac{3S_{ij}}{2\sigma_{vM}}

  4. If the stress-strain relationship of a material in simple tension is given by:

    \varepsilon = \varepsilon^e + \varepsilon^p = \frac{\sigma}{E} + \frac{\sigma^{m+1}}{E_t}

    Where m, E_t, and E are material constants.Using the incremental Von-Mises theory, show the following expression for the case of pure shear:

    \frac{d \tau}{d \gamma} = \frac{1}{\frac{1}{G} + \frac{3(m+1)}{E_t}(\sqrt{3} \tau)^m}

    Where \tau=\sigma_{12} is the applied shear stress, and \gamma=2\varepsilon_{12}^p+2\varepsilon_{12}^e is the total shear strain. This example is from the book “Plasticity for Structural Engineers” by Wai-Fah Chen.

  5. A pipe with capped ends (condition 1 in the figure below) and another infinitely long pipe (condition 2 in the figure below) are subjected to an internal pressure P. The two pipes are made out of a steel material with Young’s modulus E and Poisson’s ratio \nu. The true stress versus true strain relationship obtained from a uniaxial test follows a bilinear relationship with the initial slope equal to E and the strain hardening slope equal to E_p. The slope of the true stress versus plastic strain relationship is designated E_t. If the yield stress at the end of the initial elastic portion of the curve is equal to \sigma_Y then:

    • For EACH pipe, find the pressure P_Y (the pressure at the onset of yield) in terms of the material parameters given.

    • Find the relationship between E_t, E, and E_p.

    • Assuming that the pressure inside the pipes increases to P_Y+\Delta P, find the associated three dimensional plastic strain increments in EACH of the two cases.

  6. A long steel thin-walled tube with diameter D and wall thickness t is subjected to an interior pressure P_1 and an external pressure P_2. The ends of the tube are closed. The external pressure is assumed not to affect the axial stress component of the tube. Assuming that the material obeys the Von-Mises isotropic hardening theory and using the following relationship between the plastic strain parameter and the yield stress obtained from uniaxial experiments (a>0 is a material parameter):

    \overline{\varepsilon^p} = \alpha \sigma_{Y}^3

    Assuming that the axial direction is the direction of the basis vector e_1 and the circumferential direction is the direction of the basis vector e_2 show that the plastic strain components \varepsilon_{11}^p and \varepsilon_{22}^p developed through the following load paths are:

    (P_1, P_2) = (0,0) \to (P,1.5P) \hspace{5mm} (\varepsilon_{11}^p, \varepsilon_{22}^p) = (9aF/2, -9aF/2)

    (P_1,P_2) = (0,0) \to (0,1.5P) \to (P,1.5P) \hspace{5mm} (\varepsilon_{11}^p, \varepsilon_{22}^p) = (27aF/2, -27aF)

    (P_1,P_2) = (0,0) \to (P,0) \to (P,1.5P) \hspace{5mm} (\varepsilon_{11}^p, \varepsilon_{22}^p) = (0.9aF/2)

    Where F = \left(\frac{DP}{4t}\right)^3 Also, sketch the yield surfaces and the loading paths in the \sigma_{11},\sigma_{22} space and explain the results obtained. Hint: \sigma_{11}=\frac{P_1 D}{4t}, \sigma_{22} =\frac{(P_1-P_2 )D}{2t}, \sigma_{33}=0. This example is from the book “Plasticity for Structural Engineers” by Wai-Fah Chen

  7. In a confined compression test a compressive \sigma_{22} is applied while keeping \varepsilon_{11}, and \varepsilon_{33} equal to zero. If a metal specimen that follows the von Mises yield criterion and the associated flow rule is subjected to a confined compression test, as \sigma_{22} is increased, will the metal ever deform plastically or not. Justify your answer.

  8. A cylindrical metal piece is subjected to a lateral stress of -75MPa and a tensile vertical stress of 300MPa. Assume that the magnitudes of the stresses increase linearly together. If the yield stress is 200MPa and the relationship between the von Mises stress in MPa and the equivalent plastic strain \overline{\varepsilon^p} is linear with a slope of 10000:

    • Calculate the stresses at the onset of yielding.

    • Find the total strains at the onset of yielding.

    • Find the total strains when the stresses reach the maximum value.

  9. A piece of metal obeys the nonlinear kinematic hardening model with a constant size of the yield surface and the Armstrong-Frederick evolution law for the backstress. The following is the engineering stress-strain curve obtained when the material is stretched in a uniaxial stress state.

    Strain (%)00.0850.172.065.029.952027.07
    Stress(MPa)0178.4356.4382428.5463.8472.5456.2

    If a specimen of this material is subjected to the following load steps:

    • Step 1, a uniaxial stress component \sigma_{11} such that \varepsilon_{11}^p=0.015.

    • Step 2, the stress is removed.

    • Step 3, a uniaxial stress component \sigma_{22}.

    Find the true stress \sigma_{22} versus the “apparent” plastic strain component \varepsilon_{22}^p relationship obtained in step 3 using the nonlinear kinematic hardening model of the incremental plasticity with a constant size of the yield surface.

  10. The true stress true strain relationship of a metal in simple tension is given by:

    \varepsilon = \begin{cases} \frac{\sigma}{E} & (\sigma < \sigma_0) \\ \frac{\sigma}{E} + \frac{\sigma - \sigma_0}{E_p} & (\sigma \ge \sigma_0) \\ \end{cases}

    which is an elastic-linear hardening plastic stress strain relationship with \sigma_0 being the initial yield stress. The material follows the von Mises isotropic hardening plasticity model. Find the total strains in the following two different paths:

    • \sigma_{11} increases from 0 to \sigma_0 and then remains constant. When \sigma_{11} reaches \sigma_0, \sigma_{12} increases from 0 to \frac{\sigma_0}{\sqrt{3}}

    • \sigma_{11} and \sigma_{12} increase linearly together from (0,0) to (\sigma_0,\frac{\sigma_0}{\sqrt{3}})

  11. Show that the stress-strain relationship of a metal that follows the Armstrong-Frederick plasticity material model with a constant size of the yield surface is given by:

    \varepsilon = \begin{cases} \frac{\sigma}{E} & (\sigma < \sigma_0) \\ \frac{\sigma}{E} - \frac{1}{\gamma} \ln{1-\frac{\gamma(\sigma - \sigma_0)}{C}} & (\sigma \ge \sigma_0) \\ \end{cases}

    Where \gamma and C are the Armstrong-Frederick material constants and \sigma_0 is the initial yield stress.

  12. A metal follows the combined nonlinear kinematic-isotropic hardening plasticity model such that the kinematic hardening model follows the Armstrong-Frederick evolution law with C=3000 and \gamma=10 while the isotropic hardening model follows the relationship

    \sigma_{Y \alpha} = 200 + K \overline{\varepsilon^p}^n

    in units of MPa with K=100 and n=0.3.

    • Show that the relationship between the applied uniaxial stress \sigma_{11} and \varepsilon_{11}^p is given by:

      \sigma_{11} = 200 + K(\varepsilon_{11}^p)^n + \frac{C}{\gamma} (1-e^{-\gamma \varepsilon_{11}^p})

      Note that the above relationship between \sigma_{11} and \varepsilon_{11}^p implies that the value of the initial yield stress \sigma_0=200 MPa.

    • Find the total strains and the components of the backstress tensor at the end of the following path: \sigma_{11} increases from 0 to 200MPa and then remains constant. When \sigma_{11} reaches 200MPa, \sigma_{12} increases from 0 to 150MPa.

Leave a Reply

Your email address will not be published.