3
$\begingroup$

I'm trying to solve Euler-Bernoulli beam equation with simply supported edges.$\frac{\partial^2} {\partial x^2} [ E I \frac{\partial^2 w} {\partial x^2}] + \rho S \frac{\partial^2 w} {\partial t^2} = F_\nu(x,t),$ where $ F_\nu(x,t) = P_f \cos(-\frac{\partial w} {\partial x}) \delta(x-v),$ and $\delta$ is the Dirac delta.With boundary and initial conditions: $w=\frac{\partial^2 w} {\partial x^2}=0, x=0,L $ and $w=\frac{\partial w} {\partial t} =0, t=0$

tau = 1;
L = 2;
Elastic = 1;
Imoment = 1;
rho = 1;
S = 1;
Pf = 0.002;
v = L/20;

a = 10^-4;
del[x_] := 1/(3.14 a) Exp[-(x/a)^2]
Fnu[x_, t_] := Pf Cos[-D[w[x, t], x]] del[x - v]
eqEB1 := D[Elastic*Imoment*D[w[x, t], {x, 2}], {x, 2}] + 
   S*rho*D[w[x, t], {t, 2}] - Fnu[x, t];

Boundary and initial conditions will be :

bc = {w[0, t] == w[L, t] == w[x, 0] == 0, 
  Derivative[2, 0][w][0, t] == Derivative[2, 0][w][L, t] == 
   Derivative[0, 1][w][x, 0] == 0}

When i tried to solve it numerically by using NDSolve , it showed me an error:

solution = 
 NDSolveValue[{D[Elastic*Imoment*D[w[x, t], {x, 2}], {x, 2}] + 
     S*rho*D[w[x, t], {t, 2}] - Fnu[x, t] == 0, 
   w[0, t] == w[L, t] == w[x, 0] == 0, 
   Derivative[2, 0][w][0, t] == Derivative[2, 0][w][L, t] == 
    Derivative[0, 1][w][x, 0] == 0}, {w[x, t]}, {x, 0, L}, {t, 
   0, tau}, Method -> {"FiniteElement"}]

NDSolveValue::femcmsd: The spatial derivative order of the PDE may not exceed two.

I've tried to rewrite it as a system of two second order equations as shown there.And another error occurs:

NDSolve::femnonlinear: Nonlinear coefficients are not supported in this version of NDSolve.

But when I change Fnu, it works just fine. For instance :

Fnu[x_, t_] := Sin[3.14 x] Sin[3.14 t]

Any suggestions will be helpful. Thanks in advance.

$\endgroup$
  • 3
    $\begingroup$ "NDSolve::femnonlinear: Nonlinear coefficients are not supported in this version of NDSolve." This means that Mathematica's FEM facilities cannot treat any PDE for w in which the coefficients depend on w at the moment. Try a different setting for Method, for example, Method -> "BoundaryValues" -> {"Shooting"}. See the section "Method" in the documentation of NDSolve. $\endgroup$ – Henrik Schumacher Feb 2 '19 at 11:47
  • 1
    $\begingroup$ See mathematica.stackexchange.com/questions/187448/… $\endgroup$ – Alex Trounev Feb 2 '19 at 19:42
  • 1
    $\begingroup$ In version 12.0 the rewrite to a system of 2 second order equations and the nonlinear Fnu should work. $\endgroup$ – user21 Apr 17 '19 at 12:15
4
$\begingroup$

We use "MethodOfLines"

tau = 1;
L = 2;
Elastic = 1;
Imoment = 1;
rho = 1;
S = 1;
Pf = 0.002;
v = L/20;

a = 10^-2;
del[x_] := If[x >= 5*a, 0, 1/(Pi a) Exp[-(x/a)^2]]
Fnu[x_, t_] := Pf Cos[-D[w[x, t], x]] del[x - v]
eqEB1 := D[Elastic*Imoment*D[w[x, t], {x, 2}], {x, 2}] + 
   S*rho*D[w[x, t], {t, 2}] - Fnu[x, t];
bc = {w[0, t] == w[L, t] == w[x, 0] == 0, 
   Derivative[2, 0][w][0, t] == Derivative[2, 0][w][L, t] == 
    Derivative[0, 1][w][x, 0] == 0};
