Your boundary conditions seem to be not quite correct according to the mechanical problem. Sorry, I don't have the time to go through your code today, but I got a version running, although this will take some time and might be an overkill, since it is based on the full 3D theory. I have to go home now, I will try to take a look at your code again tomorrow, if nobody else finds the error.
I treated the problem as follows (length in x direction L1, in virtual y direction L2 and in z direction L3)
First, let's get a reference solution of the 1D theory:
(*Geometry - in m*)
L1 = 2;
L2 = 0.1;
L3 = 0.2;
Iy = L2*L3^3/12;
(*Force and densities - in N*)
F = 10;
qA = F/(L1*L2); (*area density - for 3D FEM*)
ql = qA*L2; (*line density - for 1D theory*)
(*Material parameters*)
Em = 2.1*10^9; (*Young's modulus*)
nu = 0.3;(*Poisson's ration*)
(*1D theory*)
wsol1D = DSolveValue[{
Em*Iy*D[w[x], {x, 4}] == ql
, (w[0]) == 0, (w[L1]) == 0
, (w''[0]) == 0, (w''[L1]) == 0
}, w, x];
My = -Em*Iy*wsol1D''[x];(*Moment*)
sig = My/Iy*z;(*normal stress*)
GraphicsRow[{
Plot[wsol1D[x], {x, 0, L1}, AxesLabel -> {"x", "w(x)"}]
, Plot[sig /. x -> L1/2, {z, -L3/2, L3/2},
AxesLabel -> {"z", "\[Sigma](x=L1/2,z)"}]
}
, ImageSize -> Large
]
Now, let's get the full 3D FEM solution (takes 1.4 seconds for me) with a area force density
(*FEM solution*)
Needs["NDSolve`FEM`"]
(******************************)
(*Region definition*)
reg = Cuboid[{0, -L2/2, -L3/2}, {L1, L2/2, L3/2}];
(******************************)
(*Isotropic material stiffness - fourth-order tensor*)
(*Identities*)
I2 = IdentityMatrix@3;
IdI = TensorProduct[I2, I2];
I4 = TensorTranspose[IdI, {1, 3, 2, 4}];
IS = (I4 + TensorTranspose[I4, {1, 2, 4, 3}])/2;
(*Isotropic projectors*)
P1 = 1/3*IdI;
P2 = IS - P1;
(*Isotropic stiffness*)
Ciso = l1*P1 + l2*P2;
l1 = 3*Km;
l2 = 2*Gm;
Km = 1/3*Em/(1 - 2*nu);
Gm = 1/2*Em/(1 + nu);
(******************************)
(*Equations*)
eq = Table[
Inactive[Div][
Ciso[[i, ;; , 1, ;;]].Inactive[Grad][u[x, y, z], {x, y, z}], {x,
y, z}]
+ Inactive[Div][
Ciso[[i, ;; , 2, ;;]].Inactive[Grad][v[x, y, z], {x, y, z}], {x,
y, z}]
+ Inactive[Div][
Ciso[[i, ;; , 3, ;;]].Inactive[Grad][w[x, y, z], {x, y, z}], {x,
y, z}]
, {i, 3}
];
(******************************)
(*BCs*)
(*Dirichlet*)
bcD = {
DirichletCondition[{u[x, y, z] == 0, v[x, y, z] == 0,
w[x, y, z] == 0}, x == 0 && z == 0]
, DirichletCondition[{v[x, y, z] == 0, w[x, y, z] == 0},
x == L1 && z == 0]
};
(*Neumann*)
bcN = {0, 0, NeumannValue[-qA, z == -L3/2]};
(******************************)
(*Solution*)
{usol, vsol, wsol} =
NDSolveValue[{eq == bcN, bcD}, {u, v, w}, Element[{x, y, z}, reg],
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.0001,
"MeshOrder" -> 2}}}
]; // AbsoluteTiming
{1.4435, Null}
You can take a look at the deformed mesh if you want
mesh = usol["ElementMesh"];
Show[{
mesh["Wireframe"]
, ElementMeshDeformation[mesh, {usol, vsol, wsol},
"ScalingFactor" -> 10^4][
"Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]
}, Axes -> True, AxesLabel -> {x, y, z}]
The FEM solution (FEM solution as red points) is in good accordance with the analytical 1D theory
Plot[wsol1D[x], {x, 0, L1},
Epilog -> {PointSize -> Medium, Red,
Point[Table[{x, wsol[x, 0, 0]}, {x, 0, L1, L1/10}]]}]
You can get the stress distribution at any point with the full 3D Hooke's law $\sigma_{ij} = C_{ijkl} u_{k,l}$ (FEM solution as red points) (remark: you dont need to symmetrize the displacement grandient in my code in order to obtain the infinitesimal strain, since the stiffness $C_{ijkl}$ I used symmetrizes automatically the mapped tensor). Let's get $\sigma_{xx} = \sigma_{11}$
uv[x_, y_, z_] := {usol[x, y, z], vsol[x, y, z], wsol[x, y, z]}
eps[xs_, ys_, zs_] :=
D[uv[x, y, z], {{x, y, z}, 1}] /. {x -> xs, y -> ys, z -> zs}
(*linear map of second order tensor B over fourth-order tensor A*)
lm[A_, B_] := TensorContract[TensorProduct[A, B], {{3, 5}, {4, 6}}]
(*Get Cauchy stress sigma_xx = sigma[[1,1]], at x=L1/2 depending on z with 3D Hooke's law*)
sigloc = lm[Ciso, eps[L1/2, 0, z]][[1, 1]];
siglocdata = Table[{zi, sigloc /. z -> zi}, {zi, -L3/2, L3/2, L3/10}];
Plot[sig /. x -> L1/2, {z, -L3/2, L3/2},
AxesLabel -> {"z", "\[Sigma](x=L1/2,z)"},
Epilog -> {PointSize -> Medium, Red, Point@siglocdata}]
