Solid Mechanics: Plasticity

Alessandro Mastrofini, Wolfram Research

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

ABSTRACT

This talk will highlight FEM modeling for Plasticity in Solid Mechanics, focusing on key updates in material behavior simulation, including plastic deformation and hardening . Attendees will gain insight into current FEM capabilities for structural analysis and be equipped with the tools to start their own simulations .

Agenda

Introduction to plasticity

Computational plasticity

Implementation details

Cook’s membrane application example

Still in development

Disclaimer:

this code is still under development

Material properties

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

Theoretical background

Tresca’s memoir in 1864

Theoretical background

Continuum approach

Continuum approach

Modeling elasto-plastic materials

Yield criteria

Flow rules

Hardening models

Modeling elasto-plastic materials

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

Modeling elasto-plastic materials

Yield criteria

Flow rules

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

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

basePlot=Show[
Graphics[{
Arrow[{{0,0},{0,2}}],
Arrow[{{0,0},{2,0}}],
Text["(Subscript[σ, 1])",{2.1,0}],Text["(Subscript[σ, 2])",{0,2.1}],
{Brown,Text["(Subscript[σ, ij])",{0.6,0.5}]}
}],
Region[Style[ImplicitRegion[Sqrt[(1/2)((((x-y))^2)+(((0-x))^2)+(((y-0))^2))]==(1),{x,y}],Blue]]
];

GraphicsRow[{Show[basePlot,
Graphics[{
{Brown,Arrow[{{0,0},{1,1}}]},
{Red,Line[{{0.8,1.2},{1.2,0.8}}],
Text["((dF (σ))/(Subscript[σ, ij]))",{1.3,1}],
Thick,Arrow[{{1,1},{1.4,1.4}}]}
}],
Rotate[Region[Style[ImplicitRegion[Sqrt[(1/2)((((x-y))^2)+(((0-x))^2)+(((y-0))^2))]==(1),{x,y}],Red,Dashed]],-π/20,{1,1}]
],
SwatchLegend[{Blue,Red},{"Yield surface","Potential surface"}]
}]

Modeling elasto-plastic materials

Yield criteria

Flow rules

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

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

Show[basePlot,
Graphics[{
{Red,Line[{{0.9,1.15},{1.2,0.7}}],Thick,Arrow[{{1.6,1.2},{1.05,0.9}}]},{Brown,Arrow[{{0,0},{1.6,1.2}}]}
}]]

Modeling elasto-plastic materials

Yield criteria

Flow rules

Hardening models

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

Small strain plasticity

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.

Small strain plasticity

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)),  }})}})

Small strain plasticity

