# Stress calculations using finite elements

A standard engineering problem is to calculate stresses in a structure due to applied forces. With the inclusion of the finite element method in version 10 this question attempts to investigate how this may be done. Normal, Subscript[σ, x] Subscript[σ, y] and shear Subscript[σ, x y] stresses are calculate from displacements (in two dimensions) by These stresses are gradients of the displacements u(x,y) and v(x,y) in the x and y directions. When these equations are combined with the equilibrium conditions we get the differential equations for plain stress Unfortunately these equations cannot be entered directly into NDSolveValue. I tried in this post and the ever helpful user21 showed that the way the equation is entered makes a significant difference. I feel I should be able to construct the required equation from mine but it is beyond me. A secondary question is can someone write a parser that could interpret textbook equations in the needed manner? The question here is how to best extract stresses from a finite element analysis. Also, how can one conveniently put in stresses (or forces) on the boundaries? The required version of the differential equation has been provided by user21. My first attempt is to input the above equation for stress together with the required version of the differential equation.

Needs["NDSolveFEM"];
ps = {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}]};

L = 1;
h = 0.125;
ss = 5; (* Shear stress on beam *)
reg = Rectangle[{0, -h}, {L, h}];
mesh = ToElementMesh[reg];
mesh["Wireframe"] {uif, vif, σxif, σyif, σxyif} = NDSolveValue[{
ps == {0, NeumannValue[ss, x == L]},

σx[x, y] ==
Y/(1 - ν^2) (D[u[x, y], x] + ν D[v[x, y], y] ),
σy[x, y] ==
Y/(1 - ν^2) (D[v[x, y], y] + ν D[u[x, y], x] ),
σxy[x, y] ==
Y/(1 - ν^2) (1 - ν)/2 (D[u[x, y], y] + D[v[x, y], x] ),

DirichletCondition[u[x, y] == 0, x == 0],
DirichletCondition[v[x, y] == 0, x == 0]

} /. {Y -> 10^3, ν -> 33/100},
{u, v, σx, σy, σxy},
{x, y} ∈ mesh];


The displacement results are

dmesh = ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 1];
Show[{
mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}] The stress results look good

Plot3D[σxif[x, y], {x, y} ∈ mesh,
BoxRatios -> {4, 1, 1}, PlotRange -> All]
Plot3D[σxyif[x, y], {x, y} ∈ mesh,
BoxRatios -> {4, 1, 1}, PlotRange -> All] As a check one can use the approximate results of Euler-Bernoulli beam theory which is good away from the ends where the finite element method is more correct.

Plot[{σxif[L/2, y], -((ss 12 2 h  0.5)/(2 h)^3) y }, {y, -h,
h}, PlotLegends -> LineLegend[{"Calculated", "Theory"}]]
Plot[{σxyif[L/2, y], 6  ss/(8 h^3) 2 h (h^2 - y^2)}, {y, -h,
h}, PlotLegends -> LineLegend[{"Calculated", "Theory"}]] These results are very good.

As these equations work I was hoping that the stresses on the boundaries could be entered directly as a DirichletCondition rather than NeumannValue

{uif, vif, σxif, σyif, σxyif} = NDSolveValue[{
ps == {0, 0},

σx[x, y] ==
Y/(1 - ν^2) (D[u[x, y], x] + ν D[v[x, y], y] ),
σy[x, y] ==
Y/(1 - ν^2) (D[v[x, y], y] + ν D[u[x, y], x] ),
σxy[x, y] ==
Y/(1 - ν^2) (1 - ν)/2 (D[u[x, y], y] + D[v[x, y], x] ),

DirichletCondition[u[x, y] == 0, x == 0],
DirichletCondition[v[x, y] == 0, x == 0],
DirichletCondition[σxy[x, y] == ss, x == L],
DirichletCondition[σx[x, y] == 0, x == L],
DirichletCondition[σy[x, y] == 0, x == L]

} /. {Y -> 10^3, ν -> 33/100},
{u, v, σx, σy, σxy},
{x, y} ∈ mesh];

dmesh = ElementMeshDeformation[mesh, {uif, vif},
"ScalingFactor" -> 10];
Show[{
mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}] Plot3D[σxyif[x, y], {x, y} ∈ mesh,
BoxRatios -> {4, 1, 1}, PlotRange -> All] However, this approach does not work. So what is the best way of putting stresses into, and getting stresses out of, the finite element method? Thanks

