# Euler-Bernoulli beam equation

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.

• "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. – Henrik Schumacher Feb 2 at 11:47
• – Alex Trounev Feb 2 at 19:42
• In version 12.0 the rewrite to a system of 2 second order equations and the nonlinear Fnu should work. – user21 Apr 17 at 12:15

## 1 Answer

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] 