# 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))),
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 -
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"][[1]];
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[3]/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? Apr 30, 2018 at 19:11
• I need the mathematica implementation. I will look to see if I can find the analytical expression.
– Hugh
Apr 30, 2018 at 19:12
• @MarcoB I have added the equations which I hope will give you a start.
– Hugh
Apr 30, 2018 at 20:09
• I have a question how do you deal with the boundary condition of shear stress? Mar 10, 2020 at 1:52
• @PleaseCorrectGrammarMistakes You can put in tangential as well as normal stresses. However, watch the directions very carefully.
– Hugh
Mar 10, 2020 at 12:51

This question has been around unanswered for quite a while. In the mean time much has happened and the SolidMechanicsPDEComponent can be used to set up axisymmetric solid mechanics models. I'll show the model first, and then I'll add a few comments at the end.

Create the geometry:

Needs["NDSolveFEM"];
geometricParameters = <|L -> 2, h -> 1|>;
region =
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]]] /. geometricParameters;


Create the mesh:

(meshAxi = ToElementMesh[region])["Wireframe"]


varsAxi = {{u[r, z], 0, w[r, z]}, {r, \[Theta], z}};
parsAxi =
Join[<|"RegionSymmetry" -> "Axisymmetric", "YoungModulus" -> 10^4,
"PoissonRatio" -> 33/100|>, geometricParameters];


Note the specific way I set up the axisymmetric case in that I used 0 as the displacement in the azimuthal direction and theta as a independent variable. There are other ways to set this up. But this allows us to get the 0 displacement in the v-direction.

Solve the PDE:

displacementAxi =
NDSolveValue[{SolidMechanicsPDEComponent[varsAxi, parsAxi] ==
parsAxi, <|"Pressure" -> {10  (h - z), 0, 0}|>] +
SolidBoundaryLoadValue[z == h && 3 L/4 <= r <= L, varsAxi,
parsAxi, <|"Pressure" -> {0, 0, -10}|>]
, SolidDisplacementCondition[z == 0 && r < L/10, varsAxi,
parsAxi, <|"Displacement" -> {0, None, 0}|>]
, SolidDisplacementCondition[z == h && r < L/10, varsAxi,
parsAxi, <|"Displacement" -> {0, None, 0}|>]
}, varsAxi[[1]], {r, z} \[Element] meshAxi
]


Visualize the solution:

Show[
meshAxi["Wireframe"["MeshElement" -> "BoundaryElements"]],
ElementMeshDeformation[meshAxi, displacementAxi[[{1, 3}]],
"ScalingFactor" -> 5]["Wireframe"]]


Note that I changed the boundary conditions a bit (and the Young modulus value). For an axisymmetric case you need to be careful with the boundary conditions at r==0 since that is actually within the rotational region and you can not really clamp the object there. So I choose to hold it fixed in a small part at the top and bottom. Hope that makes sense.

Compute and visualize the vonMises stress:

strain = SolidMechanicsStrain[varsAxi, parsAxi, displacementAxi];
stress =
SolidMechanicsStress[varsAxi, parsAxi, strain, displacementAxi];
vms = Head[PDEModelsVonMisesStress[varsAxi, parsAxi, stress]];
ContourPlot[vms[r, z], {r, z} \[Element] meshAxi, PlotRange -> All,
`

AspectRatio -> Automatic, Contours -> 25]