Edit only for those interested in large deflections of beams
I discovered a mistake in the equations of the original question (below): in the normal force (compression/traction) n[s_] = EA*u1'[s];
the measure of stretch should not be u1'
. The answers were very instructive in terms of numerical methods, nonetheless.
Actually, I thought the equations of beams under large deflections would be easy to derive as an ODE. I now believe that in the general case, there is no simple explicit ODE to solve; instead the weak form can be projected directly on a mesh.
However, for thin beams undergoing mostly bending (no stretching), the problem can be solved pretty simply (the equations are nicely derived in "Large deflection states of Euler-Bernoulli slender cantilever beam subjected to combined loading" by Žiga Gosar and Franc Kosel, for instance). Numerically speaking, it seems to be much simpler to solve for the rotation field first, and then for the displacement field, rather than seeking the displacement field directly as I did below.
So, for those interested in large deflections of thin beam without stretching, these are some equations you could use:
(* Governing ODE for the rotations *)
eq = theta''[s] == q/EI*s*Cos[theta[s]] - (Q + q*L)/EI*Cos[theta[s]];
thetasol = First@NDSolveValue[{eq}~Join~{theta'[1] == 0, theta[0] == 0}, {theta}
, {s, 0, L}];
(* Computation of the displacement field from the rotations *)
{xsol, ysol} = NDSolveValue[{x'[s] == Cos[thetasol[s]], y'[s] == Sin[thetasol[s]]
, x[0] == 0, y[0] == 0}, {x, y}, {s, 0, L}];
(* Plot of the deformed shape *)
ParametricPlot[{xsol[s], -ysol[s]}, {s, 0, L}, PlotRange -> Full
, AspectRatio -> Automatic]
Now, back to the original question
A lot of questions relate to solving the Euler-Bernoulli beam equation, mostly in dynamics. Actually, they mostly tackle the governing PDEs of the form $$\dfrac{\partial^2 w}{\partial t^2} + \dfrac{\partial^4 w}{\partial w^4}=0$$ which corresponds to a linearized beam equation.
Here, I would like to find the shape of a clamped-free beam (for instance) with large deflection due to gravity; consider a sheet of paper with one clamped edge for example.
This question also tries to address large deflection, but in my case, gravity couples axial and transverse displacements fields, plus I don't have a constraint on length.
So, let's write the equations in the local frame attached to the beam (ft
for the force density in the tangential direction, fn
for the force density in the normal direction):
eqs = {n'[s] - v[s]*kappa[s] + ft[s] == 0, (* local equilibrium, tang. direction *)
v'[s] + n[s]*kappa[s] + fn[s] == 0, (* local equilibrium, transverse direction *)
m'[s] + v[s] == 0} (* local equilibrium, moment *)
The beam, initially straight along the $x$ axis (between $x=0$ and $x=1$) has a deformed shape given by the parametric equation: $$(s+u_1(s), u_2(s))$$
The corresponding curvature and local frame are given by:
{{kappa[s_]}, {tvec[s_], nvec[s_]}} = FrenetSerretSystem[{s + u1[s], u2[s]}, s];
Then, with Euler-Bernoulli kinematics, the internal tangential force field n
and internal bending moment field m
are given by:
EA = EI = 1000;
n[s_] = EA*u1'[s];
m[s_] = EI*kappa[s];
Then, the gravity is projected into the local frame:
gravity = {0, -10};
ft[s_] = gravity.tvec[s]
fn[s_] = gravity.nvec[s]
The third equation in eqs
can be used to eliminate v
:
v[s_] = v[s] /. (Solve[eqs[[3]], v[s]] // Last // Last) // Simplify;
eqs = eqs[[1 ;; 2]] // Simplify;
Along with the following boundary conditions ($u_1(0)=u_2(0)=0$, $u_2'(0)=0$ for the clamped end, $u_1''(1) = u_2''(1) = 0$, $u_1'(1) = 0$):
cls = {u1[0] == 0, u2[0] == 0, u2'[0] == 0, u1''[1] == 0, u2''[1] == 0, u1'[1] == 0}
Finally:
NDSolve[eqs~Join~cls, {u1, u2}, {s, 0, 1}]
returns two successive errors:
NDSolve::ntdvdae: Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations.
NDSolve::bvdae: Differential-algebraic equations must be given as initial value problems.
Any idea on how to solve this system?
Note It might be reasonable to neglect some terms (especially the squares of first derivatives) but NDSolve
returns the same error.
xzczd suggested using his function pdftoae
but I did not manage to make it work for my system of ODEs.
NDSolve
can't put the system in a first order form (Solve
method), withResidual
it can't solve BVPs, and it fails withMassMatrix
. Should I understand that this can't be solved with MMA? I'm quite surprised (at least for the last one, which is rather simple). $\endgroup$FiniteElement
method may be a choice. A more straightforward (at least for me) approach is FDM, this can be easily implemented with mypdetoae
. $\endgroup$EA = EI = 100
(for the realistic picture). Alsopdetoae
is a different tool thenpdetoode
I am used in my answer. $\endgroup$u1''[1] == 0
is duplicated incls
. What is the correct expression forcls
? $\endgroup$