# Euler Bernoulli Beam "solver"

The Euler-Bernoulli equation describes the relationship between the beam's deflection and the applied load

$$\frac{d^2}{dx^2}\left(EI\frac{d^2w}{dx^2}\right) = q \enspace .$$

The curve $w(x)$ describes the delection of the beam at some point $x$, $q$ is a distributed load.

This equation cannot be solve in this form in Sympy. Nevertheless, we can "trick" it to do it for us. Let us rewrite the equation as two equations

\begin{align} &-\frac{d^2 M}{dx^2} = q \enspace ,\\ &- \frac{d^2w}{dx^2} = \frac{M}{EI} \enspace , \end{align}

where $M$ is the bending moment in the beam. We can, then, solve the two equation as if they have source terms and then couple the two solutions.

The code below do this

In [1]:
from sympy import *

In [2]:
%matplotlib notebook
init_printing()

In [3]:
x = symbols('x')
E, I = symbols('E I', positive=True)
C1, C2, C3, C4 = symbols('C1 C2 C3 C4')
w, M, q, f = symbols('w M q f', cls=Function)
EI = symbols('EI', cls=Function, nonnegative=True)

In [4]:
M_eq = -diff(M(x), x, 2) - q(x)

M_eq

Out[4]:
$\displaystyle - q{\left(x \right)} - \frac{d^{2}}{d x^{2}} M{\left(x \right)}$

In [5]:
M_sol = dsolve(M_eq, M(x)).rhs.subs([(C1, C3), (C2, C4)])

M_sol

Out[5]:
$\displaystyle C_{3} + \int \left(C_{4} - \int q{\left(x \right)}\, dx\right)\, dx$

In [6]:
w_eq = f(x) + diff(w(x),x,2)
w_eq

Out[6]:
$\displaystyle f{\left(x \right)} + \frac{d^{2}}{d x^{2}} w{\left(x \right)}$

In [7]:
w_sol = dsolve(w_eq, w(x)).subs(f(x), M_sol/EI(x)).rhs

w_sol

Out[7]:
$\displaystyle C_{1} + \int \left(C_{2} - \int \frac{C_{3} + \int \left(C_{4} - \int q{\left(x \right)}\, dx\right)\, dx}{\operatorname{EI}{\left(x \right)}}\, dx\right)\, dx$

We want to be sure that this solution is ok. We replaced known values for $E$, $I$ and $q$ to check it.

## Cantilever beam with end load

In [8]:
sub_list = [(q(x), 0), (EI(x), E*I)]
w_sol1 = w_sol.subs(sub_list).doit()

In [9]:
L, F = symbols('L F')
# Fixed end
bc_eq1 = w_sol1.subs(x, 0)
bc_eq2 = diff(w_sol1, x).subs(x, 0)
# Free end
bc_eq3 = diff(w_sol1, x, 2).subs(x, L)
bc_eq4 = diff(w_sol1, x, 3).subs(x, L) + F/(E*I)

In [10]:
[bc_eq1, bc_eq2, bc_eq3, bc_eq4]

Out[10]:
$\displaystyle \left[ C_{1}, \ C_{2}, \ - \frac{C_{3} + C_{4} L}{E I}, \ - \frac{C_{4}}{E I} + \frac{F}{E I}\right]$

In [11]:
constants = solve([bc_eq1, bc_eq2, bc_eq3, bc_eq4], [C1, C2, C3, C4])
constants

Out[11]:
$\displaystyle \left\{ C_{1} : 0, \ C_{2} : 0, \ C_{3} : - F L, \ C_{4} : F\right\}$

In [12]:
w_sol1.subs(constants).simplify()

Out[12]:
$\displaystyle \frac{F x^{2} \left(3 L - x\right)}{6 E I}$

## Cantilever beam with uniformly distributed load

In [13]:
sub_list = [(q(x), 1), (EI(x), E*I)]
w_sol1 = w_sol.subs(sub_list).doit()

In [14]:
L = symbols('L')
# Fixed end
bc_eq1 = w_sol1.subs(x, 0)
bc_eq2 = diff(w_sol1, x).subs(x, 0)
# Free end
bc_eq3 = diff(w_sol1, x, 2).subs(x, L)
bc_eq4 = diff(w_sol1, x, 3).subs(x, L)

In [15]:
constants = solve([bc_eq1, bc_eq2, bc_eq3, bc_eq4], [C1, C2, C3, C4])

In [16]:
w_sol1.subs(constants).simplify()

Out[16]:
$\displaystyle \frac{x^{2} \left(6 L^{2} - 4 L x + x^{2}\right)}{24 E I}$

In [17]:
sub_list = [(q(x), exp(x)), (EI(x), E*I)]
w_sol1 = w_sol.subs(sub_list).doit()

