Finite element boundary breaking

Question

I am trying to model a cantilever beam which can vibrate. On the left the beam is clamped. I am largely following user21 here and also the example in help. I start by doing a static beam which works then I work out the vibration modes (eigenvectors) of a beam with no boundary constraint, this works, and then it all goes wrong when I try to do the eigenvectors of the clamped beam. I start with some helper modules.

Needs["NDSolve`FEM`"];
(*plane stress equations two versions with and without time *)
ClearAll[planeStress];
planeStress::usage = 
  "planeStress[u,v,x,y,Y,\[Nu]] (no time) or planeStress[u,v,t,x,y,Y,\
\[Nu]] (with time) Y is modulus of elasticity and \[Nu] is Poission \
ratio. ";
planeStress[u_, v_, x_, y_, 
  Y_, \[Nu]_] := {Inactive[
     Div][{{0, -((Y*\[Nu])/(1 - \[Nu]^2))}, {-(Y*(1 - \[Nu]))/(2*(1 - \
\[Nu]^2)), 0}}.Inactive[Grad][v[x, y], {x, y}], {x, y}] + 
   Inactive[
     Div][{{-(Y/(1 - \[Nu]^2)), 
       0}, {0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}}.Inactive[Grad][
      u[x, y], {x, y}], {x, y}], 
  Inactive[Div][{{0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}, \
{-((Y*\[Nu])/(1 - \[Nu]^2)), 0}}.Inactive[Grad][u[x, y], {x, y}], {x, 
     y}] + Inactive[
     Div][{{-(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2)), 
       0}, {0, -(Y/(1 - \[Nu]^2))}}.Inactive[Grad][
      v[x, y], {x, y}], {x, y}]}
planeStress[u_, v_, t_, x_, y_, 
  Y_, \[Nu]_] := {Inactive[
     Div][{{0, -((Y*\[Nu])/(1 - \[Nu]^2))}, {-(Y*(1 - \[Nu]))/(2*(1 - \
\[Nu]^2)), 0}}.Inactive[Grad][v[t, x, y], {x, y}], {x, y}] + 
   Inactive[
     Div][{{-(Y/(1 - \[Nu]^2)), 
       0}, {0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}}.Inactive[Grad][
      u[t, x, y], {x, y}], {x, y}], 
  Inactive[Div][{{0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}, \
{-((Y*\[Nu])/(1 - \[Nu]^2)), 0}}.Inactive[Grad][
      u[t, x, y], {x, y}], {x, y}] + 
   Inactive[
     Div][{{-(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2)), 
       0}, {0, -(Y/(1 - \[Nu]^2))}}.Inactive[Grad][
      v[t, x, y], {x, y}], {x, y}]}

ClearAll[finiteElementData];
finiteElementData::usate = 
  "finiteElementData[start] Extract the finite element data";
