# Geometrically nonlinear beam deflection

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' == 0, theta == 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, y == 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[], 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, u2 == 0, u2' == 0, u1'' == 0, u2'' == 0, u1' == 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.

• 1. It's better to include the system in this post, external link may be broken in future. 2. Have you read this?: mathematica.stackexchange.com/a/158519/1871 – xzczd Sep 26 at 12:35
• @xzczd 1. Done. 2. Not very reassuring: NDSolve can't put the system in a first order form (Solve method), with Residual it can't solve BVPs, and it fails with MassMatrix. 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). – anderstood Sep 26 at 13:09
• If you can rewrite the system as 2nd order ODE system, then the FiniteElement method may be a choice. A more straightforward (at least for me) approach is FDM, this can be easily implemented with my pdetoae. – xzczd Sep 26 at 13:23
• @anderstood For gravity force of 10 it is better to use EA = EI = 100 (for the realistic picture). Also pdetoae is a different tool then pdetoode I am used in my answer. – Alex Trounev Sep 26 at 23:24
• u1'' == 0 is duplicated in cls. What is the correct expression for cls? – bbgodfrey Sep 27 at 0:14

Actually this system can be solved with NDSolve with some efforts. We use 3 equation:

eqs = {n'[s] - v[s]*kappa[s] + ft[s] == 0,
v'[s] + n[s]*kappa[s] + fn[s] == 0,m'[s] + v[s] == 0};
{{kappa[s_]}, {tvec[s_], nvec[s_]}} =
FrenetSerretSystem[{s + u1[s], u2[s]}, s]; EA = 1000; EI = 1000;
n[s_] = EA*u1'[s];
m[s_] = EI*kappa[s]; gravity = {0, -10};
ft[s_] = gravity.tvec[s];
fn[s_] = gravity.nvec[s];


Now define function dependent on 3 parameters

solp[x_?NumericQ, y_?NumericQ, z_?NumericQ] :=
Module[{p1 = x, p2 = y, p3 = z},
sol = NDSolve[
Flatten[{eqs, {u1 == 0, u2 == 0, u2' == 0, u1' == p1,
u2'' == p2, v == p3}}], {u1, u2, v}, {s, 0, 1},
Method -> {"EquationSimplification" -> "Residual"}]; sol[]];


With this function we calculate initial data at s=1

U1[x_?NumericQ, y_?NumericQ, z_?NumericQ] :=
u1'' /. solp[x, y, z];
U2[x_?NumericQ, y_?NumericQ, z_?NumericQ] := u2'' /. solp[x, y, z];
U3[x_?NumericQ, y_?NumericQ, z_?NumericQ] := u1' /. solp[x, y, z]

init = {u1'', u2'', u1'} /. solp[0, 0, 0];

solf =
FindRoot[{U1[x, y, z] == 0, U2[x, y, z] == 0,
U3[x, y, z] == 0}, {{x, init[]}, {y, init[]}, {z, init[]}}]

(*Out[]= {x -> -7.52634*10^-10, y -> -0.00166661, z -> -6.66661}*)


Finally we visualize solution and compare with pdetoae solution

{Plot[Evaluate[u1[s] /. (solp[x, y, z] /. solf)], {s, 0, 1},
AxesLabel -> {"s", "u1"}],
Plot[Evaluate[u2[s] /. (solp[x, y, z] /. solf)], {s, 0, 1},
AxesLabel -> {"s", "u2"}],
Plot[Evaluate[v[s] /. (solp[x, y, z] /. solf)], {s, 0, 1},
AxesLabel -> {"s", "v"}]} • In v12.1.1 the solf is {x -> -7.52634*10^-10, y -> -0.00166661, z -> -6.66661}, with these parameters I obtain the same solution as shown in my answer. – xzczd Sep 27 at 14:37
• @xzczd Ah, sorry, it is my parameters were differ. Now corrected. Thank you! – Alex Trounev Sep 27 at 14:51
• It's interesting that retaining 3rd equation seems to be critical here. – xzczd Sep 29 at 3:06
• @xzczd Yes, it leads to less order system of equations and it is easy with NDSolve to get right solution. – Alex Trounev Sep 29 at 10:29
• @AlexTrounev No your answer addresses the OP. I think you do not need to change anything. But I wanted to add something to my question so that people who look for "large deflections of beam" do not find my wrong equations :) – anderstood Oct 15 at 11:24

Let me add a solution based on finite difference method (FDM). I'll use pdetoae for the generation of difference equations.

