EXAMPLE : Buckling load and modes for Euler-Bernouilli beam - Type 2

from IPython.display import Image
from IPython.core.display import HTML
Image(filename = "pictures/Euler-buckling.jpg", width=700)

(c) 2024 : Hans Welleman

Given:

- cross section with bending stiffness $EI$

- constant compressive force $F$

- length $l$

---

Flexural Buckling for an Euler-Bernoulli model: \begin{equation} EI\frac{\mathrm{d}^4 w(x)}{\mathrm{d} x^4}+F \frac{\mathrm{d}^2 w(x)}{\mathrm{d} x^2} = 0 \quad and \quad Sz(x)=-EI\frac{\mathrm{d}^3 w(x)}{\mathrm{d} x^3}-F \frac{\mathrm{d} w(x)}{\mathrm{d} x} \end{equation} General Solution: \begin{equation} w(x)=C_1+C_2x+C_3sin( \alpha x) +C_4 cos( \alpha x) \end{equation} with: \begin{equation} \beta^2=\frac{F}{EI} \quad and \quad Sz(x)=-FC_2=-\beta^2 EI C_2 \end{equation}

---

Question:
Find the lowest four buckling loads and corresponding buckling modes.

import numpy as np
import sympy as sym
import scipy.optimize as so
from sympy.plotting import plot
# --------------------------------------------------------------------
# NOTE: check
#       - sympy is up to date
#       - matplotlib is up to date
#       - scipy is up to date (pip install scipy)
#       - install ipympl (e.g. pip install ipywidgets==8.0.4)
#
# units kN, m
#
# 2024 (c) Hans Welleman
# --------------------------------------------------------------------
W   = sym.symbols('W', cls=sym.Function)
C1, C2, C3, C4 = sym.symbols('C1 C2 C3 C4')
x, F= sym.symbols('x F', real=True, positive=True)
EI, beta, L = sym.symbols('EI beta L',real=True, positive=True)
# fixed values for parameters:
# beta^2 = F/EI

ODE-solution

# ODE per field:
DV = sym.Eq(W(x).diff(x,4)+(beta**2)*W(x).diff(x,2),0)
W  = sym.dsolve(DV,W(x))
display(W)
W = W.rhs
# general solution:
display(W)

$\displaystyle W{\left(x \right)} = C_{1} + C_{2} x + C_{3} \sin{\left(\beta x \right)} + C_{4} \cos{\left(\beta x \right)}$

$\displaystyle C_{1} + C_{2} x + C_{3} \sin{\left(\beta x \right)} + C_{4} \cos{\left(\beta x \right)}$

# needed derivatives for BC or IC:
phi   = -sym.diff(W , x)
kappa =  sym.diff(phi , x)
M     = EI * kappa
Sz    =  sym.diff(M , x)-beta**2*EI*sym.diff(W, x)

# Boundary and/or Interface Conditions - system of equations -
# type 4 buckling (see picture at the top)
eq1  = sym.Eq(  W.subs(x , 0) , 0)
eq2  = sym.Eq(phi.subs(x , 0) , 0)
eq3  = sym.Eq( Sz.subs(x , L) , 0)
eq4  = sym.Eq(phi.subs(x , L) , 0)
display(eq1, eq2, eq3, eq4)

$\displaystyle C_{1} + C_{4} = 0$

$\displaystyle - C_{2} - C_{3} \beta = 0$

$\displaystyle - EI \beta^{2} \left(C_{2} + C_{3} \beta \cos{\left(L \beta \right)} - C_{4} \beta \sin{\left(L \beta \right)}\right) + EI \left(C_{3} \beta^{3} \cos{\left(L \beta \right)} - C_{4} \beta^{3} \sin{\left(L \beta \right)}\right) = 0$

$\displaystyle - C_{2} - C_{3} \beta \cos{\left(L \beta \right)} + C_{4} \beta \sin{\left(L \beta \right)} = 0$

# create matrix notation and find determinant:
A, b = sym.linear_eq_to_matrix((eq1 , eq2 , eq3 , eq4) ,
               (C1  , C2  , C3  , C4))
display(A)
display(b)
kar_vgl = sym.simplify( A.det()) / (EI)
display(kar_vgl)

$\displaystyle \left[\begin{matrix}1 & 0 & 0 & 1\\0 & -1 & - \beta & 0\\0 & - EI \beta^{2} & 0 & 0\\0 & -1 & - \beta \cos{\left(L \beta \right)} & \beta \sin{\left(L \beta \right)}\end{matrix}\right]$

$\displaystyle \left[\begin{matrix}0\\0\\0\\0\end{matrix}\right]$

$\displaystyle - \beta^{4} \sin{\left(L \beta \right)}$

Lval = 1
# find first four roots
# plot characteristic equations:
%matplotlib widget
AA = kar_vgl.subs(L,Lval),(beta,0,20)
plot(AA,title='Characteristic equation',ylim=[-1,1],xlabel=r'$\beta$L',ylabel=r' char_eq($\beta$L)',adaptive=False);

Figure

Buckling Loads

# Find roots of the characteristic equation (=buckling load)
print("determinant:");
display(kar_vgl)
f = sym.lambdify(beta, kar_vgl.subs(L,Lval))
# numerical solution with 4 start values:
sol = so.fsolve(f,[3,6,9,13])
x1 = sol[0]
x2 = sol[1]
x3 = sol[2]
x4 = sol[3]
print('\N{greek small letter beta}L1, \N{greek small letter beta}L2, \N{greek small letter beta}L3, \N{greek small letter beta}L4')
print(x1, x2, x3, x4)
# buckling loads
Fk1 = x1**2*EI/L**2
Fk2 = x2**2*EI/L**2
Fk3 = x3**2*EI/L**2
Fk4 = x4**2*EI/L**2
print('Fk1, Fk2, Fk3, Fk4')
print(Fk1, Fk2, Fk3, Fk4)

determinant:

$\displaystyle - \beta^{4} \sin{\left(L \beta \right)}$

βL1, βL2, βL3, βL4
3.141592474350653 6.283185307396782 9.424777960786972 12.566370614360128
Fk1, Fk2, Fk3, Fk4
9.86960327489666*EI/L**2 39.4784176070868*EI/L**2 88.8264396101358*EI/L**2 157.913670417454*EI/L**2

Buckling Modes

# find ratio  of C1, C2 and C3 with respect to C4 which belongs to the cos(alpha x) term
sol = sym.solve((eq1,eq3,eq4) ,(C1  , C2  , C3))
# display(sol)
# substitute this in W and scale C4 to 2
MyMode=W.subs(sol).subs(C4,2)
# substitute each beta*L per buckling load:
w1=MyMode.subs(beta,x1/L).subs(L,1),(x,0,1)
w2=MyMode.subs(beta,x2/L).subs(L,1),(x,0,1)
w3=MyMode.subs(beta,x3/L).subs(L,1),(x,0,1)
w4=MyMode.subs(beta,x4/L).subs(L,1),(x,0,1)

# plot the buckling modes:
p=plot(w1,w2,w3,w4,title='Type 2, Buckling Modes 1 - 4',ylabel='w(x)',legend=True,adaptive=False, show=False);
p[0].label='mode 1' 
p[1].label='mode 2'
p[2].label='mode 3'
p[3].label='mode 4'
p.show()

Figure