In [18]:
L = symbols('L')
# Fixed end
bc_eq1 = w_sol1.subs(x, 0)
bc_eq2 = diff(w_sol1, x).subs(x, 0)
# Free end
bc_eq3 = diff(w_sol1, x, 2).subs(x, L)
bc_eq4 = diff(w_sol1, x, 3).subs(x, L)

In [19]:
constants = solve([bc_eq1, bc_eq2, bc_eq3, bc_eq4], [C1, C2, C3, C4])

In [20]:
w_sol1.subs(constants).simplify()

Out[20]:
$\displaystyle \frac{- \frac{x^{3} e^{L}}{6} + \frac{x^{2} \left(L - 1\right) e^{L}}{2} - x + e^{x} - 1}{E I}$

## Load written as a Taylor series and constant EI

We can prove that the general function is written as

In [21]:
k = symbols('k', integer=True)
C = symbols('C1:4')
D = symbols('D', cls=Function)

In [22]:
w_sol1 = 6*(C1 + C2*x) - 1/(E*I)*(3*C3*x**2 + C4*x**3 -
6*Sum(D(k)*x**(k + 4)/((k + 1)*(k + 2)*(k + 3)*(k + 4)),(k, 0, oo)))

w_sol1

Out[22]:
$\displaystyle 6 C_{1} + 6 C_{2} x - \frac{3 C_{3} x^{2} + C_{4} x^{3} - 6 \sum_{k=0}^{\infty} \frac{x^{k + 4} D{\left(k \right)}}{\left(k + 1\right) \left(k + 2\right) \left(k + 3\right) \left(k + 4\right)}}{E I}$

## Uniform load and varying cross-section

In [23]:
Q, alpha = symbols("Q alpha")
sub_list = [(q(x), Q), (EI(x), E*x**3/12/tan(alpha))]
w_sol1 = w_sol.subs(sub_list).doit()

In [24]:
M_eq = -diff(M(x), x, 2) - Q

M_eq

Out[24]:
$\displaystyle - Q - \frac{d^{2}}{d x^{2}} M{\left(x \right)}$

In [25]:
M_sol = dsolve(M_eq, M(x)).rhs.subs([(C1, C3), (C2, C4)])

M_sol

Out[25]:
$\displaystyle C_{3} + C_{4} x - \frac{Q x^{2}}{2}$

In [26]:
w_eq = f(x) + diff(w(x),x,2)
w_eq

Out[26]:
$\displaystyle f{\left(x \right)} + \frac{d^{2}}{d x^{2}} w{\left(x \right)}$

In [27]:
w_sol1 = dsolve(w_eq, w(x)).subs(f(x), M_sol/(E*x**3/tan(alpha)**3)).rhs

w_sol1 = w_sol1.doit()

In [28]:
expand(w_sol1)

Out[28]:
$\displaystyle C_{1} + C_{2} x - \frac{C_{3} \tan^{3}{\left(\alpha \right)}}{2 E x} + \frac{C_{4} \log{\left(x \right)} \tan^{3}{\left(\alpha \right)}}{E} + \frac{Q x \log{\left(x \right)} \tan^{3}{\left(\alpha \right)}}{2 E} - \frac{Q x \tan^{3}{\left(\alpha \right)}}{2 E}$

In [29]:
limit(w_sol1, x, 0)

Out[29]:
$\displaystyle - \infty \operatorname{sign}{\left(C_{3} \tan^{3}{\left(\alpha \right)} \right)}$

In [30]:
L = symbols('L')
# Fixed end
bc_eq1 = w_sol1.subs(x, L)
bc_eq2 = diff(w_sol1, x).subs(x, L)
# Finite solution
bc_eq3 = C3

In [31]:
constants = solve([bc_eq1, bc_eq2, bc_eq3], [C1, C2, C3, C4])

In [32]:
simplify(w_sol1.subs(constants).subs(C4, 0))

Out[32]:
$\displaystyle \frac{Q \left(L - x \left(\log{\left(L \right)} + 1\right) + x \log{\left(x \right)}\right) \tan^{3}{\left(\alpha \right)}}{2 E}$

The shear stress would be

In [33]:
M = -E*x**3/tan(alpha)**3*diff(w_sol1.subs(constants).subs(C4, 0), x, 2)
M

Out[33]:
$\displaystyle - \frac{Q x^{2}}{2}$

In [34]:
diff(M, x)

Out[34]:
$\displaystyle - Q x$

In [35]:
w_plot = w_sol1.subs(constants).subs({C4: 0, L: 1, Q: -1, E: 1, alpha: pi/9})
plot(w_plot, (x, 1e-6, 1));

In [ ]:

In [36]:
from IPython.core.display import HTML
def css_styling():