domain = {0, 1}; points = 50; difforder = 8;
grid = Array[# &, points, domain];
(* Definition of pdetoae isn't included in this post,
ptoafunc = pdetoae[{u1, u2}[s], grid, difforder];
ae1 = ptoafunc@eqs[] // Delete[#, {{1}, {2}, {-1}}] &;
ae2 = ptoafunc@eqs[] // Delete[#, {{1}, {-2}, {-1}}] &;
aebc = cls // ptoafunc;
guess[_, x_] := 0
sollst = Partition[#, points] &@
FindRoot[{ae1, ae2, aebc} // Flatten,
Table[{var[x], guess[var, x]}, {var, {u1, u2}}, {x, grid}] //
Flatten[#, 1] &][[All, -1]];

solfunclst = ListInterpolation[#, grid, InterpolationOrder -> difforder] & /@ sollst

ListLinePlot /@ sollst Error check:

Subtract @@@ cls /. Thread[{u1, u2} -> solfunclst]
(* {2.06795*10^-23, 5.29396*10^-23, 9.7917*10^-19,
-7.22304*10^-15, -7.42942*10^-15, -1.96557*10^-17} *)

• What equations and parameters do you use in your code? – Alex Trounev Sep 27 at 13:53
• @alex Those given in the question, of course. – xzczd Sep 27 at 14:09
• Please, check my answer. The curves look similar but different in scale with yours. May be you take gravity as {0,-1}? – Alex Trounev Sep 27 at 14:14
$Version (* "12.1.1 for Microsoft Windows (64-bit) (June 19, 2020)" *)  NDSolve "cannot solve to find an explicit formula for the derivatives", because only one of the two ODEs is fourth order, as can be seen by determining the positions of {u1''''[s], u2''''[s]}. Position[eqs, u1''''[s]] (* {{2, 1, 3, 4, 3, 1, 3}, {2, 1, 3, 4, 3, 2, 2, 3, 3, 2, 2}} *) Position[eqs, u2''''[s]] (* {{2, 1, 3, 4, 3, 3, 3, 2, 2}, {2, 1, 3, 4, 3, 4, 3, 4, 2}} *)  Indeed, there are no fourth derivatives in eqs[]. Nonetheless, some progress can be made. For convenience, define eq1 = Subtract @@ (eqs[]); eq2 = Subtract @@ (eqs[]);  which moves all terms to the left side of the equations and then discards == 0. Next, obtain the highest order derivatives in each expression. eq1h = Simplify[Collect[eq1, {u1'''[s], u2'''[s]}, Simplify][[-2 ;; -1]]] (* ((u2'[s]*u1''[s] - (1 + u1'[s])*u2''[s])*(u2'[s]*u1'''[s] - (1 + u1'[s])*u2'''[s])) /(1 + 2*u1'[s] + u1'[s]^2 + u2'[s]^2)^3 *) eq2h = Simplify[Collect[eq2, {u1''''[s], u2''''[s]}, Simplify][[-2 ;; -1]]] (* (u2'[s]*u1''''[s] - (1 + u1'[s])*u2''''[s]) /(1 + 2*u1'[s] + u1'[s]^2 + u2'[s]^2)^(3/2) *)  The similarity of these two terms indicates that the fourth derivatives can be eliminated from eq2, as follows. rat = Simplify[eq2h/eq1h (u2'[s] u1'''[s] - (1 + u1'[s]) u2'''[s])/ (u2'[s] u1''''[s] - (1 + u1'[s]) u2''''[s])] (* (1 + 2*u1'[s] + u1'[s]^2 + u2'[s]^2)^(3/2)/ (u2'[s]*u1''[s] - (1 + u1'[s])*u2''[s]) *) eq21 = Collect[eq2 - D[rat*eq1, s], {u1''''[s], u2''''[s]}, Simplify];  Although the resulting expression for eq21 is too long to reproduce here, inspecting it using {Coefficient[eq21, u1''''[s]], Coefficient[eq21, u2''''[s]]} (* {0, 0} *)  verifies that the the fourth derivatives terms indeed are gone. Moreover, Flatten@Solve[{eq1 == 0, eq21 == 0}, {u1'''[s], u2'''[s]}]  gives explicit expressions for {u1'''[s], u2'''[s]}. So, NDSolve can in principle integrate {eq1 == 0, eq21 == 0}. To do so requires specifying six boundary conditions. Presumably, {u1''' == 0, u2''' == 0} should be dropped from cls. In addition, as noted in my comment, u1'' == 0 is duplicated in cls. Let us assume that the OP meant one of the duplicates to be u2'' == 0. With these changes, cls = {u1 == 0, u2 == 0, u2' == 0, u1' == 0, u2'' == 0, u1'' == 0}  At this point, NDSolve[{eq1 == 0, eq21 == 0, cls}, {u1[s], u2[s]}, {s, 0, 1}]  runs for a while without error but eventually crashes as it searches for a boundary value solution. Having a rough guess for the solution probably would yield an exact solution. • Mmh. On my MMA 12.1.1, eq21 still has fourth order derivatives: Position[eq21, u1''''[s]] return {{2, 4}} for instance. And the NDSolve returns the same answer as in the OP, probably for that reason. I'll try to find what's wrong in the elimination of the 4th order terms. – anderstood Sep 27 at 7:12 • Found a workaround: eq21 = Collect[ eq2 /. Solve[D[eq1, s] == 0, u1''''[s]], {u1''''[s], u2''''[s]}, Simplify];. Now eq21 is a 3d-order ODE. I'll now investigate NDSolve which no longer returns the previous error. Thanks for the nice observation. – anderstood Sep 27 at 12:19 • My outputs agree with yours up to rat, which is -(1 + 2 Derivative[u1][s] + Derivative[u1][s]^2 + Derivative[u2][s]^2)^3 on my notebook. Maybe you multiplied the equations by the denominator? I get eqs by running the code from my answer. With EA=EI=1 (which I changed in the OP), the equations should be: eqs = Import["https://pastebin.com/raw/hm1tcE6X"]. Anyway, the point is that the fourth-order terms can be eliminated, now I have to wait for NDSolve to output something. – anderstood Sep 27 at 12:29 • @anderstood My expression for rat is precisely equal to eq2h[]/eq1h[[1 ;; 2]]. Perhaps, the order of terms in eq1h and eq2h differs between us. I am using $Version 12.1.1 for Microsoft Windows (64-bit) (June 19, 2020). – bbgodfrey Sep 27 at 12:48