sol = NDSolveValue[{D[Elastic*Imoment*D[w[x, t], {x, 2}], {x, 2}] + 
     S*rho*D[w[x, t], {t, 2}] - Fnu[x, t] == 0, 
   w[0, t] == w[L, t] == w[x, 0] == 0, 
   Derivative[2, 0][w][0, t] == Derivative[2, 0][w][L, t] == 
    Derivative[0, 1][w][x, 0] == 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]

Plot3D[sol[x, t], {x, 0, L}, {t, 0, tau}, PlotRange -> All, 
 AxesLabel -> {"x", "t", ""}, Mesh -> None, ColorFunction -> Hue]

fig1

Update 1. If we want to determine the frequencies that are excited, then we must increase tau to 10. Unfortunately this algorithm is unstable at tau = 10, in the end we have a message:

NDSolveValue::eerr: Warning: scaled local spatial error estimate of 16657.48584541172` at t = 10.` in the direction of independent variable x is much greater than the prescribed error tolerance. Grid spacing with 100 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options.

Therefore, we use a different algorithm that allows us to find a steady solution:

AbsoluteTiming[
 sol1 = NDSolveValue[{D[Elastic*Imoment*D[w[x, t], {x, 2}], {x, 2}] + 
       S*rho*D[w[x, t], {t, 2}] - Fnu[x, t] == 0, 
     w[0, t] == w[L, t] == w[x, 0] == 0, 
     Derivative[2, 0][w][0, t] == Derivative[2, 0][w][L, t] == 
      Derivative[0, 1][w][x, 0] == 0}, w, {x, 0, L}, {t, 0, tau}, 
    Method -> {"MethodOfLines", 
      "DifferentiateBoundaryConditions" -> False, 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "MaxPoints" -> 100, "MinPoints" -> 100, 
        "DifferenceOrder" -> 2}}, MaxSteps -> 10^6, 
    EvaluationMonitor :> (currentTime = t;)];]

Numerical solution visualization Figure 2
Now we will check how many points are used or this solution:

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

Map[Length, InterpolatingFunctionCoordinates[sol1]]

Out[]= {100, 26}

These 26 points are not enough to find the frequencies, so we will add an option to increase the number of points

AbsoluteTiming[
 sol2 = NDSolveValue[{D[Elastic*Imoment*D[w[x, t], {x, 2}], {x, 2}] + 
       S*rho*D[w[x, t], {t, 2}] - Fnu[x, t] == 0, 
     w[0, t] == w[L, t] == w[x, 0] == 0, 
     Derivative[2, 0][w][0, t] == Derivative[2, 0][w][L, t] == 
      Derivative[0, 1][w][x, 0] == 0}, w, {x, 0, L}, {t, 0, tau}, 
    Method -> {"MethodOfLines", 
      "DifferentiateBoundaryConditions" -> False, 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "MaxPoints" -> 100, "MinPoints" -> 100, 
        "DifferenceOrder" -> 2}}, MaxStepSize -> 0.05, 
    EvaluationMonitor :> (currentTime = t;)];]

Here we see a periodic solution with a period of 2.5: Figure 3

Now we check number of points

Map[Length, InterpolatingFunctionCoordinates[sol1]]

Out[]= {100, 210}
| improve this answer | |
$\endgroup$
  • $\begingroup$ A great answer. $\endgroup$ – A little mouse on the pampas Apr 1 at 6:00
  • 1
    $\begingroup$ @PleaseCorrectGrammarMistakes It is not a problem to determine frequencies. Will we continue this topic or start a new one? $\endgroup$ – Alex Trounev Apr 3 at 10:25
  • 1
    $\begingroup$ @PleaseCorrectGrammarMistakes See update to my post. $\endgroup$ – Alex Trounev Apr 3 at 13:23
  • 1
    $\begingroup$ @PleaseCorrectGrammarMistakes This is a non-linear problem - the strength depends on Cos[-D[w[x, t], x]]. In this case, the frequencies are determined based on a numerical solution. This can be compared with linear modes. $\endgroup$ – Alex Trounev Apr 4 at 11:29
  • 1
    $\begingroup$ @PleaseCorrectGrammarMistakes There you discuss vibration modes. They gave you an answer. This question about this code can be discussed in another topic. $\endgroup$ – Alex Trounev Apr 6 at 11:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.