finiteElementData[state_] := 
 Module[{femdata, initBCs, methodData, initCoeffs, vd, sd, 
   discretePDE, discreteBCs, load, stiffness, damping, spans, dimDM},
  femdata = state["FiniteElementData"];
  initBCs = femdata["BoundaryConditionData"];
  methodData = femdata["FEMMethodData"];
  initCoeffs = femdata["PDECoefficientData"];
  vd = methodData["VariableData"];
  sd = NDSolve`SolutionData[{"Space" -> ToNumericalRegion[mesh]}];
  (*Discretize the PDE and the boundary conditions:*)
  discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
  discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];
  (*Extract the system matrices:*)
  load = discretePDE["LoadVector"];
  stiffness = discretePDE["StiffnessMatrix"];
  damping = discretePDE["DampingMatrix"];
  DeployBoundaryConditions[{load, stiffness, damping}, discreteBCs];
  spans = 
   Span @@@ 
    Transpose[{Most[# + 1], Rest[#]} &[methodData["IncidentOffsets"]]];
  dimDM = Dimensions[discreteBCs["DirichletMatrix"]];
  {load, stiffness, damping, methodData, spans, dimDM}
  ]


ClearAll[solveEigensystem];
solveEigensystem::usage = 
  "solveEigensystem[n,mesh,data] where n is the number of \
eigenvalues, and data is the output from finiteElementData";
solveEigensystem[n_, 
  mesh_, {load_, stiffness_, damping_, methodData_, spans_, dimDM_}] := 
 Module[{nDiri, numEigen, mm, res, eigenValues, eigenVectors, 
   eigenVectorsIF},
  nDiri = If[Length[#] > 0, First[#], 0] &[dimDM];
  numEigen = n + nDiri;
  mm = LinearSolve[damping, stiffness];
  res = Eigensystem[mm, -numEigen];
  res = Reverse /@ res;
  eigenValues = res[[1, nDiri + 1 ;; Abs[numEigen]]];
  eigenVectors = res[[2, nDiri + 1 ;; Abs[numEigen]]];
  eigenVectors = 
   Transpose[
    Developer`ToPackedArray[eigenVectors[[All, #]] & /@ spans], {2, 1,
      3}];
  eigenVectorsIF = Table[{}, {n}, {Length[spans]}];
  Do[eigenVectorsIF[[i, j]] = 
    NDSolve`FEM`ElementMeshInterpolation[{mesh}, 
     eigenVectors[[i, j]]], {i, n}, {j, Length[spans]}];
  {eigenValues, eigenVectorsIF}
  ]

Here is the beam

(* make mesh *)
\[CapitalOmega] = ImplicitRegion[True, {x, y}];
mesh = ToElementMesh[\[CapitalOmega], {{0, 5}, {0, 0.25}}, 
   "MaxCellMeasure" -> 0.005];
mesh["Wireframe"]

Mathematica graphics

Starting with the static case -all is good

(* Set up boundary conditions and PDE*)
bcs = DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, x == 0];
pde2D = planeStress[u, v, x, y, 10^3, 33/100] == {0, 
    NeumannValue[0.1, x == 5]};
{state} = 
  NDSolve`ProcessEquations[{pde2D, bcs}, {u, v}, {x, y} \[Element] 
    mesh,
   Method -> {"PDEDiscretization" -> {"FiniteElement"}}];

Extract the data, solve and plot.

{load, stiffness, damping, methodData, spans, dimDM} = 
  finiteElementData[state];
solution = LinearSolve[stiffness, load];
uif = ElementMeshInterpolation[{mesh}, solution[[spans[[1]]]]];
vif = ElementMeshInterpolation[{mesh}, solution[[spans[[2]]]]];
dmesh = ElementMeshDeformation[mesh, {uif, vif}];
Show[{
  mesh["Wireframe"],
  dmesh["Wireframe"[
    "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]

Mathematica graphics

Next for vibration of the beam with no supports

(* Set up boundary conditions and PDE *)
bcs = Sequence[];
pde2D = {D[u[t, x, y], t], D[v[t, x, y], t]} + 
    planeStress[u, v, t, x, y, 10^3, 33/100] == {0, 0};
{state} = 
  NDSolve`ProcessEquations[{pde2D, bcs, u[0, x, y] == 0, 
    v[0, x, y] == 0}, {u, v}, {t, 0, 1}, {x, y} \[Element] mesh,
   Method -> {"PDEDiscretization" -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"FiniteElement"}}}];

fed = finiteElementData[state];
{eigenValues, eigenVectorsIF} = solveEigensystem[8, mesh, fed];

The theoretical values can be calculated assuming a slender beam. The first three are zero corresponding to rigid body modes i.e. moving horizontally, vertically and rocking with no flexure. The next four are transverse flexural vibration. The last value corresponds to longitudinal vibration.

theory = {0, 0, 0, 22/L^2 Sqrt[(Y d^2)/(12 1)], 
    61.7/L^2 Sqrt[(Y d^2)/(12 1)], 121/L^2 Sqrt[(Y d^2)/(12 1)], 
    200/L^2 Sqrt[(Y d^2)/(12 1)], \[Pi]/L Sqrt[Y/1.]} /. {Y -> 10^3, 
    d -> 0.25, L -> 5};
TableForm[Transpose[{Sqrt[Abs[eigenValues]], theory}], 
 TableHeadings -> {Automatic, { "Calculated", "Theory"}}]

Mathematica graphics

The eigenvectors plot as

Column[Table[uif = eigenVectorsIF[[n, 1]];
  vif = eigenVectorsIF[[n, 2]];
  dmesh = 
   ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 10];
  Show[{
    mesh["Wireframe"],
    dmesh[
     "Wireframe"[
      "ElementMeshDirective" -> 
       Directive[EdgeForm[Red], FaceForm[]]]]}],
  {n, 8}]]

