| Home |

Mathematics and fly fishing

Numerical solution to the Euler-Bernouilli equation

The Euler-Bernouilli equation forms the point of departure for all considerations in the technical theory of beams.

K.E.Bisshopp and D.C.Drucker have found a solution for large deflection of a cantilever beam with constant modulus of elasticity and moment of inertia for a concentrated vertical load at the free end. The solution is limited to calculating the coordinates and the slope for the end point.

The simple numerical solution described below is a very powerful tool. It solves the problem of finding the deflection curve for a cantilever beam with variable modulus of elasticity, variable moments of inertia and with no restrictions as to have the moment is applied. A problem that up till now have required a Finite Element Analysis.

The Euler-Bernouilli equation is usually presented as:

(1)

M = moment
E = modulus of elasticity
I = moment of inertia

In Equation (1) set dy/dx= tan q to obtain

For a small but final interval Dx we have:

Now we need a relation between Dq and R.

Assuming that the curve, when looked upon at small intervals, can be considered made up of n arcs of circles and L being the length of the beam, then we have:

And the solution must have the form:

for i = 1 to n
q_i must be a value inside the interval [i-1, i].
This looks very simple, but the snag is that we do not know q_i.

The idea to the numerical solution came from the observation that the Euler-Bernouilli equation can be written as:

The left-hand side of the equation is the curvature K of the deflection curve.

So if the right hand side of the equation is known then the solution of the complicated differential equation is transformed into a simple geometrical problem.

But the seemingly tough problem is that the moments are changing as the beam deflects and we are not able to calculate the moments unless we know the deflection curve. Using the successive approximation method overcomes that hurdle.

Nodes are placed at the centroid of the beam, dividing the beam into elements.
A curve defined by its curvature is not constrained by a coordinate system. We may place the beam in a Cartesian coordinate system with the clamped end located at the origin and the free end pointing to the right. The nodes are numbered from the clamped end starting with 1.

We impose the condition, that the tangent to the deflection curve shall be parallel to the axis of abscissas at the origin.

When the beam deflects the centerline of each element is replaced by an arc of a circle.

The first arc starts at the origin of the coordinate system, node i = 1 and ends at node i = 2. The first circle has its center on the y-axis at y = - R₁.

Choosing R₁=1/K₁ we get a deflection that is too small, choosing R₂=1/K₂ we get a deflection that is too large.
The radius that gives the best fit is somewhere between, and we choose the mean value of K to calculate the radius:

To calculate the location of point 2, the end point of the first arc we notice that the length of the arc = L/n equals the radius times the top angle.

The top angle for the arc of the circle from node i to i+1 is denoted a_i+1. It is equal to the increment of the slope of the tangent from node i to node i+1.
a₁is not defined, but it simplifies the programming if we let a₁ = 0.

The slope of tangent at point i is the summation of a from a₁ to a_i .

The q_i thus found is calculated at the node points and therefore not identical with the q_i in the equations found by the more formal analysis.

Knowing the top angle the chord for the circle segment from node i to node i+1 can be calculated

chord_i = 2*R_i*sin(a_i+1/2)

The slope of the chord between two nodes is the average of the slope of the tangents at the end points. A consequence of the curve being an arc of a circle.

b_i = (q_i + q_i+1)/2

Which also can be expressed as: the slope of the chord equals the slope of the curve at the first point plus half the angle of the tangent's increment between the points.

b_i = q_i + a_i+1/2

Now we are ready to start the calculation of the deflection curve:

x₁= 0 ; x_i= x_i-1 + chord_i-1 * cos(b_i-1)

y₁= 0 ; y_i= y_i-1 + chord_i-1 * sin(b_i-1)

As a first approximation we calculate the curve using the moments that we would have in the undeflected beam.

The curve is not the deflection curve we are seeking.
But using the coordinates found we can calculate new moments. The calculations continue until the absolute value of the difference between e.g. K_n-1 in two consecutive runs is less than say 1E-05.

A proof that the iteration will converge is not found.

The solution works for deflection up to 30% of the beam length.

For large deflections the iterations may not be able to achieve a difference of
1E-05 while for small deflections it possible to achieve a difference of 1E-10.

When programming the calculation it is convenient to have the program monitoring the number of runs and the difference between consecutive K_n-1

The geometrical approach results in using the chord while the analytic solution uses the arc of the circle element.

Lets look at the ratio between the chord and the arc for a circle sector for q ® 0.

So for small angles (R large compared to the length of the element) the difference is insignificant.

The accuracy of the calculation.

The successful generation of the original coordinates for a curve from the curvature and the distance between the points is backbone of the solution.

To determine the accuracy of the solution a simple parabola was chosen.

The length of the parabola from (0, 0) to (1200, 1698) is 2159 mm = 85 inches.

The coordinates for points located 31.75 mm apart measured on the curve were calculated.
Then the curvature at the same points were calculated.

When calculating the curve from the curvatures no iteration is needed as we have the correct curvatures. Best result is obtained when staring the calculation from the end with the smallest curvature.
In order to compare the original and the calculated curve we have to move and rotate one of the curves.

The results are shown in the diagram 1, 2 and 3 for 18, 35 and 69 points.
The error decreases with the number of nodes.

Diagram 1

Diagram 2

Diagram 3

For constant E and I, the deflection of the free end matches the results from Bisshopp and Drucker with an error of less than 0.05 %.

Further applications:
It is possible to calculate a beam with a specified stress distribution, fx a bamboo fly rod, by using the difference between the calculated stress and desired stress to modify the moments of inertia.
Variation in the modulus of elasticity and the moments of inertia along the beam do not complicate the calculation.