({{{, ({{(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,  }})}})

Predictor - corrector algorithm

Implementation

Solve elastic - problem

Evaluate Yielding state

Plastic flow equations

Iterative solution (come back to 1)

Implementation

1. Solve elastic problem

NDSolveValue[{SolidMechanicsPDEComponent[vars,pars]=={0,0},BCs},vars[[1]],vars[[-1]]∈mesh]

u, (Subscript[ϵ, el]), σ

Implementation

2. Evaluate Yielding state

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

(Subscript[σ, vm])=VonMisesStress[vars,pars,stress];

Implementation

2. Evaluate Yielding state

Manipulate[
(*σ0=5*(10^8);*)
ϵ0=σ0/Y;

Show[
Plot[Y e/100,{e,0,3},PlotRange->{{0,3},{0,6*(10^8)}},PlotStyle->{Gray,Dashed}],
ListLinePlot[{{0,0},{100*ϵ0,σ0},{50*100*ϵ0,σ0}}],
Graphics[{
{
If[Y l>σ0,Gray,Red]
,PointSize[Large],Point[{l*100,Y l}]},If[l>ϵ0,{Arrow[{{l*100,Y l},{l*100,σ0}}],{Red,PointSize[Large],Point[{l*100,σ0}]}}]}]
]
,
{{l,1.8/100,"Strain ϵ [%]"},0,2.8/100},
{{Y,2*(10^10),"Young Modulus"},(10^10),5*(10^10)},
{{σ0,(10^8),"Yiedling stress (Subscript[σ, y0])"},8*(10^7),5*(10^8)}
]

Implementation

3. Plastic flow equations

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

Implementation

3. Plastic flow equations

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

Where F=F(σ,(Subscript[σ, yo]),q)=(Subscript[σ, VM])-(Subscript[σ, y0])+α H

Implementation

3. Plastic flow equations

({{{, ({{(Subscript[ϵ, p])-(Subscript[ϵ, p,(n)])-(λ-(Subscript[λ, (n)]))((Sqrt[3])/(2 (Sqrt[(Subscript[J, 2])])))Dev[σ]=0
α-(Subscript[α, (n)])-(λ-(Subscript[λ, (n)]))=0
(1/Y)((Sqrt[3(Subscript[J, 2])])-(Subscript[σ, y0])+α H)=0,  }})}})

Where Dev[T]=T-(1/3)Tr[T] is the deviatoric contribution:

Implementation

3. Plastic flow equations

Stress evaluation

({{σ=(({{(Subscript[σ, 11]), (Subscript[σ, 12]), (Subscript[σ, 13])}, {(Subscript[σ, 21]), (Subscript[σ, 22]), (Subscript[σ, 23])}, {(Subscript[σ, 31]), (Subscript[σ, 32]), (Subscript[σ, 33])}}))⟶}})σ=σ[(Subscript[ϵ, el])]=σ[ϵ-(Subscript[ϵ, p])]

Implementation

3. Plastic flow equations

ϵ//MatrixForm

(σ=SolidMechanicsStress[,<|"PoissonRatio"->ν,"YoungModulus"->Y|>,ϵ])//MatrixForm

Implementation

3. Plastic flow equations

(Subscript[ϵ, tot])=(Subscript[ϵ, el])+(Subscript[ϵ, p])

ϵelastic=ϵ-ϵp;
ϵelastic//MatrixForm

(σ=SolidMechanicsStress[,<|"PoissonRatio"->ν,"YoungModulus"->Y|>,ϵelastic])//MatrixForm

Implementation

3. Plastic flow equations

({{{, ({{(Subscript[ϵ, p])-(Subscript[ϵ, p,(n)])-(γ-(Subscript[γ, (n)]))((Sqrt[3])/(2 (Sqrt[(Subscript[J, 2])])))Dev[σ]=0
α-(Subscript[α, (n)])-(γ-(Subscript[γ, (n)]))=0
(1/Y)((Subscript[σ, VM])-(Subscript[σ, y0])+α H)=0,  }})}})

Implementation

3. Plastic flow equations

eqns={
ϵp-ϵpn-(λ-λn)(3/(2 vonMisesStress[σ]))deviatoric[σ]==0,
α-αn-(λ-λn)==0,
(1/Y)(vonMisesStress[σ]-σy0-α*H)==0};

Implementation

3. Plastic flow equations

subEqns=Table[ϵp[[i,j]]-ϵpn[[i,j]]-(λ-λn)*Sqrt[3]*(Sqrt[3]/(2*vonMisesStress[σ]))*(deviatoric[σ])[[i,j]]==0,{i,1,3},{j,1,3}];

eqns=Flatten@{
Extract[subEqns,#]&/@{{1,1},{2,2},{3,3},{2,3},{1,3},{1,2}},
α-αn+(λ-λn)==0,
(1/Y)(vonMisesStress[σ]-σy0-α*H)==0};

Column[Style[#,Medium]&/@eqns]

Implementation

3. Plastic flow equations

ϵtotNumerics={{0.1,0,0},{0.1,0,0},{0.1,0,0}};
replacements=Flatten@{
Y->1.2*(10^9),
σy0->(10^6),
H->0,
Table[ϵpn[[i,j]]->0,{i,1,3},{j,1,3}],
Table[ϵ[[i,j]]->ϵtotNumerics[[i,j]],{i,1,3},{j,1,3}],
αn->0,
λn->0,
ν->0.3
};

Implementation

3. Plastic flow equations

numericsEqns=FullSimplify[eqns/.replacements];
Column[Style[#,Medium]&/@numericsEqns]

Implementation

3. Plastic flow equations

sol=FindRoot[numericsEqns,{{ϵp11,0},{ϵp22,0},{ϵp33,0},{ϵp23,0},{ϵp13,0},{ϵp12,0},{α,0},{λ,0}}];

(Normal[ϵp])/.sol//MatrixForm

Implementation

4. Iterative solution

pars["InitialStrain"]=ϵp;

SolidMechanicsPDEComponent[vars,pars];

Implementation

Practical example:

Implementation

Solve elastic - problem

Evaluate Yielding state

Plastic flow equations

Iterative solution (come back to 1)

height=1;
width=2.5;
thickness=0.01;
mesh=ToElementMesh[Rectangle[{0,0},{width,height}],"MaxCellMeasure"->0.25,"MeshElementBlocks"->3];
Show[mesh["Wireframe"] ,Graphics[{
Red,{Dashed,Line[{{width(1.1),0},{width(1.1),1.2height}}]},Line[{{width+0.01,(1+0.2)height},{width+0.01,(1)height}}],Arrow[{{width+0.01,(1+0.1)height},{width(1.1),(1+0.1)height}}],
(*Thickness[0.02],*)
{FaceForm[Transparent],EdgeForm[Black],Black,
Triangle[{{-0.1 height,0},{0,0.1height},{-0.1 height,0.2 height}}],
Table[{Triangle[{{-0.1 height,0+i},{0,0.1height+i},{-0.1 height,0.2 height+i}}],
Circle[{-0.13 height,0.05height+i},0.025height],
Circle[{-0.13 height,0.15height+i},0.025height]
},{i,0,height(1-0.2),height/4}
],
Table[{Triangle[{{0+i,-0.1 height},{0.1height+i,0},{0.2 height+i,-0.1 height}}],
Circle[{0.05height+i,-0.13 height},0.025height],
Circle[{0.15height+i,-0.13 height},0.025height]
},{i,0,width(1-0.1),height/4}
]


}

}]]

Implementation

1. Solve elastic - problem

mesh=ToElementMesh[Rectangle[{0,0},{2.5,1}],"MaxCellMeasure"->0.5];
vars={{u[x,y],w[x,y]},{x,y}};
pars=<|"PoissonRatio"->0.3,"YoungModulus"->2.1*10^9,"Thickness"->1,"YieldStress"->8*(10^7),"Hardening"->0|>;

Implementation

1. Solve elastic - problem

(Subscript[Γ, disp])=SolidDisplacementCondition[x==2.5,vars,pars,<|"Displacement"->{0.1*2.5,None}|>];
(Subscript[Γ, rollerX])=SolidDisplacementCondition[x==0,vars,pars,<|"Displacement"->{0,None}|>];
(Subscript[Γ, rollerY])=SolidDisplacementCondition[y==0,vars,pars,<|"Displacement"->{None,0}|>];

Implementation

1. Solve elastic - problem

displacement=NDSolveValue[{SolidMechanicsPDEComponent[vars,pars]=={0,0},(Subscript[Γ, disp]),(Subscript[Γ, rollerX]),(Subscript[Γ, rollerY])},vars[[1]],vars[[-1]]∈mesh]

Implementation

1. Solve elastic - problem

Show[
mesh["Wireframe"["MeshElementStyle" -> Directive[EdgeForm[LightGray], FaceForm[LightGray]]]],
  ElementMeshDeformation[mesh, displacement, "ScalingFactor" -> 2]["Wireframe"["MeshElementStyle" -> EdgeForm[Blue]]]]
SwatchLegend[{LightGray, Blue}, {"Undeformed", "Elastic"}]

Implementation

2. Evaluate Yielding state

strain=SolidMechanicsStrain[vars,pars,displacement];
stress=SolidMechanicsStress[vars,pars,strain];
σVM=VonMisesStress[vars,pars,stress]

Implementation

2. Evaluate Yielding state

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

σVM=VonMisesStress[vars,pars,stress]

DensityPlot[σVM,{x,y}∈mesh,]

Implementation

2. Evaluate Yielding state

σVMaverage=NIntegrate[σVM,{x,y}∈mesh]/NIntegrate[1,{x,y}∈mesh]

pars["YieldStress"]

GraphicsRow[{Show[
Plot[pars["YoungModulus"] ϵ,{ϵ,0,0.2},PlotLabels->{"Hooke's law"},PlotRange->{All,{(10^1),All}},AxesLabel->{"ϵ","σ"}],
ListPlot[{{0.1,σVMaverage}},PlotStyle->Directive[PointSize[Large],Blue]],
ListPlot[{{pars["YieldStress"]/pars["YoungModulus"],pars["YieldStress"]}},PlotStyle->Directive[PointSize[Large],Black],PlotLabels->Placed["Yielding point",Above]],ListLinePlot[{{0,pars["YieldStress"]},{0.2,pars["YieldStress"]}},PlotStyle->Dashed,PlotLabels->{"Perfect plasticity"}]],
SwatchLegend[{ Blue, Red}, { "Elastic", "Elastoplastic"}]}]

Implementation

3. Plastic flow equations

In a moment ...

Implementation

4. Iterative solution

While[data["Δ"]<(10^(-8)),

{plasticStrain,λ,α,data}=plasticStrainStepCompiled[vars,updatePars,epDispl,data];

updatePars=ReplacePart[updatePars,
{"InitialStrain"->plasticStrain,
"PlasticMultiplier"->λ,
"PlasticAlpha"->α}];

epDispl=NDSolveValue[eqn8[vars,updatePars],vars[[1]],vars[[-1]]∈mesh];

]

Implementation

4. Iterative solution

epDispl = {InterpolatingFunction[][x,y],InterpolatingFunction[][x,y]}

Implementation

GraphicsRow[{Show[
mesh["Wireframe"["MeshElementStyle" -> Directive[EdgeForm[LightGray], FaceForm[LightGray]]]],
  ElementMeshDeformation[mesh, displacement, "ScalingFactor" -> 2]["Wireframe"["MeshElementStyle" -> EdgeForm[Blue]]]
 ,ElementMeshDeformation[mesh, epDispl, "ScalingFactor" -> 2]["Wireframe"["MeshElementStyle" -> EdgeForm[Red]]]],
SwatchLegend[{LightGray, Blue, Red}, {"Undeformed", "Elastic", "Elastoplastic"}]}]

Implementation

σVM2=InterpolatingFunction[][x,y];
σVMplastic=NIntegrate[σVM2,{x,y}∈mesh]/NIntegrate[1,{x,y}∈mesh];

Show[
Plot[pars["YoungModulus"] ϵ,{ϵ,0,0.2},PlotLabels->{"Hooke's law"},PlotRange->{All,{(10^1),All}},AxesLabel->{"ϵ","σ"}],ListPlot[{{0.1,σVMplastic}},PlotStyle->Directive[PointSize[Large],Red]],
ListPlot[{{0.1,σVMaverage}},PlotStyle->Directive[PointSize[Large],Blue]],
ListPlot[{{pars["YieldStress"]/pars["YoungModulus"],pars["YieldStress"]}},PlotStyle->Directive[PointSize[Large],Black],PlotLabels->Placed["Yielding point",Above]],ListLinePlot[{{0,pars["YieldStress"]},{0.2,pars["YieldStress"]}},PlotStyle->Dashed,PlotLabels->{"Perfect plasticity"}],
Graphics[Arrow[{{0.1,σVMaverage},{0.1,σVMplastic}}]]]
SwatchLegend[{ Blue, Red}, { "Elastic", "Elastoplastic"}]

Implementation

Discretized solution

({{{, ({{(Subscript[ϵ, p])-(Subscript[ϵ, p,(n)])-(γ-(Subscript[γ, (n)]))((Sqrt[3])/(2 (Sqrt[(Subscript[J, 2])])))Dev[σ]=0
α-(Subscript[α, (n)])-(γ-(Subscript[γ, (n)]))=0
(1/Y)((Sqrt[3(Subscript[J, 2])])-(Subscript[σ, y0])+α H)=0,  }})}})

Discretized solution

({{{, ({{(Subscript[ϵ, p])-(Subscript[ϵ, p,(n)])-(γ-(Subscript[γ, (n)]))((Sqrt[3])/(2 (Sqrt[(Subscript[J, 2])])))Dev[σ]=0
α-(Subscript[α, (n)])-(γ-(Subscript[γ, (n)]))=0
(1/Y)((Sqrt[3(Subscript[J, 2])])-(Subscript[σ, y0])+α H)=0,  }})}})

{displU,displV}=displacement

Head[displU]["ValuesOnGrid"]

Implementation

mesh=ToElementMesh[Rectangle[],"MaxCellMeasure"->0.5];
mesh["Wireframe"]

Implementation

order=2;
intOrder=4;
isf=IntegratedShapeFunction[QuadElement,order,intOrder];
coords=mesh["Coordinates"];
nodefun=Map[ElementMeshInterpolation[mesh,#]&,IdentityMatrix[Length[mesh["Coordinates"]]]];
eleCoords=GetElementCoordinates[coords,ElementIncidents[mesh["MeshElements"]][[1]]];
mappedCoords=(isf[[All,1]].#)&/@eleCoords;
mappedCoords3D=(Append[#,0])&/@Flatten[mappedCoords,1][[1;;8]];
coords3D=(Append[#,0]&/@mesh["Coordinates"][[#[[{1,2,3,4}]]]])&/@mesh["MeshElements"][[1,1]];
nodeFunIndex={1,2,5,4,13,10};
Grid[Partition[
Table[
Show[
Plot3D[nodefun[[nodeFunIndex[[i]]]][x,y],Element[{x,y},mesh],Mesh->False,PlotRange->All]
,
Graphics3D[{Red,Point[{#[[1]],#[[2]],0}]&/@Flatten[mappedCoords,1], FaceForm[],EdgeForm[Black],Polygon[#]&/@coords3D,
Black,Point[#]&/@coords3D,
Line[{#,Flatten@{#[[1]],#[[2]],nodefun[[nodeFunIndex[[i]]]][#[[1]],#[[2]]]}}]&/@mappedCoords3D
,Blue,
Point[Flatten@{#[[1]],#[[2]],nodefun[[nodeFunIndex[[i]]]][#[[1]],#[[2]]]}]&/@mappedCoords3D}]],{i,1,Length[nodeFunIndex]}]
,3]]

Implementation

Discretization

Graphics[{
Line[{{0,0},{1,1}}],
Point[Table[{i,i},{i,0,1,0.1}]],
Arrow[BezierCurve[{{0.55,0.5},{1,0.2},{1.5,0.4}}]],
Line[{{2,0.5},{4,0.5}}],
PointSize[Large],
Point[{{2,0.5},{3,0.5},{4,0.5}}]
}]

Implementation

Needs["NumericalDifferentialEquationAnalysis`"]
gp=GaussianQuadratureWeights[2,-1,1]

pts={x1->-1,x2->0,x3->1};
phi[x_,x1_,x2_,x3_]:=(((x2-x)(x3-x))/((x2-x1)(x3-x1)));
GraphicsRow[{
Show[Plot[{(((x2-x)(x3-x))/((x2-x1)(x3-x1)))/.pts,
(((x1-x)(x3-x))/((x1-x2)(x3-x2)))/.pts,
(((x1-x)(x2-x))/((x1-x3)(x2-x3)))/.pts
}
,{x,-1,1},Filling->{1->Axis},Axes->{True,False}],Graphics[{PointSize[Large],Point[{{-1,0},{0,0},{1,0}}]}]],Show[
Plot[phi[x,-1,0,1],{x,-1,1},Axes->{True,False}],
ListPlot[{#[[1]],phi[#[[1]],-1,0,1]}&/@gp,PlotStyle->PointSize[Large],Filling->Axis]
,Graphics[{PointSize[Large],Point[{{-1,0},{0,0},{1,0}}]}]]}]

Implementation

u[x_]:=0.1*(x^2)(3*10)
du=D[u[x],x];
uDiscrete=Table[{i,u[i]},{i,1,10,2}];
udiff=Flatten[Table[{{uDiscrete[[i,1]],(-uDiscrete[[i,2]]+uDiscrete[[i+1,2]])/(uDiscrete[[i+1,1]]-uDiscrete[[i,1]])},{uDiscrete[[i+1,1]],(-uDiscrete[[i,2]]+uDiscrete[[i+1,2]])/(uDiscrete[[i+1,1]]-uDiscrete[[i,1]])}},{i,1,Length[uDiscrete]-1}],1];
udiffTrue=Prepend[udiff,{udiff[[1,1]],0}];
PrependTo[udiffTrue,{0,0}];
GraphicsColumn[{GraphicsRow[{
Show[
Plot[u[x],{x,0,9},AxesLabel->{"x","u[x]"},PlotLabel->"displacement",PlotStyle->Directive[Dashed,Brown]],
ListLinePlot[uDiscrete,PlotStyle->Directive[Orange]],
ListPlot[uDiscrete,PlotStyle->Directive[PointSize[Large],Black]
]
],
Show[
Plot[du,{x,0,10},AxesLabel->{"x","u'[x]"},PlotLabel->"\!\(\*FractionBox[\(∂\\ displacement\), \(∂x\)]\)",PlotStyle->Directive[Dashed,Blue]],
ListLinePlot[udiffTrue,PlotStyle->Directive[Gray]],
ListPlot[Table[{udiffTrue[[i,1]]+(udiffTrue[[i+1,1]]-udiffTrue[[i,1]])/2,du/.x->(udiffTrue[[i,1]]+(udiffTrue[[i+1,1]]-udiffTrue[[i,1]])/2)},{i,3,Length[udiffTrue]-1,2}],PlotStyle->{PointSize[Large],Red}],
ListPlot[Table[{k,du/.{x->k-1}},{k,1,Length[udiffTrue]-1,2}],PlotStyle->{PointSize[Large],Black}]
]
}]
,
GraphicsRow[{
LineLegend[{ {Dashed,Brown},Orange,{Dashed,Blue},Gray}, { "f[x]", "discretized f[x]","f'[x]","discretized f'[x]"}],
PointLegend[{ Black,Red}, { "Nodal evaluation","Superconvergence points"}]}]
}]

Implementation

order=2;
intOrder=4;
isf=IntegratedShapeFunction[QuadElement,order,intOrder];
sf=ElementShapeFunction[TriangleElement,order];
(*sf[r,s]*)
ip=ElementIntegrationPoints[TriangleElement,intOrder];
isf[[All,1]]===(sf@@@ip);
coords=mesh["Coordinates"];
eleCoords=GetElementCoordinates[coords,ElementIncidents[mesh["MeshElements"]][[1]]];
mappedCoords=(isf[[All,1]].#)&/@eleCoords;
Show[Graphics[{PointSize[Large],Red,Point/@mappedCoords}],mesh["Wireframe"]]

Iterative solution

Iterative solution

While[data["Δ"]<(10^(-8)),

{plasticStrain,λ,α,data}=plasticStrainStepCompiled[vars,updatePars,epDispl,data];

updatePars=ReplacePart[updatePars,
{"InitialStrain"->plasticStrain,
"PlasticMultiplier"->λ,
"PlasticAlpha"->α}];

epDispl=NDSolveValue[eqn8[vars,updatePars],vars[[1]],vars[[-1]]∈mesh];

]

Compiled functions

compiledNR[cf_,cjf_,cyield_]:=Compile[{{αp0,_Real,2},{λp0,_Real,2},{ϵp0,_Real,2},{ϵtot,_Real,2},{αpn,_Real,2},{λpn,_Real,2},{ϵpn,_Real,2}},
Module[{},
nsteps=50;
xNR=Flatten[{αp0[[1,1]],λp0[[1,1]],ϵp0[[1,1]],ϵp0[[2,2]],ϵp0[[3,3]],ϵp0[[2,3]],ϵp0[[1,3]],ϵp0[[1,2]]}];

Do[
xn=xNR;
αp=xn[[1]]; λp=xn[[2]];
ϵp=(({{xn[[1+2]], xn[[6+2]], xn[[5+2]] }, {xn[[6+2]], xn[[2+2]], xn[[4+2]]}, {xn[[5+2]], xn[[4+2]], xn[[3+2]]}}));

f=cf[ϵtot,αp,λp,ϵp,αpn[[1,1]],λpn[[1,1]],ϵpn];
jf=cjf[ϵtot,αp,λp,ϵp,αpn[[1,1]],λpn[[1,1]],ϵpn];

Δx=-LinearSolve[jf,f];

xNR=xn+Δx;
If[Norm[Δx]<tol, Break[]; ]

,{i,1,nsteps}];

xNR
]]

Compiled functions

cfComp[vars_,pars_]:=Block[{ϵtot,αp,λp,ϵp,αpn,λpn,ϵpn},
With[{body=plasticFlowSimbolic[vars,pars][[1]]},

cf=Compile[{{ϵtot,_Real,2},{αp,_Real},{λp,_Real},{ϵp,_Real,2},{αpn,_Real},{λpn,_Real},{ϵpn,_Real,2}},

Block[{},
{ϵp11,ϵp22,ϵp33,ϵp23,ϵp13,ϵp12}={ϵp[[1,1]],ϵp[[2,2]],ϵp[[3,3]],ϵp[[2,3]],ϵp[[1,3]],ϵp[[1,2]]};

{ϵpn11,ϵpn22,ϵpn33,ϵpn23,ϵpn13,ϵpn12}={ϵpn[[1,1]],ϵpn[[2,2]],ϵpn[[3,3]],ϵpn[[2,3]],ϵpn[[1,3]],ϵpn[[1,2]]};

{ϵtot11,ϵtot22,ϵtot33,ϵtot23,ϵtot13,ϵtot12}={ϵtot[[1,1]],ϵtot[[2,2]],ϵtot[[3,3]],ϵtot[[2,3]],ϵtot[[1,3]],ϵtot[[1,2]]};

body
]]];
cf
]

Compiled functions

plasticFlowSimbolic[vars_,parsIN_]:=Module[{},
pVars=Flatten[{αp,λp,ϵpVoigth[[idx]]}];

σ=applyStiffMatrix8[vars,pars,ϵtot[[1;;dim,1;;dim]]]-applyStiffMatrix8[vars,pars,ϵp[[1;;dim,1;;dim]]];

σy0=pars["YieldStress"]; H=pars["Hardening"];  σVM=vmStress8[σ];
σdev=PDEModels`MatrixToVoigtVector[deviatoric8[σ]];

subsys=
Table[(ϵpVoigth[[idx[[i]]]]-PDEModels`MatrixToVoigtVector[ϵpn][[idx[[i]]]]-(λp-λpn)*(3/(2*σVM))*σdev[[idx[[i]]]]),{i,1,Length[idx]}];
Y=pars["YoungModulus"];
dFdq=1;
	f=Flatten@{
	subsys,
	αp-αpn +(λp-λpn)*dFdq,
	(1/Y)(σVM-σy0+(αp)*H)
};

jf=Transpose[D[f,#]&/@pVars];

{f,jf}
]

Plastic flow step

(* [...] *)
globalSol=compiledNR[
			Flatten[αpguess,{1,3,4}],
			Flatten[λpguess,{1,3,4}],
			Transpose[#,{1,3,4,2}]&/@ϵpguess,
			Transpose[#,{1,3,4,2}]&/@ϵtot,
			Flatten[αpn,{1,3,4}],
			Flatten[λpn,{1,3,4}],
			Transpose[#,{1,3,4,2}]&/@ϵpn
];

(* reshape *)
(* [...] *)

(* return: *)
{plasticStrain,λ,α,data}

Examples

Cook’s membrane

Cook's membrane

bmesh=ToBoundaryMesh[Polygon[{p1={0,0},p2={48,44},p3={48,44+16},p4={0,44}}],"MaxCellMeasure"->5];

Show[
bmesh["Wireframe"],
Graphics[{Red,
Table[Arrow[{{48+1,i},{48+1,i+2.5}}],{i,44,44+16,3}],
{Black,Table[Line[{{-3,i-2},{0,i}}],{i,0,44,2}]},
Black,Thickness[0.01],Line[{{0,0},{0,44}}],
(*PointSize[Large],*)Point[{p1,p2,p3,p4}],
Gray,Thin,
Line[{{0,0},{48,0}}],
Line[{{48,2},{48,-2}}],
Text["48",{24,-2}],
Line[{{-4,0},{-4,44}}],
Line[{{-6,44},{-2,44}}],
Text["44",{-6,22}],
Line[{{48+4,44},{48+4,44+16}}],
Line[{{48+4-2,44},{48+4+2,44}}],
Line[{{48+4-2,44+16},{48+4+2,44+16}}],
Text["16",{44+10,44+8}]
}]
]

Cook's membrane

vars={{u[x,y],v[x,y]},{x,y}};
mesh=ToElementMesh[Polygon[{{0,0},{48,44},{48,44+16},{0,44}}],"MaxCellMeasure"->10];
mesh["Wireframe"]

Cook's membrane

eqn8[vars_,pars_]:=Module[{Γforce,Γwall},

Γforce=SolidBoundaryLoadValue[x==48,vars,pars,<|"Force"->{0,1500}|>];
Γwall=SolidDisplacementCondition[x==0,vars,pars,<|"Displacement"->{0,0}|>];

{SolidMechanicsPDEComponent[vars,pars]==Γforce,Γwall}

]

Cook's membrane

pars=<|"YoungModulus"->7*(10^3),
"PoissonRatio"->0.3,
"Thickness"->1,
"YieldStress"->2.5*(10^2),
"Hardening"->0,
"SolidMechanicsModelForm"->"PlaneStrain"|>;

Cook's membrane

elasticDispl=NDSolveValue[eqn8[vars,pars],vars[[1]],vars[[-1]]∈mesh]

Cook's membrane

Show[mesh["Wireframe"],  ];

Scale factor = 3

Cook's membrane

strain=SolidMechanicsStrain[vars,pars,elasticDispl];
stress=SolidMechanicsStress[vars,pars,strain];
stress=PDEModels`SolidMechanicsExtendedStress[vars,pars,displ,strain,stress];
vonMises=VonMisesStress[vars,pars,stress];

Cook's membrane

vonMises = InterpolatingFunction[][x,y]

Cook's membrane

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

Cook's membrane

yieldingFuncQ=If[vonMises<pars["YieldStress"],0,1];

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

Cook's membrane

While[data["Δ"]>tolerance,

{plasticStrain,λ,α,data}=plasticStrainStepCompiled[vars,updatePars,epDispl,data];

updatePars=ReplacePart[updatePars,{"InitialStrain"->plasticStrain,"PlasticMultiplier"->λ,"PlasticAlpha"->α}];

epDispl=NDSolveValue[eqn8[vars,updatePars],vars[[1]],vars[[-1]]∈mesh];

];

Cook's membrane

While[data["Δ"]>tolerance,

(* [... ] *)

];

epDispl =

Cook's membrane

sc=3;
Show[mesh[],
  ElementMeshDeformation[mesh, elasticDispl, "ScalingFactor" -> sc][],
ElementMeshDeformation[mesh, epDispl, "ScalingFactor" -> sc]["Wireframe"["MeshElementStyle" -> EdgeForm[Red]]],Graphics[{Red,
Table[Arrow[{{48+1,i},{48+1,i+2.5}}],{i,44,44+16,3}],
Black,Thickness[0.02],Line[{{0-0.6,0},{0-.6,44}}]}]]

Scale factor = 3

Cook's membrane

Current implementation

Future implementation

Integrating several Yielding criteria

Implementation of Kinematic Hardening

Implementation of different hardening laws

Needs["NDSolve`FEM`"]
deviatoric[T_]:=T-(1/3)Tr[T]IdentityMatrix[Length[T]];
vonMisesStress[σ_]:=Sqrt[( (((σ[[1,1]]-σ[[2,2]]))^2)+ (((σ[[2,2]]-σ[[3,3]]))^2)+(((σ[[3,3]]-σ[[1,1]]))^2)+6(((σ[[1,2]])^2)+((σ[[1,3]])^2)+((σ[[2,3]])^2)))/2];

ϵ=SymmetrizedArray@PDEModels`VoigtVectorToMatrix[{ϵ11,ϵ22,ϵ33,ϵ23,ϵ13,ϵ12}];
ϵp=SymmetrizedArray@PDEModels`VoigtVectorToMatrix[{ϵp11,ϵp22,ϵp33,ϵp23,ϵp13,ϵp12}];
ϵpn=SymmetrizedArray@PDEModels`VoigtVectorToMatrix[{ϵpn11,ϵpn22,ϵpn33,ϵpn23,ϵpn13,ϵpn12}];
ϵel=SymmetrizedArray@PDEModels`VoigtVectorToMatrix[{ϵel11,ϵel22,ϵel33,ϵel23,ϵel13,ϵel12}];