# FEM stress analysis by considering bending moment

I found an example of using FEM for stress analysis in the forum, the code provided by @user21, which is very friendly for beginners.

Needs["NDSolveFEM"];

(*Physical Parameters*)
L = 1;
h = 0.1;
ss = 10;(*Shear stress on beam. This stress is applied to the boundary, and its direction seems to be upward rather than perpendicular to the centerline of the beam.*)

(*Finite element meshing*)
reg = Rectangle[{0, -h}, {L, h}];
mesh = ToElementMesh[reg];
materialParameters = {Y -> 10^3, ν -> 33/100};

(*Elasticity assumptions*)
ps = {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))}, {-((Y (1 - ν))/(2 \
(1 - ν^2))), 0}}.Inactive[Grad][v[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}]};

(*Boundary condition And mechanical equilibrium equation calculation*)
{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} ∈ mesh];

(*Force-displacement visualization*)
dmesh = ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 1];
Show[{mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]


It calculated the displacement of the beam. Based on the linear elastic assumption, the beam should be able to undergo large deformation.

So how can I change the boundary conditions applied at the right endpoint so that the direction of force is perpendicular to the centerline and achieve a bending moment effect as shown in the figure below? I think modifying boundary condition ps == {0, NeumannValue[ss, x == L]} can achieve this, but after trying it out, it seems that we cannot adjust the direction of applied force here. Does this involve coordinate transformation? The direction of the applied moment torque will change with deformation(In a fixed coordinate system). But the stress in this code is upward along the y-axis. If the elastic body is compressible, it will be stretched very long but cannot show the effect of bending.

Here is another example of adding distributed force on a certain interval, which is somewhat similar to my problem. Modifying the boundary conditions {planeStrainOperator[10^3, 33/100] == {0, NeumannValue[-5., y == 1 && 2 <= x <= 3]} in the example can solve this kind of problem.

I checked the official documentation, but did not find how to incorporate this special stress requirement into finite element calculations. The purpose of studying this example is to consider dynamic problems in the future, that is, analyzing its motion under external forces. This requires solving partial differential equations with respect to time. Dynamic problems are more complex than stress analysis, but the FEM calculation process is similar.

• @Nasser, thank you for your help in editing to increase the readability of the code. I forgot to use the M.SE editor plug-in, which led to a complex display in this code.
• Please, add Needs["NDSolveFEM"] in the beginning of the code. Also note, this is not beam you try to describe, this is plate. The beam deformation model discussed for example on mathematica.stackexchange.com/questions/250978/… Commented Jun 13, 2023 at 3:00
• @AlexTrounev Beams, plates, and shells are always confusing to me. Do you have any simple explanation to distinguish these concepts? Because the original post I read described Euler beams, so I took this example as a beam(I think the large deformation is based on the linear elastic assumption, it is OK to modify the ps equation in code, so I don't pay much attention to this aspect.). I'll go read the post you gave me and see if I misunderstood it. Thank you for your comment!