Solid Mechanics: Plasticity

Alessandro Mastrofini, Wolfram Research

#WolframTechConf

ABSTRACT

This talk will showcase elasto-plastic materials in Solid Mechanics, highlighting advances in simulating material behavior such as permanent deformation and hardening. With examples from metal and soil plasticity, attendees will discover the latest FEM capabilities for structural analysis and gain practical insight to apply them effectively.

Agenda

Introduction to plasticity

How to solve an elasto-plastic model

Implementation details

Yield criteria

Linear hardening

Permanent deformation

Soil plasticity

Still in development

Disclaimer:

this code is still under development.

It will be available in the next release.

Plasticity Monograph

Material properties

In materials science plasticity is the ability of a solid material to undergo permanent deformation.

r=0.75;
δ=1.37;
Graphics[
{
Rotate[{
FaceForm[LightGray],EdgeForm[Black],
Table[Disk[{i,0},r],{i,4,11,2r}],
Table[Disk[{i+r,δ},r],{i,4,11,2r}],
Table[Disk[{i,2*δ},r],{i,1,13,2r}],
Table[Disk[{i+r,3*δ},r],{i,1,13,2r}],
Table[Disk[{i,4*δ},r],{i,1,15,2r}],
Table[Disk[{i+r,5*δ},r],{i,1,17,2r}],
Table[Disk[{i,6*δ},r],{i,1,17,2r}],
Table[Disk[{i+r,7*δ},r],{i,1,5,2r}],
{Gray,Thickness[0.02],Arrowheads[0.1],Arrow[{{6,11},{12,11}}],Arrow[{{12,-4},{6,-4}}]}
},π/10]
,
{Red,Thickness[0.01],Line[{{-r-0.3,7*δ}-{0.4,0},{15r,-2}-{0.4,0}}]},
{Blue,Thickness[0.01],Line[{{18r-0.5,9*δ+0.3}-{0.3,0},{13r+0.2,-3}-{0.3,0}}]}
,
Translate[Rotate[{Red,Arrowheads[0.03],Arrow[{{0,6.5δ},{0,6.5δ}+{0,2}}],Arrow[{{0,6.5δ},{0,6.5δ}+{2,0}}]},-(π/4-0.04)],{-0.3,-0.2}]
,
Blue,Text[Style["Slip plane",16],{13,13},Left,Automatic],
Black,Text[Style["Macroscopic load",16],{13,-4},Left,Automatic]
}
,ImageSize->Large

]

Theoretical background

Tresca’s memoir in 1864

ListLinePlot[{Most[Transpose[{strainData,stressData}]],Transpose[{strainData/2.2,1.1*stressData}][[1;;175]]},Frame->True,FrameLabel->{{"Stress [MPa]",None},{"Strain [%]","Ductile vs. Brittle Material Stress-Strain Curve"}},Filling->Axis,PlotLegends->{"Ductile","Brittle"},PlotRange->{{0,30},{0,1100}},ImageSize->Full]

Theoretical background

Continuum approach

ListLinePlot[{Most[Transpose[{strainData,stressData}]],Transpose[{strainData/2.2,1.1*stressData}][[1;;175]]},Frame->True,FrameLabel->{{"Stress [MPa]",None},{"Strain [%]","Ductile vs. Brittle Material Stress-Strain Curve"}},Filling->Axis,PlotLegends->{"Ductile","Brittle"},PlotRange->{{0,38},{0,1600}},
Epilog->{
Style[Rectangle[{-1,-1},{3,550}],EdgeForm[Black],FaceForm[Opacity[0.15]]],
Line[{{1.5,550},{17,1400}}],
Line[{{3,270},{26,700}}],
Style[Rectangle[{17,700},{37,1500}],EdgeForm[Black],FaceForm[Opacity[0]]],
Style[Rectangle[{16.8,700},{17,1500}],Gray,Opacity[0.3]],
Style[Rectangle[{16.8,690},{37,700}],Gray,Opacity[0.3]],
Style[{Arrowheads[0.02],Arrow[{{21,1260},{24,1400}}]},Dashed,Red],
Style[Text["Linear elasticity",{26,1410}],Red],
Inset[q,{18,750},{Left,Bottom},18]
},ImageSize->Full]

Continuum approach and macroscale plasticity

Linear elasticity