Mathematica graphics

The first three mode shapes need normalizing to give the three rigid body modes corresponding to each coordinate. Now things go wrong - I put in the boundary conditions

bcs = DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, x == 0];
pde2D = {D[u[t, x, y], t], D[v[t, x, y], t]} + 
    planeStress[u, v, t, x, y, 10^3, 33/100] == {0, 0};
{state} = 
  NDSolve`ProcessEquations[{pde2D, bcs, u[0, x, y] == 0, 
    v[0, x, y] == 0}, {u, v}, {t, 0, 1}, {x, y} \[Element] mesh,
   Method -> {"PDEDiscretization" -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"FiniteElement"}}}];
fed = finiteElementData[state];
{eigenValues, eigenVectorsIF} = solveEigensystem[5, mesh, fed];

The theoretical solutions are as follows (note that the fourth solution is for longitudinal vibration)

 theory = {3.52 Sqrt[(Y d^2)/(12 L^4)], 22 Sqrt[(Y d^2)/(12 L^4)], 
    61.7 Sqrt[(Y d^2)/(12 L^4)], \[Pi] /2 Sqrt[Y/L^2], 
    121 Sqrt[(Y d^2)/(12 L^4)]} /. {Y -> 10^3, d -> 0.25, L -> 5.};
TableForm[Transpose[{Sqrt[Abs[eigenValues]], theory}], 
 TableHeadings -> {Automatic, { "Calculated", "Theory"}}]

Mathematica graphics

Column[Table[uif = eigenVectorsIF[[n, 1]];
  vif = eigenVectorsIF[[n, 2]];
  dmesh = 
   ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 10];
  Show[{
    mesh["Wireframe"],
    dmesh[
     "Wireframe"[
      "ElementMeshDirective" -> 
       Directive[EdgeForm[Red], FaceForm[]]]]}],
  {n, 5}]]

Mathematica graphics

The first mode frequency and mode shape are wrong. The solution is badly distorted at the left hand end. If we try with a mesh that is fine at the left hand end we get

Lx = 5; Ly = 0.25;
inc = 0.01;
pts = Flatten[
   Join[{{{0, 0}, {Lx, 0}, {Lx, Ly}, {0, Ly}}, 
     Table[{0, y}, {y, Ly - inc, inc, -inc}]}], 1];
boundarymesh = ToBoundaryMesh[Polygon[pts]];
mesh = ToElementMesh[boundarymesh, "MaxCellMeasure" -> 0.005];
mesh["Wireframe"]

Mathematica graphics

(* Set up boundary conditions and PDE *)
bcs = DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, x == 0];
pde2D = {D[u[t, x, y], t], D[v[t, x, y], t]} + 
    planeStress[u, v, t, x, y, 10^3, 33/100] == {0, 0};
{state} = 
  NDSolve`ProcessEquations[{pde2D, bcs, u[0, x, y] == 0, 
    v[0, x, y] == 0}, {u, v}, {t, 0, 1}, {x, y} \[Element] mesh,
   Method -> {"PDEDiscretization" -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"FiniteElement"}}}];
fed = finiteElementData[state];
{eigenValues, eigenVectorsIF} = solveEigensystem[5, mesh, fed];
TableForm[Transpose[{Sqrt[Abs[eigenValues]], theory}], 
 TableHeadings -> {Automatic, { "Calculated", "Theory"}}]

Mathematica graphics

 Column[Table[uif = eigenVectorsIF[[n, 1]];
      vif = eigenVectorsIF[[n, 2]];
      dmesh = 
       ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 10];
      Show[{
        mesh["Wireframe"],
        dmesh[
         "Wireframe"[
          "ElementMeshDirective" -> 
           Directive[EdgeForm[Red], FaceForm[]]]]}],
      {n, 5}]]

Mathematica graphics

Again the first mode is wrong. It is clear that the left hand end elements are badly distorted. What is going wrong? Thanks

Community · Accepted Answer · 2017-04-13 12:55:40Z

I am going to show you where the code that you refer to goes south and how to fix that. However, I'd like to stress that using NDEigensystem does this correctly; and it does more that is not considered here and the other post. A second point I'd like to make is that one can compare Euler beams with what is done here, but one should be aware that the theories are different and a discrepancy in the result is expected. So, it's OK if the results do not match 100%. Let's look at NDEigensystem first.

