# Stress Operators for Finite Element Analysis

For stress analysis in the finite element method you need a stress operator. For two dimensional plane stress it can be found here. For two dimensional plain strain it can be found here. I need the two dimensional case where there is axial symmetry so that in cylindrical coordinates this is displacements and stress in terms of {r,z}. Getting these stress operators is not so straight forward since one needs a careful application of Inactive. I also need the calculation of stress from the displacement fields. Can you help please?

For stress in 3D the stress operator may be found here.

Here are the stress operators and the method for calculating stress for plane stress and plane strain. I also give vonMises stress for plane stress

ClearAll[planeStressOperator, Y, ν];
planeStressOperator[{u_, v_}, {x_,
y_}, {Y_, ν_}] := {Inactive[
Div][({{0, -((Y ν)/(1 - ν^2))}, {-((Y (1 - ν))/(2 (1 \
- ν^2))), 0}}.Inactive[Grad][v[x, y], {x, y}]), {x, y}] +
Inactive[
Div][({{-(Y/(1 - ν^2)),
0}, {0, -((Y (1 - ν))/(2 (1 - ν^2)))}}.Inactive[Grad][
u[x, y], {x, y}]), {x, y}],
Inactive[Div][({{0, -((Y (1 - ν))/(2 (1 - ν^2)))}, {-((Y \
ν)/(1 - ν^2)), 0}}.Inactive[Grad][u[x, y], {x, y}]), {x, y}] +
Inactive[
Div][({{-((Y (1 - ν))/(2 (1 - ν^2))),
0}, {0, -(Y/(1 - ν^2))}}.Inactive[Grad][
v[x, y], {x, y}]), {x, y}]};

ClearAll[planeStrainOperator];
planeStrainOperator[{u_, v_}, {x_,
y_}, {Y_, ν_}] := {Inactive[
Div][({{0, -((Y ν)/((1 -
2 ν) (1 + ν)))}, {-(Y/(2 (1 + ν))),
0}}.Inactive[Grad][v[x, y], {x, y}]), {x, y}] +
Inactive[
Div][({{-((Y (1 - ν))/((1 - 2 ν) (1 + ν))),
0}, {0, -(Y/(2 (1 + ν)))}}.Inactive[Grad][
u[x, y], {x, y}]), {x, y}],
Inactive[
Div][({{0, -(Y/(2 (1 + ν)))}, {-((Y ν)/((1 -
2 ν) (1 + ν))), 0}}.Inactive[Grad][
u[x, y], {x, y}]), {x, y}] +
Inactive[
Div][({{-(Y/(2 (1 + ν))),
0}, {0, -((Y (1 - ν))/((1 -
2 ν) (1 + ν)))}}.Inactive[Grad][
v[x, y], {x, y}]), {x, y}]};