• As far as I remember, that is not possible (to define stresses as Dirichlet conditions), since in FEM the stress condition are always classified as Neumann condition, due to their definition in terms of the gradient of the unknown function. – Mauricio Fernández Dec 21 '15 at 8:46
• @MauricioLobos I agree that it is a Neumann condition on displacements and A Dirichlet condition on stress. With my stress equations I feel I am putting the boundary conditions in correctly. Hence the Dirichlet conditions should feed the stress values which in turn feed the displacements. If I was doing this by hand that would be how it works. It seems that the finite element approach needs to be given the equations in a particular form and I don't seem to be able to work my equations into that form. Should I be trying to express my equations in the weak form? Can MM help me to do this? – Hugh Dec 21 '15 at 10:16
• @kuba Thanks for formatting my Greek. Is there a way I can do this myself? – Hugh Dec 21 '15 at 10:18
• @Hugh I'm using: meta.mathematica.stackexchange.com/q/1043 – Kuba Dec 21 '15 at 10:19
• @Hugh I dont think you have to formulate your equation in the weak form, the Mathematica FEM does that for you automatically. But only for equations defined as in the documentation. If you want to treat the stresses as an independent function, may be you can just enter the differential equation for the stresses without any connection to the displacements (only div(sigma)=0, sigma being the unknown function). But I think the Mathematica FEM is not designed for that kind of equation structure (see documentation) and you dont know the Dirichlet conditions of sigma at x=0. – Mauricio Fernández Dec 21 '15 at 11:43

There are several part to you question and I'll try to tackle them one by one (not sure I'll be able to do this in one go)

How to calculate stresses

Here is the model problem.

L = 1;
h = 0.125;
(*Shear stress on beam*)
ss = 5;
reg = Rectangle[{0, -h}, {L, h}];
mesh = ToElementMesh[reg];
materialParameters = {Y -> 10^3, \[Nu] -> 33/100};


We use a plane stress formulation as you derived it and under consideration that we do not want the operator to evaluate:

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


We then solve

{uif, vif} =
NDSolveValue[{ps == {0, NeumannValue[ss, x == L]},
DirichletCondition[u[x, y] == 0, x == 0],
DirichletCondition[v[x, y] == 0, x == 0]} /.
materialParameters, {u, v}, {x, y} \[Element] mesh];


Note that I only included the plane stress formulation. We will recover the derivatives for the stress computation just now.

dmesh = ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 1];
Show[{mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}] Here is a function to allow to recover the derivatives:

ClearAll[VonMisesStress]
VonMisesStress[{uif_InterpolatingFunction,
vif_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]) &, {uif, vif}, {x, y}];
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},
Sqrt[(sxy^2) + (syy^2)+(sxx^2)]]*)
{sxx, syy, sxy}
]


If you uncomment the ElementMeshInterpolation then this will compute the vonMises stress, but now it returns the stresses.

fac = Y/(1 - \[Nu]^2)*{{1, \[Nu], 0}, {\[Nu], 1, 0}, {0,
0, (1 - \[Nu])/2}};
facm = fac /. materialParameters;
fac // Simplify
{{Y/(1 - \[Nu]^2), (Y \[Nu])/(1 - \[Nu]^2), 0}, {(Y \[Nu])/(
1 - \[Nu]^2), Y/(1 - \[Nu]^2), 0}, {0, 0, Y/(2 + 2 \[Nu])}}


The factor is for plane stress; it needs to be adjusted for plain strain accordingly.

(*vonMisesStress =VonMisesStress[{uif,vif},facm]*)
{sxx, syy, sxy} =
VonMisesStress[{uif, vif}, facm];
ifsxx = ElementMeshInterpolation[{mesh}, sxx];
ifsyy = ElementMeshInterpolation[{mesh}, syy];
ifsxy = ElementMeshInterpolation[{mesh}, sxy];


Let's compare this with your approach:

(* Plot3D[ifsxx[x, y] - \[Sigma]xif[x, y], {x, y} \[Element] mesh,
BoxRatios -> {4, 1, 1}, PlotRange -> All] *)
Plot3D[ifsxy[x, y] - \[Sigma]xyif[x, y], {x, y} \[Element] mesh,
BoxRatios -> {4, 1, 1}, PlotRange -> All] The differences at the boundary are expected, I think. Additionally, this has the advantage of not creating additional equations that need to be solved for.

Visualize the vonMises stress in the deformed beam:

ElementMeshContourPlot[Sqrt[(sxx^2) + (syy^2) + (sxy^2)],
ElementMeshDeformation[uif["ElementMesh"], {uif, vif} ],
AspectRatio -> Automatic] So, this is how I'd compute the stresses.

Now, let's consider some of your other questions.

