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
from sympy import *
%matplotlib notebook
init_printing()
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)
M_eq = -diff(M(x), x, 2) - q(x)
M_eq
M_sol = dsolve(M_eq, M(x)).rhs.subs([(C1, C3), (C2, C4)])
M_sol
w_eq = f(x) + diff(w(x),x,2)
w_eq
w_sol = dsolve(w_eq, w(x)).subs(f(x), M_sol/EI(x)).rhs
w_sol
We want to be sure that this solution is ok. We replaced known values for $E$, $I$ and $q$ to check it.
sub_list = [(q(x), 0), (EI(x), E*I)]
w_sol1 = w_sol.subs(sub_list).doit()
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)
[bc_eq1, bc_eq2, bc_eq3, bc_eq4]
constants = solve([bc_eq1, bc_eq2, bc_eq3, bc_eq4], [C1, C2, C3, C4])
constants
w_sol1.subs(constants).simplify()
sub_list = [(q(x), 1), (EI(x), E*I)]
w_sol1 = w_sol.subs(sub_list).doit()
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)
constants = solve([bc_eq1, bc_eq2, bc_eq3, bc_eq4], [C1, C2, C3, C4])
w_sol1.subs(constants).simplify()
sub_list = [(q(x), exp(x)), (EI(x), E*I)]
w_sol1 = w_sol.subs(sub_list).doit()
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)
constants = solve([bc_eq1, bc_eq2, bc_eq3, bc_eq4], [C1, C2, C3, C4])
w_sol1.subs(constants).simplify()
We can prove that the general function is written as
k = symbols('k', integer=True)
C = symbols('C1:4')
D = symbols('D', cls=Function)
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
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()
M_eq = -diff(M(x), x, 2) - Q
M_eq
M_sol = dsolve(M_eq, M(x)).rhs.subs([(C1, C3), (C2, C4)])
M_sol
w_eq = f(x) + diff(w(x),x,2)
w_eq
w_sol1 = dsolve(w_eq, w(x)).subs(f(x), M_sol/(E*x**3/tan(alpha)**3)).rhs
w_sol1 = w_sol1.doit()
expand(w_sol1)
limit(w_sol1, x, 0)
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
constants = solve([bc_eq1, bc_eq2, bc_eq3], [C1, C2, C3, C4])
simplify(w_sol1.subs(constants).subs(C4, 0))
The shear stress would be
M = -E*x**3/tan(alpha)**3*diff(w_sol1.subs(constants).subs(C4, 0), x, 2)
M
diff(M, x)
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));
from IPython.core.display import HTML
def css_styling():
styles = open('./styles/custom_barba.css', 'r').read()
return HTML(styles)
css_styling()