Test a wooden board's vibration mode

Question

Here is a wooden board, with dimensions shown on the picture below. How we can use Mathematica's newly build-in finite element analysis features to show the different modes of its vibrations. Assuming the board is made of spruce.

Answer 1

The only reason I am attempting to answer this is to perhaps get a Reversal badge. There you go...

We will go slowly and this answer is the basis for what comes next. Let's start with two dimensions. You'll see why.

We create a rectangular region:

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[FullRegion[2], {{0, 5}, {0, 1}}, "MeshOrder" -> 1, "MaxCellMeasure" -> 0.005];

We use the plane stress equations of a material sheet. For every material you need to set the Young's modulus and Poisson's ratio.

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

Make a coupled, time dependent PDE and constrain the region to not move at the left boundary.

pde2D = D[{u[t, x, y], v[t, x, y]}, t] + planeStress == {0, 0};
(* held fixed at left *)
bcs = DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, x==0];

What follows next, is pretty much the same from the other post.

(*We use NDSolve as a pre-processor:*)
{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"}}}];

(*Extract the finite element data:*)
femdata = state["FiniteElementData"];

initBCs = femdata["BoundaryConditionData"];
methodData = femdata["FEMMethodData"];
initCoeffs = femdata["PDECoefficientData"];

(*Set up the solution and variable data:*)
vd = methodData["VariableData"];
nr = ToNumericalRegion[mesh];
sd = NDSolve`SolutionData[{"Space" -> nr, "Time" -> 0.}];

(*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"];

(*Deploy the boundary conditions:*)
DeployBoundaryConditions[{load, 
  stiffness, damping}, discreteBCs]

(*Set the number of X smallest eigen values we would like to compute \
but ignore the Dirichlet positions.*)

nDiri = If[Length[#] > 0, First[#], 0] &[
   Dimensions[discreteBCs["DirichletMatrix"]]];

numEigenToCompute = 5;
numEigen = numEigenToCompute + nDiri;

Now, things become hard. We solve the eigensystem.

(*Solve the eigen system: this is how you should do it*)
res = Eigensystem[{stiffness, damping}, -numEigen];

You will need patience. Play with the "Arnoldi" method and the shift. (left as an exercise). As a bad alternative you can play with

(* this may be a bit faster but is the dark side... *)
(*
mm=LinearSolve[damping,stiffness];
res=Eigensystem[mm,-numEigen];
*)

Further down I use "FEAST" as a solve.

For each eigenvalue we now have two eigenvectors. Once in the x-direction and one in the y-direction. So we post-process:

res = Reverse /@ res;
eigenValues = res[[1, nDiri + 1 ;; Abs[numEigen]]];
eigenVectors = res[[2, nDiri + 1 ;; Abs[numEigen]]];
(*res=Null;*)
inciOffs = methodData["IncidentOffsets"];
 spans = MapThread[Span, {Most[inciOffs] + 1, Rest[inciOffs]}];
 eigenVectors = 
  Transpose[
   Developer`ToPackedArray[eigenVectors[[All, #]] & /@ spans], {2, 1, 
    3}];
 eigenVectorsIF = Table[{}, {numEigenToCompute}, {Length[spans]}];
 Do[
  eigenVectorsIF[[i, j]] = 
   NDSolve`FEM`ElementMeshInterpolation[{mesh}, 
    eigenVectors[[i, j]]] 
  , {i, numEigenToCompute}, {j, Length[spans]}];
res = {eigenValues, eigenVectorsIF};

And visualize the first 5 eigenvectos in the x-direction:

Show[NDSolve`FEM`ElementMeshPlot3D[res[[2, #, 1]]["Coordinates"][[1]],
     NDSolve`FEM`ElementMeshDirective -> 
     Directive[EdgeForm[Gray], FaceForm[]]], 
   NDSolve`FEM`ElementMeshPlot3D[res[[2, #, 1]]], Boxed -> False, 
   Axes -> False, ImageSize -> 600] & /@ Range[numEigenToCompute];

Visualize the y-direction:

Show[NDSolve`FEM`ElementMeshPlot3D[res[[2, #, 2]]["Coordinates"][[1]],
     NDSolve`FEM`ElementMeshDirective -> 
     Directive[EdgeForm[Gray], FaceForm[]]], 
   NDSolve`FEM`ElementMeshPlot3D[res[[2, #, 2]]], Boxed -> False, 
   Axes -> False, ImageSize -> 600] & /@ Range[numEigenToCompute];

And the "breathing modes" enlarged by a factor. (I just invented that name - it might mean something else)

fact = 5;
Show[{
    NDSolve`FEM`ElementMeshPlot3D[res[[2, 2, 1]]["Coordinates"][[1]], 
     NDSolve`FEM`ElementMeshDirective -> 
      Directive[EdgeForm[Gray], FaceForm[]]],
    NDSolve`FEM`ElementMeshPlot3D[
     NDSolve`FEM`ElementMeshInterpolation[
      res[[2, #, 1]]["Coordinates"], 
      fact*Sqrt[Total[#["ValuesOnGrid"]^2 & /@ res[[2, #]]]]]]
    }, Boxed -> False, ImageSize -> 600] & /@ Range[numEigenToCompute]

The 3D case. Unfortunately, my local super-computer center is closed, but here is how to do it. Create a mesh and enlarge the features for now to create not too many elements.

base = {0, 0, 0};
h1 = 5;
h2 = 5;
w1 = 40;
l1 = 76;
cw1 = 5;
cl1 = 68;
cw2 = 36;
cl2 = 5;
offset1 = base + {(w1 - cw1)/2, (l1 - cl1)/2, 0};
offset2 = base + {(w1 - cw2)/2, (l1 - cl2)/2, 0};
offset3 = base + {(w1 - cw1)/2, (l1 - cl2)/2, 0};
ClearAll[rect]
rect[base_, w_, l_, h_] := {base + {0, 0, h}, base + {w, 0, h}, 
  base + {w, l, h}, base + {0, l, h}}
coords = ConstantArray[{0., 0., 0.}, 4 + 4 + 12 + 12];
coords[[{1, 2, 3, 4}]] = rect[base, w1, l1, 0];
coords[[{5, 6, 7, 8}]] = rect[base, w1, l1, h1];
coords[[{9, 10, 15, 16}]] = rect[offset1, cw1, cl1, h1];
coords[[{19, 12, 13, 18}]] = rect[offset2, cw2, cl2, h1];
coords[[{20, 11, 14, 17}]] = rect[offset3, cw1, cl2, h1];
coords[[20 + Range[12]]] = ({0, 0, h2} + #) & /@ 
   coords[[8 + Range[12]]];
bmesh = ToBoundaryMesh["Coordinates" -> coords, 
  "BoundaryElements" -> {QuadElement[{{1, 2, 3, 4}, {1, 2, 6, 5}, {2, 
       3, 7, 6}, {3, 4, 8, 7}, {4, 1, 5, 8},
      {5, 6, 10, 9}, {6, 12, 11, 10}, {6, 7, 13, 12}, {7, 15, 14, 
       13}, {7, 8, 16, 15}, {8, 18, 17, 16}, {8, 5, 19, 18}, {5, 9, 
       20, 19},
      Sequence @@ ({{9, 10, 11, 20}, {11, 12, 13, 14}, {14, 15, 16, 
           17}, {17, 18, 19, 20}, {20, 11, 14, 17}} + 12),
      Sequence @@ (Partition[Join[Range[9, 20]], 2, 1, 
          1] /. {i1_, i2_} :> {i1, i2, i2 + 12, i1 + 12})
      }]}]

If you want to visualize the boundary structure:

Show[
 bmesh["Wireframe"],
 bmesh["Wireframe"["MeshElement" -> "PointElements", 
   "MeshElementIDStyle" -> Red]]
 ];

Create the mesh:

mesh = ToElementMesh[bmesh, "MeshOrder" -> 1, 
   "MaxCellMeasure" -> 10];
mesh["Wireframe"]

Here is the PDE stress operator in 3D:

stressOperator[
  Y_, ν_] := {Inactive[
     Div][{{0, 0, -((Y*ν)/((1 - 2*ν)*(1 + ν)))}, {0, 0, 
       0}, {-Y/(2*(1 + ν)), 0, 0}}.Inactive[Grad][
      w[t, 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[t, 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[t, 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[t, 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[t, 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[t, 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[t, 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[t, 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[t, x, y, z], {x, y, z}], {x, y, z}]}

And the 3D PDE. no boundary conditions, it seems an unconstraint analysis is wanted.

(* choose your Y and ν -- no idea what the values for spruce are - 
   is it wet, dry, old, ...? *)
pde3D = D[{u[t, x, y, z], v[t, x, y, z], w[t, x, y, z]}, t] + 
    stressOperator[100, 1/3] == {0, 0, 0};
(* unconstraint? Yes! *)
bcs = Sequence[];

If you want a constraint analysis then you'd have to set DirichletCondition on the boundary (bcs = DirichletCondition[{u[x,y,z]==0,v[x,y,z]==0,w[x,y,z]==0},True])

Now, we use

{state} = 
  NDSolve`ProcessEquations[{pde3D, bcs, u[0, x, y, z] == 0, 
    v[0, x, y, z] == 0, w[0, x, y, z] == 0}, {u, v, w}, {t, 0, 
    1}, {x, y, z} \[Element] mesh, 
   Method -> {"PDEDiscretization" -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"FiniteElement"}}}];

for pre-processing. There may be a waring about no DirichletCondition or no NeumannValue; this is save to ignore in this case.

And now the same as above - it will take a long time.... I did not want to wait. (Also I did not want to think about how to visualize this in 3D... that's for you...)

When you do this do not forget that the result need to be sorted and post processed as in the 2D example above. i.e.

res = Reverse /@ res;
eigenValues = res[[1, nDiri + 1 ;; Abs[numEigen]]];
eigenVectors = res[[2, nDiri + 1 ;; Abs[numEigen]]];

Update:

The Eigensystem solution above takes about 450 seconds on my machine. You can use

AbsoluteTiming[
 res = Eigensystem[{stiffness, damping}, -numEigen, 
    Method -> {"FEAST", "Tolerance" -> 10^-6}];]

to get it down to 45 seconds. Which is a bit better.

Here are the deformations for the eigenmodes 7 to 10 - the first 6 are zero. (I must admit that that I am not sure if that makes sense)

res[[1]]
{0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.011403583383327644`, \
0.01526089137692353`, 0.05661022352859022`, 0.07266104128273859`}

And the visualizations:

MeshRegion[
   ElementMeshDeformation[mesh, res[[2, #]], 
    "ScalingFactor" -> 100]] & /@ Range[7, numEigenToCompute]

When you run the 3D exmaple, please adjust the numEigenToCompute to be appropriate.

Anything else?

Answer 2

This is not another answer but an attempt to repeat the calculations in the other answer in version 11. There are problems at the end for which I need help (edit: now sorted and corrected below). Note that there is a more direct method using NDEigensystem which is discussed here. There are issues with the eigen solver which may take a long time; this is discussed above and in the referenced question.

I am also writing this answer as it provides a workflow for solving a vibration eigenvalue problem where details need to be extracted. I make some comments that may help those used to doing vibration calculations.

Set up the stress operator for a general structural, time dependent systems. (Note that you can't just type in the textbook equations because there needs to be appropriate Inactive commands).

Needs["NDSolve`FEM`"]


ClearAll[stressOperator, u, v, w, t, x, y, z]; 
stressOperator[
  Y_, ν_] := {Inactive[
     Div][{{0, 0, -((Y*ν)/((1 - 2*ν)*(1 + ν)))}, {0, 0, 
       0}, {-Y/(2*(1 + ν)), 0, 0}}.Inactive[Grad][
      w[t, 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[t, 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[t, 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[t, 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[t, 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[t, 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[t, 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[t, 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[t, x, y, z], {x, y, z}], {x, y, z}]}

Now make the geometry

 base = {0, 0, 0};
    h1 = 5;
    h2 = 5;
    w1 = 40;
    l1 = 76;
    cw1 = 5;
    cl1 = 68;
    cw2 = 36;
    cl2 = 5;
    offset1 = base + {(w1 - cw1)/2, (l1 - cl1)/2, 0};
    offset2 = base + {(w1 - cw2)/2, (l1 - cl2)/2, 0};
    offset3 = base + {(w1 - cw1)/2, (l1 - cl2)/2, 0};
    ClearAll[rect]
    rect[base_, w_, l_, h_] := {base + {0, 0, h}, base + {w, 0, h}, 
      base + {w, l, h}, base + {0, l, h}}
    coords = ConstantArray[{0., 0., 0.}, 4 + 4 + 12 + 12];
    coords[[{1, 2, 3, 4}]] = rect[base, w1, l1, 0];
    coords[[{5, 6, 7, 8}]] = rect[base, w1, l1, h1];
    coords[[{9, 10, 15, 16}]] = rect[offset1, cw1, cl1, h1];
    coords[[{19, 12, 13, 18}]] = rect[offset2, cw2, cl2, h1];
    coords[[{20, 11, 14, 17}]] = rect[offset3, cw1, cl2, h1];
    coords[[20 + Range[12]]] = ({0, 0, h2} + #) & /@ 
       coords[[8 + Range[12]]];
    bmesh = ToBoundaryMesh["Coordinates" -> coords, 
       "BoundaryElements" -> {QuadElement[{{1, 2, 3, 4}, {1, 2, 6, 5}, {2,
             3, 7, 6}, {3, 4, 8, 7}, {4, 1, 5, 8}, {5, 6, 10, 9}, {6, 12, 
            11, 10}, {6, 7, 13, 12}, {7, 15, 14, 13}, {7, 8, 16, 15}, {8, 
            18, 17, 16}, {8, 5, 19, 18}, {5, 9, 20, 19}, 
           Sequence @@ ({{9, 10, 11, 20}, {11, 12, 13, 14}, {14, 15, 16, 
                17}, {17, 18, 19, 20}, {20, 11, 14, 17}} + 12), 
           Sequence @@ (Partition[Join[Range[9, 20]], 2, 1, 
               1] /. {i1_, i2_} :> {i1, i2, i2 + 12, i1 + 12})}]}];

Show the geometry

Show[bmesh["Wireframe"], 
 bmesh["Wireframe"["MeshElement" -> "PointElements", 
   "MeshElementIDStyle" -> Red]]]

mesh = ToElementMesh[bmesh, "MeshOrder" -> 1, "MaxCellMeasure" -> 10];
mesh["Wireframe"]

Set up the pde and give values for modulus of elasticity and Poisson ratio. Those users familiar to vibration may be surprised that the derivative with respect to time is only first order and not second (Newton's law is second). This is because the equations are converted to first order within the code and we are just trying to get the matrices. For this reason what is normally called the mass matrix is called the damping matrix here.

Y = 100;
ν = 1/3;
pde3D = D[{u[t, x, y, z], v[t, x, y, z], w[t, x, y, z]}, t] + 
    stressOperator[Y, ν] == {0, 0, 0};

Define the boundary conditions

bcs = Sequence[];

The equations are now processed.

{state} = 
  NDSolve`ProcessEquations[{pde3D, bcs, u[0, x, y, z] == 0, 
    v[0, x, y, z] == 0, w[0, x, y, z] == 0}, {u, v, w}, {t, 0, 
    1}, {x, y, z} ∈ mesh, 
   Method -> {"PDEDiscretization" -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"FiniteElement"}}}];

The finite element data are now extracted.

femdata = state["FiniteElementData"];
initBCs = femdata["BoundaryConditionData"];
methodData = femdata["FEMMethodData"];
initCoeffs = femdata["PDECoefficientData"];

The solution and variable data are extracted

vd = methodData["VariableData"];
nr = ToNumericalRegion[mesh];
sd = NDSolve`SolutionData[{"Space" -> nr, "Time" -> 0.}];

The PDE and boundary conditions are discretized

discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];

The system matrices can now be extracted

load = discretePDE["LoadVector"];
stiffness = discretePDE["StiffnessMatrix"];
damping = discretePDE["DampingMatrix"];

The system matrices must now be modified to account for the boundary conditions. This is done "in place" so there are no new matrices

DeployBoundaryConditions[{load, stiffness, damping}, discreteBCs];

Now we can calculate the eigen system. I look for the smallest 10 frequencies (this is the -10). They are given in order of largest to smallest so they need to be reversed to get them from small to large. These are actually the frequencies squared because of the formulation as a standard eigenvalue problem. Take the square root to get the frequencies and divide by 2 Pi to get them in Hz.

AbsoluteTiming[res = Eigensystem[{stiffness, damping}, -10,
    Method -> {"FEAST", "Tolerance" -> 10^-6}];]

This takes about 25 seconds on my machine. Here are the results

TableForm[Reverse@res[[1]], TableHeadings -> {Automatic, None}]

These agree with the original post so all is good. (This is edited from my original post, see comments, because I forgot to revers.) The first six are zero, as they should be, because there are no constraints so there are six modes of rigid body vibration at zero frequency.