Define a plane stress operator:

ps = {Inactive[
       Div][{{0, -((Y*\[Nu])/(1 - \[Nu]^2))}, {-(Y*(1 - \[Nu]))/(2*(1 \
- \[Nu]^2)), 0}}.Inactive[Grad][v[x, y], {x, y}], {x, y}] + 
     Inactive[
       Div][{{-(Y/(1 - \[Nu]^2)), 
         0}, {0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}}.Inactive[Grad][
        u[x, y], {x, y}], {x, y}], 
    Inactive[
       Div][{{0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}, {-((Y*\[Nu])/(1 \
- \[Nu]^2)), 0}}.Inactive[Grad][u[x, y], {x, y}], {x, y}] + 
     Inactive[
       Div][{{-(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2)), 
         0}, {0, -(Y/(1 - \[Nu]^2))}}.Inactive[Grad][
        v[x, y], {x, y}], {x, y}]} /. {Y -> 10^3, \[Nu] -> 33/100};

Your first example is then:

{vals, funs} = 
  NDEigensystem[{ps}, {u, v}, {x, y} \[Element] 
    Rectangle[{0, 0}, {5, 0.25}], 8];

This gives:

theory = {0, 0, 0, 22/L^2 Sqrt[(Y d^2)/(12 1)], 
    61.7/L^2 Sqrt[(Y d^2)/(12 1)], 121/L^2 Sqrt[(Y d^2)/(12 1)], 
    200/L^2 Sqrt[(Y d^2)/(12 1)], \[Pi]/L Sqrt[Y/1.]} /. {Y -> 10^3, 
    d -> 0.25, L -> 5};
TableForm[Transpose[{Sqrt[Abs[vals]], theory}], 
 TableHeadings -> {Automatic, {"Calculated", "Theory"}}]

The second case is:

bcs = DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, x == 0];
{vals, funs} = 
  NDEigensystem[{ps, bcs}, {u, v}, {x, y} \[Element] 
    Rectangle[{0, 0}, {5, 0.25}], 5];

Checking:

theory = {3.52 Sqrt[(Y d^2)/(12 L^4)], 22 Sqrt[(Y d^2)/(12 L^4)], 
    61.7 Sqrt[(Y d^2)/(12 L^4)], \[Pi]/2 Sqrt[Y/L^2], 
    121 Sqrt[(Y d^2)/(12 L^4)]} /. {Y -> 10^3, d -> 0.25, L -> 5.};
TableForm[Transpose[{Sqrt[Abs[vals]], theory}], 
 TableHeadings -> {Automatic, {"Calculated", "Theory"}}]

Visualizing:

Needs["NDSolve`FEM`"]
mesh = funs[[1, 1]]["ElementMesh"];
Column[Table[uif = funs[[n, 1]];
  vif = funs[[n, 2]];
  dmesh = 
   ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 0.1];
  Show[{mesh["Wireframe"], 
    dmesh["Wireframe"[
      "ElementMeshDirective" -> 
       Directive[EdgeForm[Red], FaceForm[]]]]}], {n, 5}]]

