First to fifth order natural frequencies of a cantilever beam with known material constants

Question

We assume that the homogeneous Euler-Bernoulli cantilever with constant cross-section length is 1, the material constant EI is 1 and ρS is 1.

According to the constraints, we can establish the following differential equations:

DSolve[{D[D[ω[x, t], {x, 2}], {x, 2}] + 
    D[ω[x, t], {t, 2}] == 0, ω[0, 0] == 
   0, ω[x, 0] == 0, ω[0, t] == 0, 
  D[ω[0, t], {t, 1}] == 
   0(*Constraints on simply supported ends*), 
  D[ω[1, t], {t, 2}] == 0, 
  D[ω[1, t], {t, 3}] == 
   0(*Constraints on fixed end*)}, ω[x, t], {x, t}]

But I still can't find the only analytical solution, while the textbook can determine the natural frequencies of each stage of the cantilever according to these conditions (In addition, I find that $\rho$ and $S$ cannot be replaced effectively in the following code).

Table[(((2. i - 1) π)/(
    2 l))^2 Sqrt[(\[DoubleStruckCapitalE]*\[DoubleStruckCapitalI])/(\
ρ*S)], {i, 1, 10}] /. {l -> 
   1., \[DoubleStruckCapitalE]*\[DoubleStruckCapitalI] -> 
   1., ρ -> 1., S -> 1}

How can I solve the previous differential equation to get similar frequency results?

When I use the relevant code of this post, I can't even get the numerical solution of the cantilever vibration:

ClearAll["Global`*"]
tau = 10;
L = 1;
Elastic = 1;
Imoment = 1;
ρ = 1;
S = 1;
sol = NDSolveValue[{D[Elastic*Imoment*D[w[x, t], {x, 2}], {x, 2}] + 
     S*ρ*D[w[x, t], {t, 2}] == 0, w[0, t] == 0, 
   w[x, 0] == x^2/6 (3 - x), Derivative[0, 1][w][0, t] == 0, 
   Derivative[0, 2][w][L, t] == Derivative[0, 3][w][L, t] == 0}, 
  w, {x, 0, L}, {t, 0, tau}, 
  Method -> {"MethodOfLines", 
    "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 100},
     "SpatialDiscretization" -> {"TensorProductGrid", 
      "MaxPoints" -> 100, "MinPoints" -> 100, 
      "DifferenceOrder" -> 2}}, MaxSteps -> 10^6]

Answer 1

To find eigenmodes, time is not involved. Following Wikipeida You just need to solve an ODE and not a PDE

---

---

This below is the code to reproduce the eigenmodes shown at the above page.

I tried to use Mathematica's NDEigesystem on this, but I could not make it do it, so did it by "hand".

The trick to finding eigenvalues is not to put all 4 boundaries conditions in at once, else Mathematica will return trivial solution of course, since all BC are zero and there is no load (it is free vibration). So we put 3 BC. You pick which 3 to put in and which one to keep out. Below I kept out the last one. (the third derivative on the free end).

After you get solution, it will have one constant of integration in it. Then by inspection, we can see now the condition to make the solution non-trivial. Now we set up an equation (it will be nonlinear) to solve for eigenvalues (Letting length of beam = 1).

Now we find the eigenvalues, going back to the ODE and plugging each eigenvalue at a time and putting back now all the 4 BC in. This results in corresponding eigenmode solution. To plot them all on same plot, I normalized them all to 1.

This is the final plot

This is what Wikipedia gives

(ps. For some reason, I had to flip the sign on some of the eigenmodes below to get same plot as Wikipedia. otherwise it was upside down. I have not found out why).

Here is a diagram of the BC (diagram thanks to this page)

For a cantilevered beam, the boundary conditions are as follows: (Thanks to this page)

---

w(0)=0 . This boundary condition says that the base of the beam (at the wall) does not experience any deflection. w'(0)=0 . We also assume that the beam at the wall is horizontal, so that the derivative of the deflection function is zero at that point. w''(L)=0 . This boundary condition models the assumption that there is no bending moment at the free end of the cantilever. w'''(L)=0 . This boundary condition models the assumption that there is no shearing force acting at the free end of the beam.

