CIVL4250 Numerical Methods in Engineering

Dr Dorival Pedroso

The code for this notebook is available in https://github.com/cpmech/CIVL4250py

Euler-Bernoulli beam theory – with constant properties

The following fourth order partial differential equation

$$ \frac{\mathrm{d}^2}{\mathrm{d} x^2} \left (E\,I\,\frac{\mathrm{d}^2 u}{\mathrm{d} x^2} \right) = q(x) $$

represents the transverse deflection ($u$) of a beam according to Euler-Bernoulli beam theory (see, e.g. Euler-Bernoulli in Wikipedia. Therein, $q(x)$ is a transverse distributed load, $E$ is the Young's elasticity modulus and $I$ the moment of inertia of the cross-sectional area.

With constant properties ($E\,I$), simple loads and simple boundary conditions, the above equation can be solved in closed form. For the cantilever illustrated above (fully fixed at $x=0$, with a concentrated load $F$ at $x=L$ and distributed load given by $q(x)=x^2$), the resulting transverse displacement $u$ can be analytically deduced.

The following convention is adopted: positive values of displacements and forces are those corresponding to the $y$ axis. In this way, negative values of q and F mean downward loads.

Finite differences approximation

With constant properties, the differential equation becomes

$$ \frac{\mathrm{d}^4 u}{\mathrm{d} x^4} = f(x) $$

where $f(x) = q(x) / (E\,I)$.

This equation can be approximated with a 5-point stencil (grid) using the following fourth order central differences formula

$$ \frac{\mathrm{d}^4 u}{\mathrm{d} x^4} \approx \frac{ u_{i-2} - 4\,u_{i-1} + 6\,u_{i} - 4\,u_{i+1} + u_{i+2} } {h^4} $$

where $h = \Delta x$ and the grid indices go from $0$ to $n-1$ with $n$ being the number of points.

Substituting the above approximation into the differential equation, the following general expression for any point is obtained

$$ u_{i-2} - 4\,u_{i-1} + 6\,u_{i} - 4\,u_{i+1} + u_{i+2} = f_i \, h^4 $$

where the points outside $0\leq i\leq n-1$ are known as ghost points and must be handled after considering the boundary conditions.

Boundary conditions (cantilever)

The boundary conditions are now taken into consideration.

At $x=0$, the deflection $u_0$ at point $i=0$ is known and is equal to zero. The derivative at $i=0$ is also zero; thus, by means of a first order finite differences approximation at $i=0$,

$$ u_0 = 0 \qquad \text{and} \qquad u_{-1} = u_{1} $$

At $x=L$, although the deflection is maximum, its second derivative is zero because the bending moment is zero at this position. Note that the bending moment is proportional to the second derivative of $u$. Thus, by using a second order central differences formula

$$ \left.{\frac{\mathrm{d}^2 u}{\mathrm{d} x^2}}\right\rvert_{x=L} \approx \frac{ u_{b-1} - 2\,u_b + u_{b+1} } {h^2} = 0 $$

where $b$ is the index of the last point at $x=L$; i.e. $b=n-1$. Therefore

$$ u_{b+1} = 2\,u_b - u_{b-1} $$

The third derivative of deflection is proportional to the shear force $Q$, i.e.:

$$ Q = -E\,I\,\frac{\mathrm{d}^3 u}{\mathrm{d} x^3} $$

At $x=L$, the shear force equals the applied concentrated force $F$ (that could be zero). Therefore, with a third order central differences,

$$ \left.{\frac{\mathrm{d}^3 u}{\mathrm{d} x^3}}\right\rvert_{x=L} \approx \frac{ -0.5\,u_{b-2} + u_{b-1} - u_{b+1} + 0.5\,u_{b+2} } {h^3} = g_b $$

where $g_b = -F / (E\,I)$.

From the above approximation, the deflection at the right-most (ghost) point $b+2$ can be found as follows

$$ u_{b+2} = 2\,g_b\,h^3 + u_{b-2} - 2\,u_{b-1} + 2\,u_{b+1} $$

which, after substituting the results for $u_{b+1}$, becomes

$$ u_{b+2} = 2\,g_b\,h^3 + u_{b-2} - 4\,u_{b-1} + 4\,u_b $$

Linear system of equations

Now, the expressions for the ghost points can be inserted into the general formula for any point. Note also that the equation for point $i=0$ is not required since $u_0=0$.

For point $i=1$, from the general expression given in the previous section, the equation becomes

$$ u_{-1} - 4\,u_{0} + 6\,u_{1} - 4\,u_{2} + u_{3} = f_1 \, h^4 $$

thus, since $u_0=0$ and $u_{-1}=u_1$,

$$ \text{eq. for point i=1}:\quad 7\,u_{1} - 4\,u_{2} + u_{3} = f_1 \, h^4 $$

Similarly, at the right-hand side, for point $i=b-1$,

$$ u_{b-1-2} - 4\,u_{b-1-1} + 6\,u_{b-1} - 4\,u_{b-1+1} + u_{b-1+2} = f_{b-1} \, h^4 $$

i.e.

$$ u_{b-3} - 4\,u_{b-2} + 6\,u_{b-1} - 4\,u_b + u_{b+1} = f_{b-1} \, h^4 $$

which, after substituting the value for ghost point $b+1$, becomes

$$ \text{eq. for point i=b-1}:\quad u_{b-3} - 4\,u_{b-2} + 5\,u_{b-1} - 2\,u_b = f_{b-1} \, h^4 $$

Finally, for point $i=b$,

$$ u_{b-2} - 4\,u_{b-1} + 6\,u_b - 4\,u_{b+1} + u_{b+2} = f_b \, h^4 $$

which, after substituting the value for ghost points $b+1$ and $b+2$, becomes

$$ u_{b-2} - 4\,u_{b-1} + 6\,u_b - 4\,(2\,u_b - u_{b-1}) + (2\,g_b\,h^3 + u_{b-2} - 4\,u_{b-1} + 4\,u_b) = f_b \, h^4 $$

i.e.

$$ 2\,u_{b-2} - 4\,u_{b-1} + 2\,u_b = f_b \, h^4 - 2\,g_b\,h^3 $$

or, dividing by 2,

$$ \text{eq. for point i=b}:\quad u_{b-2} - 2\,u_{b-1} + u_b = 0.5\,f_b \, h^4 - g_b\,h^3 $$

Hence, the resulting linear system of equations, expressed in matrix form, is

$$ \begin{bmatrix} 7 & -4 & 1 & & & & & & & \\ -4 & 6 & -4 & 1 & & & & & & \\ 1 & -4 & 6 & -4 & 1 & & & & & \\ & 1 & -4 & 6 & -4 & \ddots & & & & \\ & & 1 & -4 & 6 & \ddots & 1 & & & \\ & & & 1 & -4 & \ddots & -4 & 1 & & \\ & & & & 1 & \ddots & 6 & -4 & 1 & \\ & & & & & \ddots & -4 & 6 & -4 & 1 \\ & & & & & & 1 & -4 & 5 & -2 \\ & & & & & & & 1 & -2 & 1 \end{bmatrix} \begin{Bmatrix} u_1 \\ u_2 \\ u_3 \\ u_4 \\ u_i \\ \vdots \\ u_{b-3} \\ u_{b-2} \\ u_{b-1} \\ u_b \end{Bmatrix} = \begin{Bmatrix} f_1 \, h^4 \\ f_2 \, h^4 \\ f_3 \, h^4 \\ f_4 \, h^4 \\ f_i \, h^4 \\ \vdots \\ f_{b-3} \, h^4 \\ f_{b-2} \, h^4 \\ f_{b-1} \, h^4 \\ 0.5\,f_b \, h^4 - g_b\,h^3 \end{Bmatrix} $$

where we can observe that the resulting matrix is symmetric and banded. Note also that this matrix corresponds to the $K_{11}$ matrix from the lecture notes, since the prescribed boundary conditions were already considered.

Calculation script

The Python calculation script for the cantilever beam equations is similar to the one presented for the 1D Poisson equation in the lecture notes. The complete code can be downloaded from here.

Considering the figure below

where negative values of q and F mean downward loads, the input data and FDM grid are defined as follows

Although a specialised solver considering the banded and symmetric form of the $K$ matrix could provide the fastest computations, a general sparse solver as used for the finite element code is employed here. We start then with the assemblage of the $K_{11}$ matrix using the lil_matrix structure from scipy.

Now, the system of equations can be solved for the deflections according to

and the results can be plotted as follows