Calculation of stress in two dimensions is covered in this post by the helpful user21. I am now trying to work in three dimensions. For this we need a stress operator which I have taken from Help.
ClearAll[stressOperator, u, v, w, x, y, z, Y, ν];
stressOperator[
Y_, ν_] := {Inactive[
Div][{{0, 0, -((Y*ν)/((1 - 2*ν)*(1 + ν)))}, {0, 0,
0}, {-Y/(2*(1 + ν)), 0, 0}}.Inactive[Grad][
w[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{0, -((Y*ν)/((1 - 2*ν)*(1 + ν))),
0}, {-Y/(2*(1 + ν)), 0, 0}, {0, 0, 0}}.Inactive[Grad][
v[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{-((Y*(1 - ν))/((1 - 2*ν)*(1 + ν))), 0,
0}, {0, -Y/(2*(1 + ν)), 0}, {0,
0, -Y/(2*(1 + ν))}}.Inactive[Grad][
u[x, y, z], {x, y, z}], {x, y, z}],
Inactive[Div][{{0, 0, 0}, {0,
0, -((Y*ν)/((1 -
2*ν)*(1 + ν)))}, {0, -Y/(2*(1 + ν)),
0}}.Inactive[Grad][w[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{0, -Y/(2*(1 + ν)),
0}, {-((Y*ν)/((1 - 2*ν)*(1 + ν))), 0, 0}, {0, 0,
0}}.Inactive[Grad][u[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{-Y/(2*(1 + ν)), 0,
0}, {0, -((Y*(1 - ν))/((1 - 2*ν)*(1 + ν))), 0}, {0,
0, -Y/(2*(1 + ν))}}.Inactive[Grad][
v[x, y, z], {x, y, z}], {x, y, z}],
Inactive[Div][{{0, 0, 0}, {0,
0, -Y/(2*(1 + ν))}, {0, -((Y*ν)/((1 -
2*ν)*(1 + ν))), 0}}.Inactive[Grad][
v[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{0, 0, -Y/(2*(1 + ν))}, {0, 0,
0}, {-((Y*ν)/((1 - 2*ν)*(1 + ν))), 0, 0}}.Inactive[
Grad][u[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{-Y/(2*(1 + ν)), 0, 0}, {0, -Y/(2*(1 + ν)), 0}, {0,
0, -((Y*(1 - ν))/((1 - 2*ν)*(1 + ν)))}}.Inactive[
Grad][w[x, y, z], {x, y, z}], {x, y, z}]}
Here is an example of the application.
Needs["NDSolve`FEM`"];
Len = 1; (*length *)
ht = 0.125; (* height *)
wd = 0.5; (* width *)
materialParameters = {Y -> 10^3, ν -> 33/100};
ss = 1;(*Shear stress on beam end*)
reg = Cuboid[{0, -wd/2, 0}, {Len, wd/2, ht}];
mesh = ToElementMesh[reg];
mesh["Wireframe"]
The problem is a thick beam clamped at one end with a shear stress on the other. To solve for the deflections I use
{uif, vif, wif} = NDSolveValue[{
stressOperator[Y, ν] == {0, 0, NeumannValue[ss, x == Len]},
DirichletCondition[u[x, y, z] == 0, x == 0],
DirichletCondition[v[x, y, z] == 0, x == 0],
DirichletCondition[w[x, y, z] == 0, x == 0]
} /. materialParameters, {u, v, w}, {x, y, z} ∈ mesh];
dmesh = ElementMeshDeformation[mesh, {uif, vif, wif},
"ScalingFactor" -> 1];
Show[{mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]
All this is working well. Now I would like to find stresses and vonMises stress. I began to look at the equations but thought I would be lazy and see if anyone has done this already. How do I get from the displacements to the stress? Thanks for your efforts if you have done this.
Edit - a module for stresses
I examined the answer from Mauricio Lobos (thank you) but could not put together a simple module for calculating stress. I did manages to put one together based on extending the 2D example to 3D. Here it is
ClearAll[stress3D];
stress3D::usage =
"stress3D[{u,v,w},{Y,ν}] determines interpolation functions for \
stress from the interpolation functions {u,v,w} for displacement in \
each direction. {Y,ν} are the modulus of elasticity and Poission \
ratio.
The output is a list of the three normal stress components, three \
shear stresses and finally the von Mises stress";
stress3D[{uif_InterpolatingFunction, vif_InterpolatingFunction,
wif_InterpolatingFunction}, {Y_, ν_}] :=
Block[{dd, df, mesh, coords, dv, ux, uy, uz, vx, vy, vz, wx, wy, wz,
gxy, gxz, gyz, c1, c2, c3, sxx, syy, szz, sxy, syz, szx, vonM},
dd = Outer[(D[#1[x, y, z], #2]) &, {uif, vif, wif}, {x, y, z}];(*
make derivative interpolation functions *)
df = Table[
Function[{x, y, z}, Evaluate[dd[[i, j]]]], {i, 3}, {j, 3}]; (*
make functions for each derivative *)
mesh = uif["Coordinates"][[1]]; (* extract coordinates from mesh *)
coords = mesh["Coordinates"];
dv = Table[df[[i, j]] @@@ coords, {i, 3}, {j, 3}]; (*
get values of derivatives at each node *)
ux = dv[[1, 1]]; (* extract individual strain components *)
uy = dv[[1, 2]];
uz = dv[[1, 3]];
vx = dv[[2, 1]];
vy = dv[[2, 2]];
vz = dv[[2, 3]];
wx = dv[[3, 1]];
wy = dv[[3, 2]];
wz = dv[[3, 3]];
gxy = (uy + vx); (* shear strains *)
gxz = (uz + wx);
gyz = (vz + wy);
c1 = Y/((1 - 2 ν) (1 + ν));(* Constitutive constants *)
c2 = (1 - ν);
c3 = Y/(2 (1 + ν));
sxx = c1 ( c2 ux + ν vy + ν wz); (*stresses *)
syy = c1 ( ν ux + c2 vy + ν wz);
szz = c1 ( ν ux + ν vy + c2 wz);
sxy = c3 gxy; (* shear stresses *)
syz = c3 gyz;
szx = c3 gxz;
vonM = Sqrt[
1/2 ((sxx - syy)^2 + (syy - szz)^2 + (szz - sxx)^2 +
6 (sxy^2 + syz^2 + szx^2))];
{
ElementMeshInterpolation[{mesh}, sxx],
ElementMeshInterpolation[{mesh}, syy],
ElementMeshInterpolation[{mesh}, szz],
ElementMeshInterpolation[{mesh}, sxy],
ElementMeshInterpolation[{mesh}, syz],
ElementMeshInterpolation[{mesh}, szx],
ElementMeshInterpolation[{mesh}, vonM]
}
]
This produces the normal stresses and the shear stress (6 interpolating functions) and also the von Mises stress (last interpolating function).
Continuing with the example the stresses are calculated by
{sxx, syy, szz, sxy, syz, szx, vM} =
stress3D[{uif, vif, wif}, {Y, ν} /. materialParameters];
The results can be checked by using the approximate engineering theory for the a beam. The force at the end of the beam and the second moment of area are given by
force = ss ht wd;
I2 = 1/12 wd ht^3;
The shear stress at one section can be calculated as
With[{x = 0.5},
Plot[{sxx[x, 0, z], (ht/2 - z) force (Len - x)/I2}, {z, 0, ht}]
]
Which agrees well. The shear stress along the beam can be calculated from
Plot[{sxx[x, 0, ht], -(ht /(2 I2)) force (Len - x),}, {x, 0, Len},
PlotLegends -> LineLegend[{"FE", "Approx."}]]
Again this agrees well except at the root where the FE is more correct. The shear stress in the z direction on the cross-section is given by
With[{x = 0.5},
Plot[{szx[x, 0, z], (ss wd ht)/(2 I2) (ht^2/4 - (z - ht/2)^2)}, {z,
0, ht}, PlotLegends -> LineLegend[{"FE", "Approx."}]]
]
Which is slightly off. This may be because the approximate theory is only good for cross-sections which are small compared to the length. Also, we could do with more grid. Finally looking at the shear stress where it is applied to the end of the beam gives.
Plot[{szx[Len, 0, z], ss}, {z, 0, ht}, PlotRange -> All,
PlotLegends -> LineLegend[{"FE", "Approx."}]]
Again this is a bit disappointing since the shear stress is applied equally over the end. Again this may be due to a poor mesh.
Overall the results look reasonable. I hope you find the module for calculating stress useful. Let me know if there are problems.