---

Here is the code. This solve the ODE shown above. Notice it has 3 BC and not 4.

ClearAll[w, wHat, x, L];
pde = D[wHat[x], {x, 4}] == w^2 wHat[x];
leftEndBc = {wHat[0] == 0, wHat'[0] == 0};
rightEndBc = {wHat''[L] == 0};
sol = wHat[x] /. First@DSolve[{pde, leftEndBc, rightEndBc}, wHat[x], x]

This is what DSolve gives

$$ -\frac{c_1 e^{-\sqrt{w} x} \left(e^{L \sqrt{w}+2 \sqrt{w} x} \sin \left(L \sqrt{w}\right)+e^{2 L \sqrt{w}+\sqrt{w} x} \sin \left(\sqrt{w} x\right)-e^{L \sqrt{w}+2 \sqrt{w} x} \cos \left(L \sqrt{w}\right)-e^{2 L \sqrt{w}+\sqrt{w} x} \cos \left(\sqrt{w} x\right)+2 e^{L \sqrt{w}+\sqrt{w} x} \cos \left(L \sqrt{w}\right) \sin \left(\sqrt{w} x\right)-2 e^{L \sqrt{w}+\sqrt{w} x} \sin \left(L \sqrt{w}\right) \cos \left(\sqrt{w} x\right)+e^{2 L \sqrt{w}}+e^{L \sqrt{w}} \sin \left(L \sqrt{w}\right)+e^{L \sqrt{w}} \cos \left(L \sqrt{w}\right)-e^{2 \sqrt{w} x}+e^{\sqrt{w} x} \sin \left(\sqrt{w} x\right)+e^{\sqrt{w} x} \cos \left(\sqrt{w} x\right)\right)}{e^{2 L \sqrt{w}}+2 e^{L \sqrt{w}} \sin \left(L \sqrt{w}\right)-1} $$

This makes up equation to satisfy the missing BC at right end (the one left out)

  eq = FullSimplify[(D[sol, {x, 3}] /. x -> L)]

$$ \frac{2 c_1 w^{3/2} \left(\cos \left(L \sqrt{w}\right) \cosh \left(L \sqrt{w}\right)+1\right)}{\sin \left(L \sqrt{w}\right)+\sinh \left(L \sqrt{w}\right)} $$

Here is the "manual" step. Looking at the above, since the above is zero (it is the last B.C.), then for non-trivial solution we do not want $c_1=0$, then only other choice to make the above zero is that

$$ \cos \left(L \sqrt{w}\right) \cosh \left(L \sqrt{w}\right)+1=0 $$

Now NSolve is used to find eigenvalues $w$ (solutions of the above equation). We can use $L=1$ now. Then the equation to solve for frequencies is the above.

frequencies = w /. NSolve[1 + Cos[Sqrt[w]] Cosh[Sqrt[w]] == 0 && 0 < w < 200, w, Reals]

gives

  {3.51602, 22.0345, 61.6972, 120.902, 199.86}

These are the first 5 frequencies. Now for each, we solve the ODE again.

pde = D[wHat[x], {x, 4}] == frequencies[[1]]^2 *wHat[x];
L = 1;
leftEndBc = {wHat[0] == 0, wHat'[0] == 0};
rightEndBc = {wHat''[L] == 0, wHat'''[L] == 0};
sol = NDSolveValue[{pde, leftEndBc, rightEndBc}, wHat, {x, 0, 1}];
max = Maximize[{Abs[sol[x]], 0 < x < 1}, x][[1]];
p1 = Plot[Callout[-sol[x]/max, "mode 1", 0.7], {x, 0, 1}, PlotStyle -> Blue]

pde = D[wHat[x], {x, 4}] == frequencies[[2]]^2 *wHat[x];
L = 1;
leftEndBc = {wHat[0] == 0, wHat'[0] == 0};
rightEndBc = {wHat''[L] == 0, wHat'''[L] == 0};
sol = NDSolveValue[{pde, leftEndBc, rightEndBc}, wHat, {x, 0, 1}];
max = Maximize[{Abs[sol[x]], 0 < x < 1}, x][[1]];
p2 = Plot[Callout[-sol[x]/max, "mode 2", .4, .4], {x, 0, 1}, PlotStyle -> Magenta]

pde = D[wHat[x], {x, 4}] == frequencies[[3]]^2 *wHat[x];
L = 1;
leftEndBc = {wHat[0] == 0, wHat'[0] == 0};
rightEndBc = {wHat''[L] == 0, wHat'''[L] == 0};
sol = NDSolveValue[{pde, leftEndBc, rightEndBc}, wHat, {x, 0, 1}];
max = Maximize[{Abs[sol[x]], 0 < x < 1}, x][[1]];
p3 = Plot[Callout[-sol[x]/max, "mode 3", .2], {x, 0, 1}, PlotStyle -> Brown]

pde = D[wHat[x], {x, 4}] == frequencies[[4]]^2 *wHat[x];
L = 1;
leftEndBc = {wHat[0] == 0, wHat'[0] == 0};
rightEndBc = {wHat''[L] == 0, wHat'''[L] == 0};
sol = NDSolveValue[{pde, leftEndBc, rightEndBc}, wHat, {x, 0, 1}];
max = Maximize[{Abs[sol[x]], 0 < x < 1}, x][[1]];
p4 = Plot[Callout[sol[x]/max, "mode 4", Below], {x, 0, 1}, PlotStyle -> {Thick, Green}]


 Show[{p1, p2, p3, p4}, PlotRange -> {{0, 1}, {-1.2, 1.2}}]

Answer 2

I think Bill Watts answer is the correct approach. Just to add a bit more that can be generalised to other boundary conditions. Note that this is an eigenvalue problem. We are looking for frequencies that make the differential equation work. Thus we ought to form a determinant and look for roots.

Starting with Bill's equation but without the greek

X[x_] := c1 Cos[b x] + c2 Sinh[b x] + c3 Sin[b x] + c4 Cosh[b x]

We now introduce the boundary conditions for each end of the beam

  eqns = {
  X[0] == 0,  (* No displacment at x = 0 *)
  (D[X[x], x] /. x -> 0) == 0, (* No slope at x = 0 *)
  (D[X[x], {x, 2}] /. x -> L) == 0, (* No bending moment at x = L *)
  (D[X[x], {x, 3}] /. x -> L) == 0 (* No shear force at x = L *)
  }

We can make this a matrix equation as follows. First we make a vector out of the unknown coefficients c1, c2, c3, c4 and then we use CoefficientArrays to get a matrix

vec = {c1, c2, c3, c4};
{rhs, mat} = Normal[CoefficientArrays[eqns, vec]];
MatrixForm[mat]

We have a matrix equation in the form mat.vec == 0. So either we have the trivial solution vec == 0 or the determinant of the matrix must be zero. This is the standard eigenvalue problem. Taking the determinant we get

det = Det[mat] // Simplify

(*2 b^6 (1 + Cos[b L] Cosh[b L]) *)

Which is the same as Bills solution. We can see that because of the cosh term the roots are going to be approximately (2 n - 1) pi/2.

First we tidy up the determinant by combining L and b.

det1 = det/b^6 /. b -> Lb/L

Now we can get the roots of the determinate using the approximate solution as a starting point.

 roots = {b -> Lb/L} /. 
  Table[FindRoot[det1 == 0, {Lb, (2 n - 1) \[Pi]/2}], {n, 5}]

    (* {{b -> 1.8751/L}, {b -> 4.69409/L}, {b -> 7.85476/L}, {b -> 10.9955/
   L}, {b -> 14.1372/L}} *)

The eignevectors are found by determining the NullSpace of the matrix when each root is substituted in. We can do this as follows

evecs = Thread[vec -> #] & /@ (NullSpace[#][[1]] & /@ (mat /. roots));

Now I join all the root and corresponding eigenvector in one set of replacement rules.

valsVecs = Join[#[[1]], #[[2]]] & /@ Transpose[{roots, evecs}];

Here are the mode shapes

Plot[Evaluate[X[L x] /. valsVecs], {x, 0, 1}]

By applying other boundary conditions one can deal with all beam configurations.

Answer 3

Slightly different way for constant EI and ρS.

Start with the pde.

pde = D[EI*D[w[x, t], x, x], x, x] + \[Rho]S*D[w[x, t], t, t] == 0

Assume sinusoidal in time.

w[x_, t_] = X[x] Sin[ω t]

which simplifies the pde

pde
(*EI*Derivative[4][X][x]*Sin[t*ω] - ρS*ω^2*X[x]*Sin[t*ω] == 0*)

dsol = DSolve[pde, X[x], x] // Flatten

Assign the constants to a new constant, β and solve for ω

ωRule = Solve[(ρS^(1/4) Sqrt[ω])/EI^(1/4) == β, ω] // Flatten

Clean up a bit

X[x_] = X[x] /. dsol /. {C[1] -> c1, C[2] -> c2, C[3] -> c3, 
     C[4] -> c4} /. ωRule // PowerExpand

(*c1 Cos[β x] + c2 E^(β (-x)) + c3 Sin[β x] + c4 E^(β x)*)

Get in Trig form.

X[x] // ExpToTrig // Collect[#, {Sinh[β x], Cosh[β x]}] &
(*c1 Cos[β x] + (c4 - c2) Sinh[β x] + (c2 + c4) Cosh[β x] + c3 Sin[β x]*)

X[x_] = % /. {c4 - c2 -> c2, c2 + c4 -> c4}
(*c1 Cos[β x] + c2 Sinh[β x] + c3 Sin[β x] + c4 Cosh[β x]*)

Apply some boundary conditions for a cantilever. The displacement is zero at the fixed end x = 0.

X[0] == 0
(*c1+c4==0*)

c4 = c4 /. Solve[%, c4][[1]]

Derivative of displacement is zero at the fixed end x = 0.

(D[X[x], x] /. x -> 0) == 0

c3 = c3 /. Solve[%, c3][[1]]

The Shear is zero at the free end.

(D[X[x], x, x, x] /. x -> L) == 0 //Simplify

c2 = c2 /. Solve[%, c2][[1]]

Bending Moment is zero at the free end.

(D[w[x, t], x, x] /. x -> L) == 0//Simplify

(β c1 Sin[t ω] (Cos[β L] Cosh[β L] + 1))/(
 Cos[β L] + Cosh[β L]) == 0

This gives us the condition for β.

Use βL so we can divide by L later.

F[βL_] = Cos[βL] Cosh[βL] + 1

Find the first five roots of βL with FindRoot. I didn't have much luck with NSolve. Plot to get starting values.

Plot[Evaluate[F[βL]], {βL, 0, 10}, PlotRange -> {-10, 10}]

n = 5;
a = FindRoot[F[α] == 0, {α, 2}]
βL[1] = α /. a // N;
a = FindRoot[F[α] == 0, {α, 4.8}]
βL[2] = α /. a // N;
For[i = 3, i <= n,
    a = FindRoot[ F[α] == 0, {α, 2 βL[i - 1] - βL[i - 2]}];
   βL[i] = α /. a // N; 
 i++]

The first five frequencies.

Table[ω /. ωRule /. β -> βL[i]/L, {i, 5}]
{(3.51602 Sqrt[EI])/(L^2 Sqrt[ρS]), (22.0345 Sqrt[EI])/(
 L^2 Sqrt[ρS]), (61.6972 Sqrt[EI])/(L^2 Sqrt[ρS]), (
 120.902 Sqrt[EI])/(L^2 Sqrt[ρS]), (199.86 Sqrt[EI])/(
 L^2 Sqrt[ρS])}

You don't have to plug in numbers for the frequencies until the end this way. We have found all the constants except for c1 which is found from the initial conditions of the beam and orthogonality. The displacement of beam w = X[x] Sin[ω t] is found by an infinite sum over all the ω`s.