ClearAll[stress2D]
stress2D[{uif_InterpolatingFunction,
vif_InterpolatingFunction}, {Y_, ν_}, ps_] :=
Block[{fac, dd, df, mesh, coords, dv, ux, uy, vx, vy, ex, ey, gxy,
sxx, syy, sxy},
dd = Outer[(D[#1[x, y], #2]) &, {uif, vif}, {x, y}];
Which[
ps == "PlaneStress",
fac = Y/(1 - ν^2)*{{1, ν, 0}, {ν, 1, 0}, {0,
0, (1 - ν)/2}},
ps == "PlaneStrain",
fac = Y/((1 + ν) (1 - 2 ν)) *{{1 - ν, ν,
0}, {ν, 1 - ν, 0}, {0, 0, (1 - 2 ν)/2}}
];

df = Table[Function[{x, y}, Evaluate[dd[[i, j]]]], {i, 2}, {j, 2}];
(*the coordinates from the ElementMesh*)
mesh = uif["Coordinates"][];
coords = mesh["Coordinates"];
dv = Table[df[[i, j]] @@@ coords, {i, 2}, {j, 2}];
ux = dv[[1, 1]];
uy = dv[[1, 2]];
vx = dv[[2, 1]];
vy = dv[[2, 2]];
ex = ux;
ey = vy;
gxy = (uy + vx);
sxx = fac[[1, 1]]*ex + fac[[1, 2]]*ey;
syy = fac[[2, 1]]*ex + fac[[2, 2]]*ey;
sxy = fac[[3, 3]]*gxy;
{ElementMeshInterpolation[{mesh}, sxx],
ElementMeshInterpolation[{mesh}, syy],
ElementMeshInterpolation[{mesh}, sxy]
}]

ClearAll[vonMisesPlaneStress];
vonMisesPlaneStress[{sxx_, syy_, sxy_}] :=
Function[{x, y}, Sqrt[
sxx[x, y]^2 - sxx[x, y] syy[x, y] + syy[x, y]^2 + 3 sxy[x, y]^2]]


Here is an example of their use. I am looking at stress concentrations near a cut-out. I include a MeshRefinementFunction at the cut-out which I am using to show that improving the mesh does not help with the stress concentration in this case.

Needs["NDSolveFEM"];
L = 2; (* Length *)
h = 1; (* Height *)
Y = 10^3; (* Modulus of elasticity *)
ν = 33/100;  (* Poisson ratio *)
stress = 10;  (* stress *)
gls = 0.1; (* Grid length scale *)
nr = 10;  (* Factor for refinement *)
r0 = 2 gls;  (* radius for refinement *)

ClearAll[cf];
cf = Compile[{{c, _Real, 2}, {a, _Real, 0}},
Block[{com, c1, c2, r1, r2, r},
com = Total[c]/3;
c1 = {0.2 L, 0.6 h};
c2 = {0.4 L, 0.6 h};
r1 = Norm[com - c1];
r2 = Norm[com - c2];
r = Min[{r1, r2}];

If[r < r0 && a > (Sqrt/4) (gls/nr + r/r0 gls)^2, True, False]

]
];

mesh = ToElementMesh[
RegionDifference[Rectangle[{0, 0}, {L, h}],
RegionUnion[Rectangle[{0.2 L, 0.6 h}, {0.4 L, h}],
Disk[{0.75 L, 0.4 h}, 0.2 h]]],
MaxCellMeasure -> {"Length" -> .1}, MeshRefinementFunction -> cf];

Show[mesh["Wireframe"],
Graphics[{
{Thickness[0.01], Line[{{0, -0.2 h}, {0, 1.2 h}}]},
Table[Arrow[{{L, y}, {L + L/2 (h - y), y}}], {y, 0, h, h/12}],
Table[Arrow[{{x, 3 h/2}, {x, h}}], {x, 3 L/4, L, L/20}]
}]
] I can solve for plane stress and work out the displacements as follows.

{uif, vif} = NDSolveValue[{
planeStressOperator[{u, v}, {x, y}, {Y, ν}] == {
NeumannValue[stress (h - y), x == L && 0 <= y <= h],
NeumannValue[-stress, 3 L/4 <= x <= L && y == h]},
DirichletCondition[u[x, y] == 0, x == 0],
DirichletCondition[v[x, y] == 0, x == 0]
}, {u, v}, {x, y} ∈ mesh];

dmesh = ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 1];
ContourPlot[uif[x, y], {x, y} ∈ mesh,
AspectRatio -> Automatic, ColorFunction -> "Temperature",
PlotRange -> All]
Show[{
mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]  Finally I calculate the stresses and the von Mises stress.

{σxx, σyy, σxy} =
stress2D[{uif, vif}, {Y, ν}, "PlaneStress"];
vm = vonMisesPlaneStress[{σxx, σyy, σxy}];

Plot3D[vm[x, y], {x, y} ∈ mesh, BoxRatios -> {L, h, h},
PlotRange -> All] All this shows that the finite element method is looking to be in good shape. What I need to complete my stress analysis is the axisymmetric stress operator for {r, z} and a method of calculating the stresses. Has anyone done this?

I have found a good reference for this problem. To help get started here are the strain-displacement and strain-stress equations for three dimensions { r, θ, z}.

$\epsilon _{\text{rr}}=\frac{\partial u_r}{\partial r}$

$\epsilon _{\theta \theta }=\frac{\partial u_{\theta }}{r\partial \theta }+\frac{u_r}{r}$

$\epsilon _{\text{r$\theta $}}=\frac{\partial u_{\theta }}{\partial r}-\frac{u_{\theta }}{r}+\frac{\partial u_r}{r\partial \theta }$

$\epsilon _{\text{zz}}=\frac{\partial u_z}{\partial z}$

$\epsilon _{\text{$\theta $z}}=\frac{\partial u_{\theta }}{\partial z}+\frac{\partial u_z}{r\partial \theta }$

$\epsilon _{\text{rz}}=\frac{\partial u_r}{\partial z}+\frac{\partial u_z}{\partial r}$

$Y \epsilon _{\text{rr}}=\sigma _{\text{rr}}-\nu \left(\sigma _{\theta \theta }+\sigma _{\text{zz}}\right)$

$Y \epsilon _{\theta \theta }=\sigma _{\theta \theta }-\nu \left(\sigma _{\theta \theta }+\sigma _{\text{zz}}\right)$

$Y \epsilon _{\text{zz}}=\sigma _{\text{zz}}-\nu \left(\sigma _{\theta \theta }+\sigma _{\text{rr}}\right))$

$\tau _{\text{r$\theta $}}=G \epsilon _{\text{r$\theta $}}$

$\tau _{\text{rz}}=G \epsilon _{\text{rz}}$

$\tau _{\text{$\theta $z}}=G \epsilon _{\text{$\theta $z}}$

Here Y, G and ν are the modulus of elasticiy, the shear modulus and Poisson ratio respectively.

Putting these into Mathematica format gives

strainDisplacement = {
err[r, θ, z] == D[ur[r, θ, z], r],
eθθ[r, θ, z] ==
D[uθ[r, θ, z], θ]/r + ur[r, θ, z]/r,
erθ[r, θ, z] ==
D[uθ[r, θ, z], r]/r - uθ[r, θ, z]/r +
D[ur[r, θ, z], θ]/r,
ezz[r, θ, z] == D[uz[r, θ, z], z],
eθz[r, θ, z] ==
D[uθ[r, θ, z], z] + D[uz[r, θ, z], θ]/
r,
erz[r, θ, z] ==
D[ur[r, θ, z], z] + D[uz[r, θ, z], z] +
D[uz[r, θ, z], r]
};
strainStress = {
Y err[r, θ, z] == σrr[r, θ,
z] - ν (σθθ[r, θ,
z] + σzz[r, θ, z]),
Y eθθ[r, θ, z] == σθθ[
r, θ,
z] - ν (σzz[r, θ,
z] + σθθ[r, θ, z]),
Y ezz[r, θ, z] == σzz[r, θ,
z] - ν (σrr[r, θ,
z] + σθθ[r, θ, z]),
τrθ[r, θ, z] == G erθ[r, θ, z],
τrz[r, θ, z] == G erz[r, θ, z],
τθz[r, θ, z] == G eθz[r, θ, z]
};


I guess that as I want a two dimensional axisymmetric {r,z} formulation all the derivatives with respect to θ will be zero. Not sure where to go then... Somehow we need to get the equations in the same form as above involving just derivatives of the two displacements ur and uz. Also we need the Inactive in the correct place.

• Hugh, I am confused whether you are asking about the fundamental physics, or how to implement them in Mathematica. Do you have an analytical expression / a starting point for the stress operator you are interested in? – MarcoB Apr 30 '18 at 19:11
• I need the mathematica implementation. I will look to see if I can find the analytical expression. – Hugh Apr 30 '18 at 19:12
• @MarcoB I have added the equations which I hope will give you a start. – Hugh Apr 30 '18 at 20:09