# Stress analysis using finite element

Is there a way to automatically extrapolate the stress distribution? I'm analysing a disk under compressive force. I only have found the displacement but I need the stress distribution.

Needs["NDSolveFEM"];

Reg1 = Annulus[{0, 0}, {14, 25}];

mesh = ToElementMesh[Reg1, "MaxCellMeasure" -> 0.6,
"MaxBoundaryCellMeasure" -> 0.09];
hld = {\[Sigma]x[x, y] ==
Y/(1 - \[Nu]^2) (D[u[x, y], x] + \[Nu] u[x, y]/x), \[Sigma]y[x,
y] == Y/(1 - \[Nu]^2) (u[x, y]/x + \[Nu] D[u[x, y], x])};
planeStress = {Inactive[
Div][{{0, -((Y \[Nu])/(1 - \[Nu]^2))}, {-((Y (1 - \[Nu]))/(2 \
(1 - \[Nu]^2))), 0}}.Inactive[Grad][v[x, y], {x, y}], {x, y}] +
Inactive[
Div][{{-(Y/(1 - \[Nu]^2)),
0}, {0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}}.Inactive[
Grad][u[x, y], {x, y}], {x, y}],
Inactive[
Div][{{0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}, {-((Y \
\[Nu])/(1 - \[Nu]^2)), 0}}.Inactive[Grad][u[x, y], {x, y}], {x, y}] +
Inactive[
Div][{{-((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2))),
v[x, y], {x, y}], {x, y}]} /. {Y -> 3416, \[Nu] -> 33/100};
{uif, vif, \[Sigma]x1, \[Sigma]y1} =
NDSolveValue[{planeStress == {0, NeumannValue[-30, y >= 24.5]}, hld,
{DirichletCondition[u[x, y] == 0, y < -24.5],
DirichletCondition[v[x, y] == 0, y < -24.5]}}, {u,
v, \[Sigma]x, \[Sigma]y}, {x, y} \[Element] mesh];
dmesh = ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 1];

• Is there something wrong with your solution??? I mean, with your code and defining something for Young's modulus, Y  and Poisson's raion, nu. But you have to use them also in your Hooke's law. You can plot your solved stress, e.g., with Plot3D[\[Sigma]x1[x, y], Element[{x, y}, Reg1]]. Dec 16, 2016 at 9:59

Here is a way how to recover the stress from the distribution. Upfront, please double check that the math is correct.

Needs["NDSolveFEM"];

Y = 1000;
\[Nu] = 1/3;

Reg1 = Annulus[{0, 0}, {14, 25}];

mesh = ToElementMesh[Reg1, "MaxCellMeasure" -> 1];
(*hld={\[Sigma]x[x,y]\[Equal]Y/(1-\[Nu]^2) (D[u[x,y],x]+\[Nu] \
u[x,y]/x),\[Sigma]y[x,y]\[Equal]Y/(1-\[Nu]^2) (u[x,y]/x+\[Nu] \
D[u[x,y],x])};*)

planeStress = {Inactive[
Div][{{0, -((Y \[Nu])/(1 - \[Nu]^2))}, {-((Y (1 - \[Nu]))/(2 \
(1 - \[Nu]^2))), 0}}.Inactive[Grad][v[x, y], {x, y}], {x, y}] +
Inactive[
Div][{{-(Y/(1 - \[Nu]^2)),
0}, {0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}}.Inactive[
Grad][u[x, y], {x, y}], {x, y}],
Inactive[
Div][{{0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}, {-((Y \
\[Nu])/(1 - \[Nu]^2)), 0}}.Inactive[Grad][u[x, y], {x, y}], {x, y}] +
Inactive[
Div][{{-((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2))),
v[x, y], {x, y}], {x, y}]} /. {Y -> 3416, \[Nu] -> 33/100};
{uif, vif} =
NDSolveValue[{planeStress == {0,
NeumannValue[-30, y >= 24.5]}, {DirichletCondition[u[x, y] == 0,
y < -24.5], DirichletCondition[v[x, y] == 0, y < -24.5]}}, {u,
v}, {x, y} \[Element] mesh];
dmesh = ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 1];
dmesh["Wireframe"] Next we write a little helper function to compute the von Mises stress:

ClearAll[VonMisesStress]
VonMisesStress[{ufun_InterpolatingFunction,
vfun_InterpolatingFunction}, fac_] := Block[
{dd, df, mesh, coords, dv, ux, uy, vx, vy, ex, ey, gxy, sxx, syy,
sxy},
dd = Outer[(D[#1[x, y], #2]) &, {ufun, vfun}, {x, y}];
df = Table[Function[{x, y}, Evaluate[dd[[i, j]]]], {i, 2}, {j, 2}];
(* the coordinates from the ElementMesh *)

mesh = ufun["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,syy,sxy}*)
ElementMeshInterpolation[{mesh}, Sqrt[( sxy^2 + syy ^2 + sxx^2 )]]
]

fac = Y/(1 - \[Nu]^2)*{{1, \[Nu], 0}, {\[Nu], 1, 0}, {0, 0, (1 - \[Nu])/2}};
vonMisesStress = VonMisesStress[{uif, vif}, fac];


The factor fac changes depending on weather you are looking at plane stress or plain strain. If I recall correctly this one is for plane stress (but check! See in the comments below. According to Wiki this is to be used Sqrt[sxx^2 - sxx*syy + syy^2 + 3*sxy^2])

Edit

If you want interpolating functions for sxx, syy and sxy un-comment the part in the VonMisesStress function and comment the next line out.

End Edit

We can then visualize the stress in the anulus:

ElementMeshContourPlot[ vonMisesStress["ValuesOnGrid"],
ElementMeshDeformation[uif["ElementMesh"], {uif, vif} ],
AspectRatio -> Automatic] If you do go through the trouble verifying this I'd appreciate if you could send me a small note. One other thing that may be interesting is to use the approach you used and compare that to this approach.

• typo : weather instead of wether. Could be subject to sarcasms :) Dec 16, 2016 at 13:20
• @andre thanks!! Dec 16, 2016 at 13:25
• Thank you for the accurate answer. Anyway what if I want to use sxx, syy and sxy in another function? How can I extrapolate these values outside of VonMisesStress function?
– user44593
Jan 12, 2017 at 11:00
• @GiovannaDike, I am not sure I understand the question, do you want to return 3 interpolating functions for {sxx, syy, sxy}? What do you mean by "extrapolate [...] outside of [....][a] function"? Jan 12, 2017 at 11:25
• I'm sorry, that's my fault. By extrapolate I mean that I want to use those values outside the function VonMisesStress. For example I want to calculate the difference between the principal stresses sigma1 and sigma2 but I need the values of sxx, syy and sxy. How can I do it? I hope it's a little bit more clear now.
– user44593
Jan 12, 2017 at 16:45