Transverse beam vibration: Dynamic Euler-Bernoulli equation

[alansaillet]

I want to simulate the transverse vibration of a beam. I used the standard Dynamic Euler-Bernoulli beam equation and clamped-clamped boundary conditions. The solution seems to be OK for small values of t (time), but becomes unstable afterwards. I modeled it successfully in Python using Euler-Implicit method, and I do not have that kind of behavior, which is not coherent with the reality.

Equation:

I used the standard Dynamic Euler-Bernoulli beam equation:

which translates in

PDE = D[y[x, t], {t, 2}] + c^2*D[y[x, t], {x, 4}] == 0

For 0<t<4, I obtain:

For 0<t<8, I obtain:

Do you have an idea of the source of the problem?

FULL CODE

(*CONSTANTS*)
c = 5;
L = 2*Pi;
T = 8;

(*EQUATION*)
PDE = D[y[x, t], {t, 2}] + c^2*D[y[x, t], {x, 4}] == 0

(*INITIAL CONDITIONS*)
IC1 = y[x,0]== Sin[x];(*Initial shape of the beam*)
IC2 = D[y[x, t],{t, 1}] == 0 /. {t -> 0};(*No initial velocity*)

(*BOUNDARY CONDITIONS*)
(*Clamped on the left side of the beam*)
BCLeft1 = y[x, t] == 0 /. {x -> 0};
BCLeft2 = D[y[x, t],{x, 1}] == 0 /. {x -> 0};
(*Clamped on the right side of the beam*)
BCRight1 = y[x, t] == 0 /. {x -> L};
BCRight2 = D[y[x, t],{x, 1}]== 0 /. {x -> L};