What went wrong in the code above. In essence you need to remove the DirichletCondition from the system of equations. You'd notice this if you use Eigensystem with the method Arnoldi option. This gives you a hint with an error message saying that the matrices are not symmetric. Here is the code modified to do that:

    ClearAll[finiteElementData];
    finiteElementData::usate = 
      "finiteElementData[start] Extract the finite element data";
    finiteElementData[state_] := 
     Module[{femdata, initBCs, methodData, initCoeffs, vd, sd, 
       discretePDE, discreteBCs, load, stiffness, damping, spans, dimDM},
      femdata = state["FiniteElementData"];
      initBCs = femdata["BoundaryConditionData"];
      methodData = femdata["FEMMethodData"];
      initCoeffs = femdata["PDECoefficientData"];
      vd = methodData["VariableData"];
      sd = NDSolve`SolutionData[{"Space" -> ToNumericalRegion[mesh]}];
      (*Discretize the PDE and the boundary conditions:*)

      discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
      discreteBCs = 
       DiscretizeBoundaryConditions[initBCs, methodData, sd];
      (*Extract the system matrices:*)

      load = discretePDE["LoadVector"];
      stiffness = discretePDE["StiffnessMatrix"];
      damping = discretePDE["DampingMatrix"];
      DeployBoundaryConditions[{load, stiffness, damping}, discreteBCs, 
       "ConstraintMethod" -> "Remove"];
      spans = 
       Span @@@ 
        Transpose[{Most[# + 1], Rest[#]} &[methodData["IncidentOffsets"]]];
      dimDM = Dimensions[discreteBCs["DirichletMatrix"]];
      {load, stiffness, damping, methodData, spans, dimDM, discreteBCs}
      ]

ClearAll[solveEigensystem];
solveEigensystem::usage = 
  "solveEigensystem[n,mesh,data] where n is the number of \
eigenvalues, and data is the output from finiteElementData";
solveEigensystem[n_, 
  mesh_, {load_, stiffness_, damping_, methodData_, spans_, dimDM_, 
   discreteBCs_}] := 
 Module[{nDiri, numEigen, mm, res, eigenValues, eigenVectors, 
   eigenVectorsIF},
  nDiri = If[Length[#] > 0, First[#], 0] &[dimDM];
  (*numEigen=n+nDiri;*)
  numEigen = n;
  (*mm=LinearSolve[damping,stiffness];
  res=Eigensystem[mm,-numEigen,Method->"Arnoldi"];*)

  res = Eigensystem[{stiffness, damping}, -numEigen, 
    Method -> "Arnoldi"];
  res = Reverse /@ res;
  eigenValues = res[[1, 1 ;; Abs[numEigen]]];
  eigenVectors = res[[2, 1 ;; Abs[numEigen]]];

  eigenVectors = 
   NDSolve`FEM`DirichletValueReinsertion[#, discreteBCs] & /@ 
    eigenVectors;

  eigenVectors = 
   Transpose[
    Developer`ToPackedArray[eigenVectors[[All, #]] & /@ spans], {2, 1,
      3}];
  eigenVectorsIF = Table[{}, {n}, {Length[spans]}];
  Do[eigenVectorsIF[[i, j]] = 
    NDSolve`FEM`ElementMeshInterpolation[{mesh}, 
     eigenVectors[[i, j]]], {i, n}, {j, Length[spans]}];
  {eigenValues, eigenVectorsIF}
  ]

What it does it uses the method option "ConstraintMethod" -> "Remove" the eliminate the Dirichlet boundary positions and the returns the data structure that contains the Dirichlet information. It's better the use Eigensystem with the Arnoldi method to get the speed in place of the LinearSolve variant. A bit of shifting the amount of the eigenvalues and vecors that need to be computed is done and then the Dirichlet conditions are then reinserted again.

(1) Many thanks, once again, to user21 for the full and helpful response. My code now works and the alternative of NDEigensystem is possibly where I should be looking. I completely agree that my theory results (included for comparison above) is only approximate and I should have made this clear in my question. The Bernoulli beam theory does not include Poission ratio effects, shear deformation and rotary inertia. The finite element method should be better in this respect. The new option value "Remove" is what was needed. Finding these options is a challenge. Thanks again — Hugh, Nov 22, 2015 at 11:30
(2) I will investigate NDEigensystem further. For the examples you give am I correct in noting that the assumption made is that the mass matrix is not only diagonal but the identity matrix? I will investigate circumstances where this is not the case and see if I can work things out. A look at the mass/damping matrix suggested that it was not diagonal. An alternative assumption would be that a uniform density of 1.0 is associated with the inertia term. — Hugh, Nov 22, 2015 at 12:01
@Hugh, no the system matrices are not lumped. The system matrices from the code above are the same as what NDEigensystem uses. You can verify that by inspecting the system matrices of NDEigensystem via the documented "VectorNormalization" method option. Now, the fact that "Remove" (and a few other other constraint deployment methods) are not documented is that they do not yet work across the board. Once they do, they will be documented. Promise. — user21, Nov 23, 2015 at 8:34

Stack Exchange Network

Finite element boundary breaking

1 Answer 1

Your Answer

Not the answer you're looking for? Browse other questions tagged
differential-equations
finite-element-method
or ask your own question.

Linked

Hot Network Questions

Finite element boundary breaking

1 Answer 1

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged differential-equationsfinite-element-method or ask your own question.

Linked

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
differential-equations
finite-element-method
or ask your own question.