Stress analysis using finite element

Question

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["NDSolve`FEM`"];

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))), 
     0}, {0, -(Y/(1 - \[Nu]^2))}}.Inactive[Grad][
    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];

Answer 1

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

Needs["NDSolve`FEM`"];

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))), 
         0}, {0, -(Y/(1 - \[Nu]^2))}}.Inactive[Grad][
        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"][[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,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.