How to specify forces/stresses on the region and boundary.

A force on the entire body (e.g. gravity)

{uif, vif} =
NDSolveValue[{ps == {0, -9.8},
DirichletCondition[u[x, y] == 0, x == 0],
DirichletCondition[v[x, y] == 0, x == 0]} /.
materialParameters, {u, v}, {x, y} \[Element] mesh];
dmesh = ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 1];
Show[{mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}] A force on the boundary:

{uif, vif} =
NDSolveValue[{ps == {0, NeumannValue[ss, x == L]},
DirichletCondition[u[x, y] == 0, x == 0],
DirichletCondition[v[x, y] == 0, x == 0]} /.
materialParameters, {u, v}, {x, y} \[Element] mesh];
dmesh = ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 1];
Show[{mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}] Combined body and boundary force: (Gravity and a load on the right)

{uif, vif} =
NDSolveValue[{ps == {0, -9.8 + NeumannValue[ss, x == L]},
DirichletCondition[u[x, y] == 0, x == 0],
DirichletCondition[v[x, y] == 0, x == 0]} /.
materialParameters, {u, v}, {x, y} \[Element] mesh];
dmesh = ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 1];
Show[{mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}] Sometimes one would like to model geometries that are under an initial stress (e.g. a pre-stressed beam that will be un-stressed once a load is applied). To pre-stress one can use:

preStress = {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][
Inactive[
Plus][({{-(Y/(1 - \[Nu]^2)),
0}, {0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}}.Inactive[
Grad][u[x, y], {x, y}]), {{0}, {0}}], {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][
Inactive[
Plus][({{-((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2))),
v[x, y], {x, y}]), {{ss}, {0}}], {x, y}]};


Note that now we have added a + {{ss},{0}} to the second equation (for completeness I added a {{0},{0}} to the first one as well)

Inspecting what is parsed one gets:

{state} =
NDSolveProcessEquations[{preStress == {0, 0},
DirichletCondition[u[x, y] == 0, x == 0],
DirichletCondition[v[x, y] == 0, x == 0]} /.
materialParameters, {u, v}, {x, y} \[Element] mesh];
state["FiniteElementData"][
{{{{0}, {0}}}, {{{5}, {0}}}}


And solving the equation:

{uif, vif} =
NDSolveValue[{preStress == {0, 0},
DirichletCondition[u[x, y] == 0, x == 0],
DirichletCondition[v[x, y] == 0, x == 0]} /.
materialParameters, {u, v}, {x, y} \[Element] mesh];
dmesh = ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 1];
Show[{mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}] Which is the same as your example. Note that no boundary conditions on the right hand side were specified. This beam is pre-stressed.

I do not know how (if at all) one can achieve this with adding the stress equations as you did. I am not saying it's not possible but I don't know how.

A parser:

You could use something like this:

ClearAll[PlaneStress];
PlaneStress[{Y_, nu_}, {u_, v_}, X : {x_, y_}] :=
Module[{pStress},
pStress = -Y/(1 -
nu^2)*{{{{1, 0}, {0, (1 - nu)/2}}, {{0, nu}, {(1 - nu)/2,
0}}}, {{{0, (1 - nu)/2}, {nu, 0}}, {{(1 - nu)/2, 0}, {0, 1}}}};

]
PlaneStress[{Y, nu}, {u[x, y], v[x, y]}, {x, y}]

{Inactive[Div][{{0, -((nu*Y)/(1 - nu^2))}, {-((1 - nu)*Y)/(2*(1 - nu^2)), 0}} .
Inactive[Grad][v[x, y], {x, y}], {x, y}] +
Inactive[Div][{{-(Y/(1 - nu^2)), 0}, {0, -((1 - nu)*Y)/(2*(1 - nu^2))}} .
Inactive[Grad][u[x, y], {x, y}], {x, y}],
Inactive[Div][{{0, -((1 - nu)*Y)/(2*(1 - nu^2))}, {-((nu*Y)/(1 - nu^2)), 0}} .
Inactive[Grad][u[x, y], {x, y}], {x, y}] +
Inactive[Div][{{-((1 - nu)*Y)/(2*(1 - nu^2)), 0}, {0, -(Y/(1 - nu^2))}} .
Inactive[Grad][v[x, y], {x, y}], {x, y}]}


Hope this helps.

• is there also a multi parallel solver version in forum? – ABCDEMMM Jul 29 '19 at 22:36
• @ABCDEMMM, I am not 100% I understand your question but there is a Domain Decomposition solver in the FEMAddOns. you can have a look at the documentation here – user21 Jul 30 '19 at 6:27