Finite ELement Model - simple beam, ends pinned, uniform load over 1/2 of beam, getting several errors

Question

here is the code I am using - having issues getting started in FEA for this simple beam model. The problem is simple - I have looked at posted examples of MMA FEA, but they are all 2D and 3D problems. If I can get this simple 1D model made and play with it, I know I can get a better understanding of the general approach for MMA FEA.

(*Define parameters*)
Ls = 10; (*Length of the beam in meters*)
Es = 210*10^9; (*Young's modulus in Pascals*)
Is = 1*10^-6; (*Moment of inertia in m^4*)
qs = 1000; (*Uniform load in N/m applied over half the beam*)

(*Load function definition*)
loadFunction[x_] := If[x <= Ls/2, qs, 0];

(*Beam equation setup using FEM*)
Needs["NDSolve`FEM`"];

(*ToElementMesh["Coordinates"\[Rule]{{0},{Ls}},"MeshElements"\[Rule]{\
LineElement[{{1,2}}]}]
beamDomain=ToElementMesh[LineElement[{0,0},{0,Ls}]];*)
beamDomain = ToElementMesh[Line[{0, 0}, {0, Ls}]];

(*Check if mesh is created properly*)
MeshRegion[beamDomain]

(*Boundary conditions for a simply supported beam*)
boundaryConditions = {w[0] == 0, w[Ls] == 0, Derivative[2][w][0] == 0,
    Derivative[2][w][Ls] == 0};

(*Beam equation*)
beamPDE = Es*Is*D[w[x], {x, 4}] + loadFunction[x] == 0;

(*Solve the beam equation using NDSolve with FEM*)
solution = 
  NDSolveValue[{beamPDE, boundaryConditions}, w, {x, 0, Ls}, 
   Method -> {"FiniteElement"}];

(*Plot the deflection diagram*)
Plot[solution[x], {x, 0, Ls}, PlotLabel -> "Deflection Diagram", 
 AxesLabel -> {"x (m)", "w(x) (m)"}]

(*Calculate shear force V(x)=-EI*d^3w/dx^3*)
shearForce = -Es*Is*D[solution[x], {x, 3}];
(*Plot the shear force diagram*)
Plot[shearForce, {x, 0, Ls}, PlotLabel -> "Shear Force Diagram", 
 AxesLabel -> {"x (m)", "V(x) (N)"}]

(*Calculate bending moment M(x)=-EI*d^2w/dx^2*)
bendingMoment = -Es*Is*D[solution[x], {x, 2}];
(*Plot the bending moment diagram*)
Plot[bendingMoment, {x, 0, Ls}, PlotLabel -> "Bending Moment Diagram",
  AxesLabel -> {"x (m)", "M(x) (N*m)"}]

Answer 1

To make the FEM-solution work you have to change the element mesh and you have to transform the pde-system to second order.

From mechanical point of view you only have to define an additional unknown moment[x] == -Es*Is*D[w[x], {x, 2}]

(*Define parameters*)Ls = 10;(*Length of the beam in meters*)Es = 
 210*10^9;(*Young's modulus in Pascals*)Is = 
 1*10^-6;(*Moment of inertia in m^4*)qs = 1000;(*Uniform load in N/m \
applied over half the beam*)(*Load function definition*)
loadFunction[x_] := If[x <= Ls/2, qs, 0];

(*Beam equation setup using FEM*)
Needs["NDSolve`FEM`"];
 
beamDomain = ToElementMesh [ImplicitRegion[0 <= x <= Ls, x]] ;

(*Boundary conditions for a simply supported beam*)
boundaryConditions = {w[0] == 0, w[Ls] == 0, moment[0] == 0, 
   moment[Ls] == 0};

(*Beam equation*)
beamPDE = { D[moment[x], {x, 2}] + loadFunction[x] == 0, 
   moment[x] == -Es*Is*D[w[x], {x, 2}]};

(*Solve the beam equation using NDSolve with FEM*)
solution = 
  NDSolveValue[{beamPDE, boundaryConditions}, {w, moment}, {x, 0, Ls},
    Method -> {"FiniteElement"}] ;
(*Plot the deflection diagram*)
Plot[solution[[1]][x], {x, 0, Ls}, PlotLabel -> "Deflection Diagram", 
 AxesLabel -> {"x (m)", "w(x) (m)"}]

(*Calculate shear force V(x)=-EI*d^3w/dx^3*)
shearForce = -Es*Is*D[solution[x], {x, 3}];
(*Plot the shear force diagram*)
Plot[solution[[2]]'[x], {x, 0, Ls}, 
 PlotLabel -> "Shear Force Diagram", 
 AxesLabel -> {"x (m)", "V(x) (N)"}]

  
(*Plot the bending moment diagram*)
Plot[solution[[2]][x], {x, 0, Ls}, 
 PlotLabel -> "Bending Moment Diagram", 
 AxesLabel -> {"x (m)", "M(x) (N*m)"}]

Hope it helps!