Show[{
pt1={1,10};
pt2={4,10};
pt3={4,14};
pt4={4,6};
pt5={2,20};
ListLinePlot[{{{0,0},pt1}},PlotStyle->Thick,AxesLabel->{"ϵ","σ"},Ticks->None],
ListLinePlot[{{pt1,pt5}},PlotStyle->Dashed,PlotLabels->"Linear elasticty"],
ListPlot[{pt1},PlotStyle->{Red,PointSize[0.025]}],
ListLinePlot[{{0,10},pt1},PlotStyle->{Dashed,Black,Thin},PlotLabels->"Elastic limit"],Graphics[Inset["\!\(\*SubscriptBox[\(σ\), SubscriptBox[\(y\), \(0\)]]\)",{-0.1`,10}]]
},PlotRange->{{-0.5`,3},All},ImageSize->Large]

Continuum approach and macroscale plasticity

Beyond elasticity

Show[{
pt1={1,10};
pt2={4,10};
pt3={4,14};
pt4={4,6};
pt5={2,20};ListLinePlot[{{pt1,pt2},{pt1,pt3},{pt1,pt4}},PlotLabels->{"Perfect Plasticity","Hardening","Softening"},PlotRange->{{0,All},{0,All}},PlotStyle->{Thick,Thick,Thick},AxesLabel->{"ϵ","σ"},Ticks->None],
ListLinePlot[{{{0,0},pt1}},PlotStyle->Thick],
ListLinePlot[{{pt1,pt5}},PlotStyle->Dashed,PlotLabels->"Linear elasticity"],
ListPlot[{pt1},PlotStyle->{Red,PointSize[0.025]}],
ListLinePlot[{{0,10},pt1},PlotStyle->{Dashed,Black,Thin}],Graphics[Inset["\!\(\*SubscriptBox[\(σ\), SubscriptBox[\(y\), \(0\)]]\)",{-0.1`,10}]]},PlotRange->{{-0.5`,All},All},ImageSize->Large]

Continuum approach and macroscale plasticity

Beyond elasticity

Show[{
pt1={1,10};
pt2={4,10};
pt3={4,14};
pt4={4,6};
pt5={2,20};ListLinePlot[{{pt1,pt2},{pt1,pt3},{pt1,pt4}},PlotLabels->{"Perfect Plasticity",Style["Hardening",Opacity[.2]],Style["Softening",Opacity[.2]]},PlotRange->{{0,All},{0,All}},PlotStyle->{{Red,Thick},{LightGray,Thick},{LightGray,Thick}},AxesLabel->{"ϵ","σ"},Ticks->None],
ListLinePlot[{{{0,0},pt1}},PlotStyle->Thick],
ListLinePlot[{{pt1,pt5}},PlotStyle->Dashed,PlotLabels->"Linear elasticity"],
ListPlot[{pt1},PlotStyle->{Red,PointSize[0.025]}],
ListLinePlot[{{0,10},pt1},PlotStyle->{Dashed,Black,Thin}],Graphics[Inset["\!\(\*SubscriptBox[\(σ\), SubscriptBox[\(y\), \(0\)]]\)",{-0.1`,10}]]},PlotRange->{{-0.5`,All},All},ImageSize->Large]

How to solve an elasto-plastic problem?

Graphics[{
Arrowheads[0.02],
FontSize->10,
Text["Solve equilibrium PDEs (σ,ϵ)",{0,0}],
Δ=0.5;
Arrow[{{0,-0.1},{0,-Δ}}],
Text["Evaluate yielding function F(σ)",{0,-0.6}],
Red,Dashed,
Arrow[{{0,-0.7},{0,-0.7-Δ}}],
Text["Over elastic limit",{0,-0.7-Δ-0.1}],
Arrow[{{0,-0.7-Δ-0.2},{0,-0.7-Δ-0.2-Δ}}],
Text["Solve plastic flow",{0,-0.7-Δ-0.3-Δ}],
Arrow[{{0,-0.7-Δ-0.4-Δ},{0,-0.7-Δ-0.4-3Δ/2},{-3,-0.7-Δ-0.4-3Δ/2},{-3,-0.7-Δ-0.4-3Δ/2},{-3,0},{-1.2,0}}],
{Gray,Text["ϵ=(Subscript[ϵ, el])+(Subscript[ϵ, p])",{-2.2,-0.1}]},
Blue,Thick,
Arrow[{{1.2,-Δ-0.1},{1.8,-Δ-0.1}}],
Text["Elastic limit OK",{2.5,-0.6}],
Arrow[{{2.5,-Δ},{2.5,2Δ},{0,2Δ},{0,Δ+0.1}}],
Text["Increase boundary load",{0,Δ}],
Arrow[{{0,Δ-0.1},{0,0.1}}],
Opacity[0.2],
Gray,Polygon[{{-1.2,((-Δ)/2)},{1.5,((-Δ)/2)},{1.5,-2Δ},{3.5,-2Δ},{3.5,1.5},{-1.2,1.5}}],
Opacity[1],Text["SolidMechanicsPDE",{1.2,1.35}]
}]

Application: steel beam

Plasticity Monograph

dl=0.04;
δ=0.15;
Graphics[{
{FaceForm[Brown],EdgeForm[Black],Rectangle[{0,0},{10,1}]},
Table[Line[{{-0.2,i},{0,0.15+ i}}],{i,-0.12,0.9,0.1}],
{Thickness[0.01],Line[{{0,0},{0,1}}]},
Red,{Thickness[0.01], Line[{{10,0},{10,1}}]},
Line[{{10,0},{10.5,0}}], Arrow[{{10.5,0},{10.5,-1.15}}]
},ImageSize->Large]

Application: steel beam

Geometry

mesh=ToElementMesh[Rectangle[{0,0},{10,1}],"MaxCellMeasure"->0.1];

mesh["Wireframe"]

Application: steel beam

Variables

mesh=ToElementMesh[Rectangle[{0,0},{10,1}],"MaxCellMeasure"->0.1];

vars={{u[x,y],v[x,y],0},{x,y,z}};

Application: steel beam

Material parameters

mesh=ToElementMesh[Rectangle[{0,0},{10,1}],"MaxCellMeasure"->0.1];

vars={{u[x,y],v[x,y],0},{x,y,z}};

pars=<|"YoungModulus"->90*(10^9),"PoissonRatio"->0.3,
	"Thickness"->0.01,"SolidMechanicsModelForm"->"ExtendedPlaneStrain"|>;

Application: steel beam

Solid Mechanics

eqn[vars_,pars_,imposedDispl_]:={SolidMechanicsPDEComponent[vars,pars]=={0,0},SolidFixedCondition[x==0,vars,pars],
SolidDisplacementCondition[x==10,vars,pars,<|"Displacement"->{None, imposedDispl}|>]}

Application: steel beam

Solid Mechanics

eqn[vars_,pars_,imposedDispl_]:={SolidMechanicsPDEComponent[vars,pars]=={0,0},SolidFixedCondition[x==0,vars,pars],
SolidDisplacementCondition[x==10,vars,pars,<|"Displacement"->{None, imposedDispl}|>]}

Style[Column[eqn[vars,pars,δBC]],Small]

Application: steel beam

Solution

δ=-1.15;

elasticDispl=NDSolveValue[eqn[vars,pars,δ],vars[[1]],{x,y}∈mesh];

Application: steel beam

Deformed mesh

δ=-1.15;

elasticDispl=NDSolveValue[eqn[vars,pars,δ],vars[[1]],{x,y}∈mesh];

VectorDisplacementPlot[elasticDispl[[1;;2]],{x,y}∈mesh,]

Application: steel beam

Stress and yielding point

strain=SolidMechanicsStrain[vars,pars,elasticDispl];
stress=SolidMechanicsStress[vars,pars,strain];
vm=Head@VonMisesStress[vars,pars,stress];

ContourPlot[vm[x,y],{x,y}∈mesh,]

Application: steel beam

Stress and yielding point

{Px,Py}={0.3,0.98};

Show[mesh["Wireframe"],Graphics[{Red,PointSize[0.05],Point[{Px,Py}]}],]

σVME=vm[Px,Py]

Application: steel beam

Stress and yielding point

σVME=vm[Px,Py];

σy0=650 *(10^6);

Show[Plot[{σy0,ϵ*pars["YoungModulus"]},{ϵ,0,ϵf=0.018},],
ListPlot[{{{σVME/pars["YoungModulus"],σVME}}},],ListPlot[{{{σy0/pars["YoungModulus"],σy0}}},]]

Application: steel beam

Plasticity

parsPlastic=pars;
parsPlastic["YieldStress"]=σy0;

Application: steel beam

Plasticity

parsPlastic=pars;
parsPlastic["YieldStress"]=σy0;

RegionPlot[{vm[x,y]>=parsPlastic["YieldStress"],vm[x,y]<parsPlastic["YieldStress"]},{x,y}∈mesh,]

Application: steel beam

Elasto-plastic solution

maxIteration=100;
elastoplasticDispl=elasticDispl; 
updatePars=PDEModels`InitializePlasticity[vars,parsPlastic,elastoplasticDispl];

AbsoluteTiming[Monitor[
i=1;
While[PDEModels`PlasticStepQ[vars,updatePars] &&i<maxIteration,
updatePars= PDEModels`PlasticStrainStep[vars,updatePars,elastoplasticDispl];
elastoplasticDispl=NDSolveValue[eqn[vars,updatePars,δ],vars[[1]],{x,y}∈mesh];
i++;],
i]]

Application: steel beam

Elasto-plastic solution

Legended[Show[
mesh["Wireframe"[]],
ElementMeshDeformation[mesh, elasticDispl[[1;;2]], "ScalingFactor" -> 1][],
  ElementMeshDeformation[mesh, elastoplasticDispl[[1;;2]], "ScalingFactor" -> 1][],]
,SwatchLegend[{LightGray, Blue, Red}, {"Undeformed", "Elastic", "Elastoplastic"}]]

Application: steel beam

Elasto-plastic solution

strainEP=SolidMechanicsStrain[vars,updatePars,elastoplasticDispl];
stressEP=SolidMechanicsStress[vars,updatePars,strainEP];
vmEP=Head[VonMisesStress[vars,updatePars,stressEP]];

RegionPlot[{vmEP[x,y]>=parsPlastic["YieldStress"],vmEP[x,y]<parsPlastic["YieldStress"]},{x,y}∈mesh,]

Application: steel beam

Elasto-plastic solution

{Px,Py}={0.3,0.98};

σVME=vm[Px,Py];
σVMEP=vmEP[Px,Py];

Show[Plot[{updatePars["YieldStress"],ϵ*updatePars["YoungModulus"]},{ϵ,0,ϵf=0.018},PlotStyle->{Dashed,Dashed},],
ListLinePlot[{{0,0},{σVMEP/updatePars["YoungModulus"],σVMEP},{ϵf,σVMEP}}],
ListPlot[{{{σVME/updatePars["YoungModulus"],σVME}},{{σVME/updatePars["YoungModulus"],σVMEP}}},],ListPlot[{{{updatePars["YieldStress"]/updatePars["YoungModulus"],updatePars["YieldStress"]}}},]]

Application: steel beam

Plastic strain

plasticStrain=PDEModels`PlasticStrain[vars,updatePars]

Application: steel beam

Plastic strain

plasticStrain=PDEModels`PlasticStrain[vars,updatePars];

equivalentPlasticStrain=PDEModels`EquivalentStrain[vars,updatePars,plasticStrain]

ContourPlot[equivalentPlasticStrain,{x,y}∈mesh,]

Implementation details

Implementation details

Solve elastic - problem

Evaluate Yielding state

Plastic flow equations

Iterative solution (come back to 1)

Modeling plasticity

Yield criteria

von Mises yielding: F(σ,(Subscript[σ, y0]))=(Sqrt[(1/2)[((((Subscript[σ, 1])-(Subscript[σ, 2])))^2)+((((Subscript[σ, 2])-(Subscript[σ, 3])))^2)+((((Subscript[σ, 3])-(Subscript[σ, 1])))^2)])-(Subscript[σ, y0])

Tresca yielding: F(σ,(Subscript[σ, y0]))=(1/2)max[|(Subscript[σ, 1])-(Subscript[σ, 2])|,|(Subscript[σ, 2])-(Subscript[σ, 3])|,|(Subscript[σ, 3])-(Subscript[σ, 1])| ]-(Subscript[σ, y0])

GraphicsRow[{
Show[
Graphics3D[{myRot[graph_]:=Rotate[Rotate[graph,0.7853981633974483`,{0,-1,0}],0.7853981633974483`,{0,0,1}];Arrowheads[0.02`],Arrow[{{0,0,0},{5,0,0}}],Text["\!\(\*SubscriptBox[\(σ\), \(3\)]\)",{5.2`,0,0}],Arrowheads[0.02`],Arrow[{{0,0,0},{0,5,0}}],Text["\!\(\*SubscriptBox[\(σ\), \(2\)]\)",{0,5.3`,0}],Arrowheads[0.02`],Arrow[{{0,0,0},{0,0,-6.5`}}],Text["\!\(\*SubscriptBox[\(σ\), \(1\)]\)",{0,0,-7}],{Opacity[0.2`],LightBlue,myRot[Cylinder[{{0,0,-5},{0,0,5}},1]]},{Opacity[0.1`,Blue],EdgeForm[Black]},
CapForm["Butt"],Opacity[1],EdgeForm[Black],Opacity[0.2`,Blue],myRot[Tube[{{0,0,-5},{0,0,5}},1]]
,{Opacity[0.1],Black,Cuboid[{0,-3,-3},{0.001,3,3}]}},Method->{"TubePoints"->6}]
,ViewVertical->{1,0,0},ViewPoint->{1,1,-3},Boxed->False,ImagePadding->{{-120,10},{-80,-20}},ImageSize->Medium],

Show[ptSide=1.4;
RegionPlot[{Sqrt[(1/2)((((x-y))^2)+(((0-x))^2)+(((y-0))^2))]<=(ptSide),(1/2)Max[{Abs[x-y],Abs[y-0],Abs[0-x]}]<=(0.5 ptSide)},{x,-2,2},{y,-2,2},Axes->False,Frame->False],
Graphics[{Text["(Subscript[σ, 1])",{0,2.6}],Text["(Subscript[σ, 2])",{2.6,0}],Arrow[{{0,-1.7},{0,2.5}}],Arrow[{{-2.5,0},{2.5,0}}],
Gray,
Arrow[BezierCurve[{{-1.3,-1.5},{-1,-2},{-0.5,-2}}]],
Text["von Mises",{0,-2}],
Arrow[BezierCurve[{{0.8,-0.6},{1.1,-0.7},{1.2,-1.2}}]],
Text["Tresca",{1.3,-1.4}],
Opacity[0.1],Black,Rectangle[{-3,-3},{3,3}]}],PlotRange->{All,All},ImageSize->Small]
},ImageSize->Full]

Modeling plasticity

Yield criteria

Flow rule

(Subscript[(Overscript[ϵ, \[Application]]), p])=λ((∂(Subscript[Q, p]))/(∂σ))

where (Subscript[Q, p]) is the plastic potential.

Modeling plasticity

Yield criteria

Flow rule

Hardening models

(Subscript[σ, ys])=(Subscript[σ, (Subscript[y, 0])])+H (Subscript[ϵ, p])

GraphicsRow[{side=2;Show[
Graphics3D[]
,ViewVertical->{1,0,0},ViewPoint->{1,1,-3},Boxed->False,ImagePadding->{{-120,10},{-80,-20}},ImageSize->Large],
Show[

Graphics[{
Arrowheads[0.03],
Arrow[{{0,0},{0,side}}],Arrow[{{0,0},{side,0}}]
}],
RegionPlot[Sqrt[(1/2)((((x-y))^2)+(((0-x))^2)+(((y-0))^2))]<=(0.5 ptSide),{x,-2,2},{y,-2,2},Axes->False],
RegionPlot[Sqrt[(1/2)((((x-y))^2)+(((Sqrt[0.25]-x))^2)+(((y-0.25))^2))]<=(0.5ptSide),{x,-2,2},{y,-2,2},Axes->False,BoundaryStyle->Brown,PlotStyle->{Opacity[0.2],Brown}],
Graphics[]
,PlotRange->{All,All}]
},ImageSize->Full]

Modeling plasticity

Yield criteria

Flow rule

Hardening models

Show[{
pt1={1,10};
pt2={4,10};
pt3={4,14};
pt4={4,6};
pt5={2,20};ListLinePlot[{{pt1,pt2},{pt1,pt3},{pt1,pt4}},PlotLabels->{"Perfect Plasticity","Hardening","Softening"},PlotRange->{{0,All},{0,All}},PlotStyle->{Thick,Thick,Thick},AxesLabel->{"ϵ","σ"},Ticks->None],
ListLinePlot[{{{0,0},pt1}},PlotStyle->Thick],
ListLinePlot[{{pt1,pt5}},PlotStyle->Dashed,PlotLabels->"Linear elasticity"],
ListPlot[{pt1},PlotStyle->{Red,PointSize[0.025]}],
ListLinePlot[{{0,10},pt1},PlotStyle->{Dashed,Black,Thin}],Graphics[Inset["\!\(\*SubscriptBox[\(σ\), SubscriptBox[\(y\), \(0\)]]\)",{-0.1`,10}]]},PlotRange->{{-0.5`,All},All},ImageSize->Large]

Modeling plasticity

Small strain

Flow rule in term of the yielding function:

(Subscript[(Overscript[ϵ, \[Application]]), p])=λ((∂(Subscript[Q, p]))/(∂σ))⟶ ({{{, ({{(Subscript[(Overscript[ϵ, \[Application]]), p])=λ((∂F)/(∂σ))
(Overscript[α, \[Application]])=λ((∂F)/(∂q)),  }})}})

Where q=-((∂W)/(∂α)) is the back stress and α the hardening parameters.

Modeling plasticity

Small strain

Discrete  time  derivative:

({{{, ({{(Subscript[(Overscript[ϵ, \[Application]]), p])=λ((∂F)/(∂σ))
(Overscript[α, \[Application]])=λ((∂F)/(∂q)),  }})}})

({{{, ({{(1/(Δt))((Subscript[ϵ, p])-(Subscript[ϵ, p,(n)]))=λ((∂F)/(∂σ)),  }, {(1/(Δt))(α-(Subscript[α, (n)]))=λ((∂F)/(∂q)),  }})}})

Modeling plasticity

Small strain

({{{, ({{(Subscript[(Overscript[ϵ, \[Application]]), p])=λ((∂F)/(∂σ))
(Overscript[α, \[Application]])=λ((∂F)/(∂q)),  }})}})

({{{, ({{(1/(Δt))((Subscript[ϵ, p])-(Subscript[ϵ, p,(n)]))=λ((∂F)/(∂σ)),  }, {(1/(Δt))(α-(Subscript[α, (n)]))=λ((∂F)/(∂q)),  }})}})

({{{, ({{(Subscript[ϵ, p])-(Subscript[ϵ, p,(n)])-(γ-(Subscript[γ, (n)]))((∂F)/(∂σ))=0
α-(Subscript[α, (n)])-(γ-(Subscript[γ, (n)]))((∂F)/(∂q))=0
F(σ,(Subscript[σ, yo]),q)=0,  }})}})

Modeling plasticity

Predictor-corrector algorithm

Graphics[{
Arrowheads[0.02],
FontSize->10,
Text["Solve equilibrium PDEs (σ,ϵ)",{0,0}],
Δ=0.5;
Arrow[{{0,-0.1},{0,-Δ}}],
Text["Evaluate yielding function F(σ)",{0,-0.6}],
Red,Dashed,
Arrow[{{0,-0.7},{0,-0.7-Δ}}],
Text["Over elastic limit",{0,-0.7-Δ-0.1}],
Arrow[{{0,-0.7-Δ-0.2},{0,-0.7-Δ-0.2-Δ}}],
Text["Solve plastic flow",{0,-0.7-Δ-0.3-Δ}],
Arrow[{{0,-0.7-Δ-0.4-Δ},{0,-0.7-Δ-0.4-3Δ/2},{-3,-0.7-Δ-0.4-3Δ/2},{-3,-0.7-Δ-0.4-3Δ/2},{-3,0},{-1.2,0}}],
{Gray,Text["ϵ=(Subscript[ϵ, el])+(Subscript[ϵ, p])",{-2.2,-0.1}]},
Blue,Thick,
Arrow[{{1.2,-Δ-0.1},{1.8,-Δ-0.1}}],
Text["Elastic limit OK",{2.5,-0.6}],
Arrow[{{2.5,-Δ},{2.5,2Δ},{0,2Δ},{0,Δ+0.1}}],
Text["Increase boundary load",{0,Δ}],
Arrow[{{0,Δ-0.1},{0,0.1}}],
Opacity[0.2],
Gray,Polygon[{{-1.2,((-Δ)/2)},{1.5,((-Δ)/2)},{1.5,-2Δ},{3.5,-2Δ},{3.5,1.5},{-1.2,1.5}}],
Opacity[1],Text["SolidMechanicsPDE",{1.2,1.35}]
},ImageSize->800]

Modeling plasticity

Predictor-corrector algorithm

While[PDEModels`PlasticStepQ[vars,updatePars] &&i<maxIteration,

updatePars= PDEModels`PlasticStrainStep[vars,updatePars,elastoplasticDispl];
elastoplasticDispl=NDSolveValue[eqn[vars,updatePars,δ],vars[[1]],{x,y}∈mesh];
i++;
]

Yield criterion

Manipulate[Show[{pt1={s/10,s};
pt2={4,s};
pt3={4,s+5 s/20};
pt4={4,s-5 s/20};
pt5={25/10,25};
ListLinePlot[{{pt1,pt2},{pt1,pt3},{pt1,pt4}},PlotLabels->{"Perfect Plasticity",Style["Hardening",Opacity[.2]],Style["Softening",Opacity[.2]]},PlotRange->{{0,All},{0,All}},PlotStyle->{{Red,Thick},{LightGray,Thick},{LightGray,Thick}},AxesLabel->{"ϵ","σ"},Ticks->None],ListLinePlot[{{{0,0},pt1}},PlotStyle->Thick],ListLinePlot[{{pt1,pt5}},PlotStyle->Dashed,PlotLabels->"Linear elasticity"],ListPlot[{pt1},PlotStyle->{Red,PointSize[0.025]}],ListLinePlot[{{0,s},pt1},PlotStyle->{Dashed,Black,Thin}],Graphics[Inset["\!\(\*SubscriptBox[\(σ\), SubscriptBox[\(y\), \(0\)]]\)",{-0.1`,s}]]},PlotRange->{{-0.5`,All},All},ImageSize->Medium],{{s,12,"\!\(\*SubscriptBox[\(σ\), \(yield\)]\)"},5,20,0.1},SaveDefinitions->True]

Yield criterion

Stress tensor

(S=Table[Symbol["s" <> ToString@i <> ToString@j],{i,3},{j,3}])//MatrixForm

Yield criterion

von Mises

VonMisesStress[,<||>,Symmetrize[S]]

Yield criterion

von Mises

VonMisesStress[,<||>,Symmetrize[S]]>= (Subscript[σ, y0])

Tresca

(1/2)Max[Abs[Subtract@@#]&/@Subsets[{σ1,σ2,σ3},{2}]]>= (Subscript[τ, y0])

Yield criterion

parsPlastic=pars;
parsPlastic["ShearYieldStress"]=325 *(10^6);

tresca=(1/2)Max[Abs[Subtract@@#]&/@Subsets[{σ1,σ2,σ3},{2}]]

Yield criterion

parsPlastic=pars;
parsPlastic["ShearYieldStress"]=325 *(10^6);

tresca=(1/2)Max[Abs[Subtract@@#]&/@Subsets[{σ1,σ2,σ3},{2}]]

{principalStress,pricipalStressDirections}=PDEModels`PrincipalEigensystem[vars,pars,stress];

trescaStress = tresca/.Thread[Rule[{σ1,σ2,σ3},principalStress]];

Yield criterion

parsPlastic=pars;
parsPlastic["ShearYieldStress"]=325 *(10^6);

tresca=(1/2)Max[Abs[Subtract@@#]&/@Subsets[{σ1,σ2,σ3},{2}]]

{principalStress,pricipalStressDirections}=PDEModels`PrincipalEigensystem[vars,pars,stress];

trescaStress = tresca/.Thread[Rule[{σ1,σ2,σ3},principalStress]];

RegionPlot[{trescaStress>=parsPlastic["ShearYieldStress"],trescaStress <parsPlastic["ShearYieldStress"]},{x,y}∈mesh,]

Yield criterion

vm=Head@VonMisesStress[vars,pars,stress];

σY0=2 parsPlastic["ShearYieldStress"];

RegionPlot[{vm[x,y]>=σY0,
trescaStress>=parsPlastic["ShearYieldStress"]},{x,y}∈mesh,]

Hardening

Hardening

Linear hardening

Show[{
pt1={1,10};
pt2={4,10};
pt3={4,14};
pt4={4,6};
pt5={2,20};ListLinePlot[{{pt1,pt2},{pt1,pt3},{pt1,pt4}},PlotLabels->{Style["Perfect Plasticity",Opacity[.2]],Style["Hardening",Opacity[1]],Style["Softening",Opacity[.2]]},PlotRange->{{0,All},{0,All}},PlotStyle->{{Red,Opacity[0.2],Thick},{Red,Thick},{LightGray,Thick}},AxesLabel->{"ϵ","σ"},Ticks->None],
ListLinePlot[{{{0,0},pt1}},PlotStyle->Thick],
ListLinePlot[{{pt1,pt5}},PlotStyle->Dashed,PlotLabels->"Linear elasticity"],
ListPlot[{pt1},PlotStyle->{Red,PointSize[0.025]}],
ListLinePlot[{{0,10},pt1},PlotStyle->{Dashed,Black,Thin}],Graphics[Inset["\!\(\*SubscriptBox[\(σ\), SubscriptBox[\(y\), \(0\)]]\)",{-0.1`,10}]]},PlotRange->{{-0.5`,All},All},ImageSize->Large]

Hardening

ListAnimate[]

Hardening

pars["YieldStress"]=(Subscript[σ, (Subscript[y, 0])]);
pars["HardeningModulus"]=H;

Hardening

Holed plate

Show[
Region[Style[RegionDifference[Rectangle[{-2,-1},{2,1}],Disk[{0,0},0.33]],{FaceForm[LightGray],EdgeForm[Directive[Dotted,Black]]}]],
Region[Style[RegionDifference[Rectangle[{0,0},{2,1}],Disk[{0,0},0.33]],{FaceForm[Gray],EdgeForm[Black]}]]
,Graphics[{FaceForm[LightGray],EdgeForm[Black],
Red,Arrow[{{2,1.2},{2.6,1.2}}],{Thick,Line[{{2,0},{2,1.2}}]},
Text["δ",{2.25,1.35}],
Blue,Dashed,Line[{{-0.1,1},{-0.1,0}}],Line[{{0,-.1},{2,-.1}}]
}],ImageSize->Large]

Plasticity Monograph

Hardening

Holed plate

l=0.04;
h=0.02;
t=0.001;

reg=RegionDifference[Rectangle[{0,0},{l,h}],Disk[{0,0},h/3]]

Hardening

Holed plate

mesh=ToElementMesh[reg,"MaxCellMeasure"->1.5 (10^(-6))];

vars={{u[x,y],v[x,y],w[x,y]},{x,y,z}};

pars=<|"YoungModulus"->70*(10^9),"PoissonRatio"->0.3,"Thickness"->t,"SolidMechanicsModelForm"->"ExtendedPlaneStrain"|>;

Hardening

Holed plate

δ=1.5*(10^(-3));

Γsym={SolidDisplacementCondition[x==0,vars,pars,<|"Displacement"->{0,None}|>],SolidDisplacementCondition[y==0,vars,pars,<|"Displacement"->{None,0}|>]};

eqn[vars_,pars_,δ_]:={SolidMechanicsPDEComponent[vars,pars]=={0,0},Γsym,SolidDisplacementCondition[x==l,vars,pars,<|"Displacement"->{δ,None}|>]}

Hardening

Elastic solution

elasticDispl=NDSolveValue[eqn[vars,pars,δ],vars[[1]],{x,y}∈mesh];

VectorDisplacementPlot[elasticDispl[[1;;2]],]

strainE=SolidMechanicsStrain[vars,pars,elasticDispl];
stressE=SolidMechanicsStress[vars,pars,strainE];
vmE=Head@ VonMisesStress[vars,pars,stressE];

Hardening

Holed plate

parsPlastic=pars;
parsPlastic["YieldStress"]=250*(10^7);
parsPlastic["HardeningModulus"]=150*(10^8);

Hardening

Yield limit

RegionPlot[{vmE[x,y]>=parsPlastic["YieldStress"],vmE[x,y]<parsPlastic["YieldStress"]},{x,y}∈mesh,]

Hardening

Elasto-plastic solution

maxIteration=100;
elastoplasticDispl=elasticDispl; 
updatePars=PDEModels`InitializePlasticity[vars,parsPlastic,elastoplasticDispl];
Monitor[
i=1;
While[PDEModels`PlasticStepQ[vars,updatePars]&&i<maxIteration,
updatePars=PDEModels`PlasticStrainStep[vars,updatePars,elastoplasticDispl];
elastoplasticDispl=NDSolveValue[eqn[vars,updatePars,δ],vars[[1]],{x,y}∈mesh];
i++;],
i]

Hardening

Elasto-plastic solution

strainEP=SolidMechanicsStrain[vars,updatePars,elastoplasticDispl];
stressEP=SolidMechanicsStress[vars,updatePars,strainEP];
vmEP=Head@ VonMisesStress[vars,updatePars,stressEP];

Hardening

Elasto-plastic solution

pt={h/2.9,h/2.9};

Show[mesh[],Graphics[{Red,PointSize[.05],Point[pt]}]]

Hardening

Elasto-plastic solution

σVMEP=vmEP@@pt

σVME=vmE@@pt

Hardening

Elasto-plastic solution

σy0=parsPlastic["YieldStress"];

ϵ0=σy0/Y;

eCurve=ϵ*Y;

epCurve={{0,0},{ϵ0,σy0},{ϵf,σy0+H*(ϵf-ϵ0)}};

Hardening

Elasto-plastic solution

σy0=parsPlastic["YieldStress"];

ϵ0=σy0/Y;

eCurve=ϵ*Y;

epCurve={{0,0},{ϵ0,σy0},{ϵf,σy0+H*(ϵf-ϵ0)}};

rep={ϵf->0.065,σy0->parsPlastic["YieldStress"], Y-> parsPlastic["YoungModulus"],H->parsPlastic["HardeningModulus"]};

pointE = {σVME/Y,σVME}/.rep;
pointEP={σVME/Y,σVMEP}/.rep;

Hardening

Elasto-plastic solution

Show[
Plot[{σy0/.rep,eCurve/.rep},{ϵ,0,ϵf/.rep},],
ListLinePlot[epCurve/.rep,],
ListPlot[{{pointE},{pointEP}},]]/.rep

Load history

Load history

ListAnimate[]

Load history

Permanent deformation

dl=0.04;
δ=0.15;
Graphics[{
{FaceForm[Brown],EdgeForm[Black],Rectangle[{0,0},{10,1}]},
Table[Line[{{-0.2,i},{0,0.15+ i}}],{i,-0.12,0.9,0.1}],
{Thickness[0.01],Line[{{0,0},{0,1}}]},
Red,{Thickness[0.01], Line[{{10,0},{10,1}}]},
Line[{{10,0},{10.5,0}}], Arrow[{{10.5,0},{10.5,-1.15}}]
},ImageSize->Large]

mesh=ToElementMesh[Rectangle[{0,0},{10,1}],"MeshElementType"->TriangleElement,"MaxCellMeasure"->0.06];

Load history

Permanent deformation

mesh=ToElementMesh[Rectangle[{0,0},{10,1}],"MeshElementType"->TriangleElement,"MaxCellMeasure"->0.06];

vars={{u[x,y],v[x,y],0},{x,y,z}};

pars=<|"YoungModulus"->90*(10^9),"PoissonRatio"->0.3,"Thickness"->0.01,"SolidMechanicsModelForm"->"ExtendedPlaneStrain"|>;

parsPlastic=pars;
parsPlastic["YieldStress"]=650 *(10^6);

Load history

Permanent deformation

eqn[vars_,pars_,imposedDispl_]:={SolidMechanicsPDEComponent[vars,pars]=={0,0},SolidFixedCondition[x==0,vars,pars],
SolidDisplacementCondition[x==10,vars,pars,<|"Displacement"->{None,imposedDispl}|>]}

Load history

Permanent deformation

eqn[vars_,pars_,imposedDispl_]:={SolidMechanicsPDEComponent[vars,pars]=={0,0},SolidFixedCondition[x==0,vars,pars],
SolidDisplacementCondition[x==10,vars,pars,<|"Displacement"->{None,imposedDispl}|>]}

δ=-1.15;

elasticDispl=NDSolveValue[eqn[vars,pars,δ],vars[[1]],{x,y}∈mesh];

VectorDisplacementPlot[elasticDispl[[1;;2]],{x,y}∈mesh,]

Load history

Permanent deformation

strainE=SolidMechanicsStrain[vars,pars,elasticDispl];
stressE=SolidMechanicsStress[vars,pars,strainE];
vmE=Head@VonMisesStress[vars,pars,stressE];

RegionPlot[{vmE[x,y]>=parsPlastic["YieldStress"],vmE[x,y]<parsPlastic["YieldStress"]},{x,y}∈mesh,]

Load history

Permanent deformation

Solve elasto-plastic problem

maxIteration=100;
elastoplasticDispl=elasticDispl; 
updatePars=PDEModels`InitializePlasticity[vars,parsPlastic,elastoplasticDispl];
Monitor[
i=1;
While[PDEModels`PlasticStepQ[vars,updatePars] &&i<maxIteration,
updatePars=PDEModels`PlasticStrainStep[vars,updatePars,elastoplasticDispl];
elastoplasticDispl=NDSolveValue[eqn[vars,updatePars,δ],vars[[1]],{x,y}∈mesh];
i++;],
i]

Load history

Permanent deformation

Zero-load deformation

noLoadDispl=NDSolveValue[eqn[vars,updatePars,None],vars[[1]],{x,y}∈mesh];

Load history

Permanent deformation

Zero-load deformation

scalingFactor=2;

elasticDeformedMesh=ElementMeshDeformation[mesh, elasticDispl[[1;;2]], "ScalingFactor" ->scalingFactor];

elastoplasticDeformedMesh=ElementMeshDeformation[mesh, elastoplasticDispl[[1;;2]], "ScalingFactor" ->scalingFactor];

plasticizedMesh=ElementMeshDeformation[mesh,noLoadDispl[[1;;2]],"ScalingFactor"->scalingFactor];

Load history

Permanent deformation

Zero-load deformation

Legended[Show[
mesh[],elasticDeformedMesh[],elastoplasticDeformedMesh[],plasticizedMesh[],],]

Load history

Permanent deformation

Stress and strain

pt={0.3,0.98};

Show[mesh[],Graphics[{Red,PointSize[.03],Point[pt]}]]

Load history

Permanent deformation

Stress and strain

strainEP=SolidMechanicsStrain[vars,updatePars,elastoplasticDispl];
stressEP=SolidMechanicsStress[vars,updatePars,strainEP];
vmEP=Head@ VonMisesStress[vars,updatePars,stressEP];

σE=vmE@@pt

σEP=vmEP@@pt

Load history

Permanent deformation

Stress-strain limit

ϵtotal=σE/pars["YoungModulus"];

ϵyield=parsPlastic["YieldStress"]/pars["YoungModulus"];

p1={ϵyield,parsPlastic["YieldStress"]};

p2={ϵtotal,σEP};

line=y-p2[[2]]==pars["YoungModulus"] (x-p2[[1]]);
unloading=SolveValues[line,y][[1]]

permanentStrain=SolveValues[unloading==0,x][[1]]

p3={permanentStrain,0};

Load history

Permanent deformation

Stress-strain curve

Show[Plot[{σEP,ϵ*pars["YoungModulus"]},{ϵ,0,ϵf=0.016},],
ListLinePlot[{{0,0},{ϵyield,σEP},{ϵf,σEP}}],
ListPlot[{{p1},{p2},{p3}},PlotStyle->{Red,Green,Blue},PlotMarkers->{"●",Medium},PlotLabels->Placed[{"P1","P2","P3"},Above]],
Graphics[Arrow[{p2,{permanentStrain,0}}]]
]

Soil plasticity

p=Cuboid[{0,0,0},{10,5,3}]; pg=Cuboid[{0,0,0},{-10,5,3}];
l=Cuboid[{0,0,3},{1,5,3.2}]; lg=Cuboid[{0,0,3},{-1,5,3.2}];
Graphics3D[{EdgeForm[],FaceForm[],p,pg, 
EdgeForm[], FaceForm[Gray],l,lg,
Opacity[0.1],Cuboid[{-0.01,-2,0},{0.01,5+2,3+2}]
},Boxed->False]

Plasticity Monograph

Soil plasticity

p2D=Rectangle[{0,0},{5,3}];
l2D=Rectangle[{0,3},{1,3.2}];
p2Dg=Rectangle[{0,0},{-5,3}];
l2Dg=Rectangle[{0,3},{-1,3.2}];

Graphics[ {FaceForm[Brown], p2D, FaceForm[Gray], l2D, {Opacity[0.2],
FaceForm[Gray], p2Dg, l2Dg },
Black, Table[Arrow[{{x,4},{x,3.25}}],{x,0.1,1,0.2}],
Gray,Line[{{-0.2,0},{-0.2,3}}],Table[Circle[{-0.1,y},0.1],{y,0.3,2.8,0.35}],Line[{{5.2,0},{5.2,3}}],Table[Circle[{5+0.1,y},0.1],{y,0.3,2.8,0.35}],Line[{{0,0},{5,0}}],Table[Line[{{x,-0.1},{x+0.1,0}}],{x,0.,5.8,0.15}],
{Dashed,Line[{{0,-0.4},{0,4.5}}]}, Gray, Arrowheads[0.025],
Arrow[{{2.5,-.5},{0.,-0.5}}],Arrow[{{2.5,-.5},{5.,-0.5}}],
Arrow[{{.6,4.5},{0.,4.5}}],Arrow[{{.6,4.5},{1.2,4.5}}],
Arrow[{{-.6,1.5},{-.6,0}}],Arrow[{{-.6,1.5},{-.6,3}}],
Text["0.6 m", {.6, 4.7}], Text["5 m", {2.4,-.7}], Text["3 m", {-1,1.5}],
{White,Text["soil", {2.5,1.5}]},
{Black,Text["footing", {1.6,3.4}]}
}]

Soil plasticity

h=3;
l=5;
footing=0.6;

mesh=ElementMeshRegionProduct[
ToGradedMesh[Line[{{0}, {l}}],],
 ToGradedMesh[Line[{{0}, {h}}],]];

mesh["Wireframe"]

Soil plasticity

vars={{u[x,y],v[x,y],0},{x,y,z}};

pars=<|"YoungModulus"->50*(10^6),"PoissonRatio"->0.45,"Thickness"->l,"SolidMechanicsModelForm"->"ExtendedPlaneStrain"|>;

eqn[vars_,pars_]:={SolidMechanicsPDEComponent[vars,pars]==SolidBoundaryLoadValue[x<=footing && y==h,vars,pars,<|"Force"->{None,-load}|>],
SolidDisplacementCondition[x==0||x==l,vars,pars,<|"Displacement"->{0,None}|>],
SolidDisplacementCondition[y==0,vars,pars,<|"Displacement"->{None,0}|>]
}

Soil plasticity

load=1.52*(10^6);

elasticDispl=NDSolveValue[eqn[vars,pars],vars[[1]],{x,y}∈mesh];

Show[mesh[],ElementMeshDeformation[mesh, elasticDispl[[1;;2]],"ScalingFactor"->20][]]

Soil plasticity

strain=SolidMechanicsStrain[vars,pars,elasticDispl];
stress=SolidMechanicsStress[vars,pars,strain];

Soil plasticity

Mohr-Coulomb: F(σ,c, cosϕ)=(Sqrt[(Subscript[J, 2])(σ)])m(θ)  + (Subscript[α, MC]) (Subscript[I, 1])(σ)-(Subscript[k, MC])

Drucker-Prager: F(σ,(Subscript[α, DP]),(Subscript[k, DP]))=(Sqrt[(Subscript[J, 2])])+((Subscript[α, DP])/3) (Subscript[I, 1])-(Subscript[k, DP])

Soil plasticity

ϕMC=30 Degree;
cMC=30 (10^3);

αMC=Sin[ϕMC]/3;
kMC=cMC Cos[ϕMC];

Soil plasticity

αDP=(2/(Sqrt[3]))((Sin[ϕMC])/(3-Sin[ϕMC]));
kDP=(6/(Sqrt[3]))((cMC Cos[ϕMC])/(3-Sin[ϕMC]));

Soil plasticity

I1=Tr[stress];

vm=VonMisesStress[vars,pars,stress];
J2Sqrt=vm/Sqrt[3];

devS=stress-(1/3)Tr[stress] IdentityMatrix[3];

arg=Clip[(((((Det[devS])/(J2Sqrt^2))))^3),{-1,1}];
theta=(1/3)ArcCos[arg];
m=(1/(Sqrt[3]))(1+Sin[ϕMC])Cos[theta] -(1-Sin[ϕMC])Cos[theta+2Pi/3];

Soil plasticity

Mohr-Coulomb: F(σ,c, cosϕ)=(Sqrt[(Subscript[J, 2])(σ)])m(θ)  + (Subscript[α, MC]) (Subscript[I, 1])(σ)-(Subscript[k, MC])

mcCriterion=J2Sqrt m+αMC *I1-kMC;

Drucker-Prager: F(σ,(Subscript[α, DP]),(Subscript[k, DP]))=(Sqrt[(Subscript[J, 2])])+((Subscript[α, DP])/3) (Subscript[I, 1])-(Subscript[k, DP])

dpCriterion=J2Sqrt+αDP*I1/3-kDP;

RegionPlot[{dpCriterion>=0,mcCriterion>=0,dpCriterion <0&&mcCriterion <0},{x,y}∈mesh,]

Soil plasticity

parsPlasticMC=pars;
parsPlasticMC["Cohesion"]=cMC;
parsPlasticMC["AngleInternalFriction"]=ϕMC;

parsPlasticMC["Plasticity"]=<|"YieldCriterion"->"MohrCoulomb"|>;

maxIteration=150;
elastoplasticDisplMC=elasticDispl; 
updateParsMC=PDEModels`InitializePlasticity[vars,parsPlasticMC,elastoplasticDisplMC];
Monitor[
i=1;
While[PDEModels`PlasticStepQ[vars,updateParsMC] &&i<maxIteration,
updateParsMC= PDEModels`PlasticStrainStep[vars,updateParsMC,elastoplasticDisplMC];
elastoplasticDisplMC=NDSolveValue[eqn[vars,updateParsMC],vars[[1]],{x,y}∈mesh];
i++;],
i]

Soil plasticity

parsPlasticDP=pars;
parsPlasticDP["CohesionParameter"]=kDP;
parsPlasticDP["PressureSensitivity"]=αDP;

parsPlasticDP["Plasticity"]=<|"YieldCriterion"->"DruckerPrager"|>;

maxIteration=150;
elastoplasticDisplDP=elasticDispl; 
updateParsDP=PDEModels`InitializePlasticity[vars,parsPlasticDP,elastoplasticDisplDP];
Monitor[
i=1;
While[PDEModels`PlasticStepQ[vars,updateParsDP] &&i<maxIteration,
updateParsDP= PDEModels`PlasticStrainStep[vars,updateParsDP,elastoplasticDisplDP];
elastoplasticDisplDP=NDSolveValue[eqn[vars,updateParsDP],vars[[1]],{x,y}∈mesh];
i++;],
i]

Soil plasticity

Show[
mesh[],
ElementMeshDeformation[mesh, elasticDispl[[1;;2]], "ScalingFactor" -> 10][],
  ElementMeshDeformation[mesh, elastoplasticDisplMC[[1;;2]], "ScalingFactor" -> 10][],
 ElementMeshDeformation[mesh, elastoplasticDisplDP[[1;;2]], "ScalingFactor" -> 10][],]
SwatchLegend[{LightGray, Blue,DarkRed, Red}, {"Undeformed", "Elastic", "Elastoplastic with Mohr-Coulomb","Elastoplastic with Drucker-Prager"}]

Soil plasticity

strainEP=SolidMechanicsStrain[vars,updateParsMC,elastoplasticDisplMC];
stressEP=SolidMechanicsStress[vars,updateParsMC,strainEP];

{{σ1,σ2,σ3},directions}=PDEModels`PrincipalEigensystem[vars,updateParsMC,stressEP];

{e1,e2,e3}=directions[[All,1;;2]];

Soil plasticity

Legended[StreamPlot[e3,{x,y}∈mesh,],]

Soil plasticity

ω = π/4 -ϕMC/2;

np=e1 Sin[ω]+e3 Cos[ω];
nm=-e1 Sin[ω]+e3 Cos[ω];

τmax=ElementMeshInterpolation[mesh,Abs[((Head[σ3]["ValuesOnGrid"]-Head[σ1]["ValuesOnGrid"])/2)]];

Show[
DensityPlot[τmax[x,y],{x,y}∈mesh,],
StreamPlot[np,{x,y}∈mesh,],
StreamPlot[nm,{x,y}∈mesh,],
]

Conclusion

Conclusion

GraphicsGrid[{
{,,},
{,,}
}]

Conclusion

paclet:PDEModels/tutorial/StructuralMechanics/Plasticity

GraphicsGrid[{
{,},
{,}}]

Auxiliaries

Needs["NDSolve`FEM`"]

steelStressStrainData=Import[FileNameJoin[{$InstallationDirectory,"SystemFiles","Components","PDEModels","SupportFiles","SteelStressStrainData.xls"}]][[1]];
loadData=Quantity[steelStressStrainData[[All,2]]*10^3,"Newtons"];
strainData=steelStressStrainData[[All,3]]*100;
specimenArea=Quantity[(6.37/2)^2*π,"Millimeters"^2];
stressData=UnitConvert[loadData/specimenArea,"MegaPascals"];

q=ListLinePlot[{Most[Transpose[{strainData,stressData}]],Transpose[{strainData/2.2,1.1*stressData}][[1;;175]],
Table[{x,3*(10^3) x},{x,0,0.2,0.1}],
Table[{x,1.3*(10^3) x},{x,0,0.4,0.1}]
},Frame->False,AxesLabel->{"Strain [%]","Stress [MPa]"},
Filling->{1->2,2->2,None,None},PlotRange->{{0,3},All},PlotStyle->{{Dashed,Gray},{Dashed,Gray},Red,Red}];

DumpSave[FileNameJoin[{NotebookDirectory[],"session.mx"}],"Global`"]

Get[FileNameJoin[{NotebookDirectory[],"session.mx"}]]