(*SOLVE*)
soln = NDSolve[{PDE,IC1,IC2,BCLeft1,BCLeft2,BCRight1,BCRight2}, y[x, t], {x, 0,L},{t,0,T},Method -> {"MethodOfLines", 
         "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 0}}][[1,1,2]];
Plot3D[soln, {x, 0,L},{t,0,T},PlotPoints->150]

Dynamics.nb (11.9 MB)

[mariuszi28cg (Mariusz Iwaniuk)]

Maybe this helps.

[alansaillet]

Thank you for the link! Very interesting to see the different approaches

In my case, I want a Fully NDSolve-based Numeric Solution
Unfortunately, the solution proposed gives the exact same results (I believe these are the default parameters or the ScaleFactor of DifferentiateBoundaryConditions is set to 0 when inconsistent boundary conditions are detected

[anon_91240867]

I do not have a thorough citation or reference book for modeling (rigid? or elastic?) beam vibration using (rigid? or elastic?) holding. But I believe your BC are the issue.

[anon_91240867]

When NDSolve checks the boundaries for u and u’ and they are not as specified, it balks. y[x, 0] == Sin[x]. But y[0,t] == 0. But D[y[x,t],x] == Cos[ x] and will not be 0 (if i read the following correclty, and my reading of it was extremely brief).

The Numerical Method of Lines—Wolfram Documentation

(I remind that Mariusz suggestion should be followed through with)

[alansaillet]

You are right, I tried various BCs, and finally made a mistake : I used Clamped-Clamped instead of Pinned-Pinned but used the Pinned-Pinned BCs. Nevertheless, it is important to keep T in the order of magnitude of 1 to see some results.

Pinned-Pinned Results: OK

(*INITIAL CONDITIONS*)
IC1 = y[x,0]== Sin[x];(*Initial shape of the beam*)
IC2 = D[y[x, t],{t, 1}] == 0 /. {t -> 0};(*No initial velocity*)
(*BOUNDARY CONDITIONS*)
(*Pinned on the left side of the beam*)
BCLeft1 = y[x, t] == 0 /. {x -> 0};
BCLeft2 = D[y[x, t],{x, 2}] == 0 /. {x -> 0};
(*Pinned on the right side of the beam*)
BCRight1 = y[x, t] == 0 /. {x -> L};
BCRight2 = D[y[x, t],{x, 2}]== 0 /. {x -> L};

Note: I still get the message

Warning: boundary and initial conditions are inconsistent.

Clamped-Clamped:

(*INITIAL CONDITIONS*)
IC1 = y[x,0]== Cos[x]-1;(*Initial shape of the beam*)
IC2 = D[y[x, t],{t, 1}] == 0 /. {t -> 0};(*No initial velocity*)

(*BOUNDARY CONDITIONS*)
(*Clamped on the left side of the beam*)
BCLeft1 = y[x, t] == 0 /. {x -> 0};
BCLeft2 = D[y[x, t],{x, 1}] == 0 /. {x -> 0};
(*Clamped on the right side of the beam*)
BCRight1 = y[x, t] == 0 /. {x -> L};
BCRight2 = D[y[x, t],{x, 1}]== 0 /. {x -> L};

For this one, I still have problems, even though BCs seem to be correct.
I will check the alternative solutions you gave. How to retrieve the articles you quoted (e.g. farlow)?

EDIT: the solution was to change the solve method to either

ppR = 25
Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", 
"MinPoints" -> ppR,"MaxPoints" -> ppR, "DifferenceOrder" -> 2}}

or

Method -> {"Adams", "MaxDifferenceOrder" -> 1}

If someone has an explanation…

[anon_91240867]

The below works for me. As I said above the 4th order equations change very rapidly in time, so it makes sense T=8 (seconds) might be enough to cause approximation errors (or very small vibrations, unmeasurable).

I’m was hoping after my “brief hint” (see above) that you’d see the inequality of y[x,0] and y’[0,t] and adjust them, you did with (Cos[x]-1). Great!

c = 5;
L = 2 Pi;
T = 1/100;

(*EQUATION*)
PDE = D[y[x, t], {t, 2}] + c^2*D[y[x, t], {x, 4}] == 0;

(*INITIAL CONDITIONS*)

IC1 = y[x, 0] == Cos[x] - 1;(*Initial shape of the beam*)
IC2 = D[y[x, t], {t, 1}] == 0 /. {t ->0};(*No initial velocity*)

(*BOUNDARY CONDITIONS*)

(*Clamped on the left side of the beam*)
BCLeft1 = y[x, t] == 0 /. {x -> 0};
BCLeft2 = D[y[x, t], {x, 1}] == 0 /. {x -> 0};
(*Clamped on the right side of the beam*)  
BCRight1 = y[x, t] == 0 /. {x -> L};
BCRight2 = D[y[x, t], {x, 1}] == 0 /. {x -> L};

(*SOLVE*)
soln = 
 NDSolve[{PDE, IC1, IC2, BCLeft1, BCLeft2, BCRight1, BCRight2}, 
   y[x, t], {x, 0, L}, {t, 0, T}][[1, 1, 2]]

Using MaxStepSize, which the Mathematica NDSolve Message suggested, I was able to plot up to T==2. I was unable to do T==8 using MaxStepSize.

     NDSolve[{PDE, IC1, IC2, BCLeft1, BCLeft2, BCRight1, BCRight2}, 
       y[x, t], {x, 0, L}, {t, 0, 2}, MaxStepSize -> .01][[1, 1, 2]]

If instead of MaxStepSize I set c = 5/(10^2) then T==8 is allowed, I think because the solutions on the broken plot are10^32 (and growing).

“lastly” I tried MaxStepFraction → 1/100, and that worked like your second choice did (the softer lobe plot, Method → {“Adams”, “MaxDifferenceOrder” → 1}).

Why? As you know the most basic ODE approx. method has the problem: the further time gets from the initial values provided the worst the error is. For PDE boundary problems the discretization of the boundary effects error and maybe time evolutions does too. I am not familiar with all the methods, nor have I hand approximated this problem, to say more.