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}$

Cantilever beam with exponential loading



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():
    styles = open('./styles/custom_barba.css', 'r').read()
    return HTML(styles)
css_styling()









    Out[36]:



In [ ]: