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["NDSolve`FEM`"];
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))),
0}, {0, -(Y/(1 - ν^2))}}.Inactive[Grad][
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
