Voorbeeld : Kniklast verend ingeklemde, uitkragende Euler-Bernoulli ligger

from IPython.display import Image
from IPython.core.display import HTML
Image(filename = "pictures/maple-python-cantilever-spring.jpg", width=400)

(c) 2024 : Hans Welleman

Gegeven:

- doorsnede, buitgstijfheid $EI$

- lengte $L$

---

Buigingsknik voor een 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 en \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} Algemene oplossing: \begin{equation} w(x)=C_1+C_2x+C_3sin( \alpha x) +C_4 cos( \alpha x) \end{equation} met: \begin{equation} \alpha^2=\frac{F}{EI} \quad en \quad Sz(x)=-FC_2 \end{equation} \begin{equation} en \quad \quad \rho=\frac{kL}{EI} \end{equation}

---

Vraag:
Vind de laagste vier kniklasten en de bijbehorende knikvormen.

import numpy as np
import sympy as sym
import scipy.optimize as so
from sympy.plotting import plot
# --------------------------------------------------------------------
# NOOT: controleer:
#       - sympy is up to date
#       - matplotlib is up to date
#       - scipy is up to date (pip install scipy)
#       - installeer ipympl (b.v. met pip install ipywidgets==8.0.4)
#
# eenheden kN, m
#
# 2024 (c) Hans Welleman
# --------------------------------------------------------------------
w   = sym.symbols('w', cls=sym.Function)
C1, C2, C3, C4 = sym.symbols('C1 C2 C3 C4')
x    = sym.symbols('x', real=True, positive=True)
EI, alpha, A, L, k, rho = sym.symbols('EI alpha A L k rho',real=True, positive=True)
# vaste waarden voor parameters:
k = rho*EI/L

DV-oplossing

# differentiaal vergelijking:
DV = sym.Eq(sym.diff(w(x),x,4)+(alpha**2)*sym.diff(w(x),x,2),0)
w  = sym.dsolve(DV,w(x))
w = w.rhs
# algemene oplossing:
display(w)

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

# afgeleide vergelijkinen voor de randvoorwaarden:
phi   = -sym.diff(w , x)
kappa =  sym.diff(phi , x)
M     = EI * kappa
Sz    = sym.simplify(sym.diff(M , x)-(alpha**2*EI)*sym.diff(w,x))
display(Sz)

$\displaystyle - C_{2} EI \alpha^{2}$

# Randvoorwaarden:
eq1  = sym.Eq(  w.subs(x , 0) , 0)
eq2  = sym.Eq(  M.subs(x,0) - k*phi.subs(x , 0) , 0)
eq3 =  sym.Eq(sym.simplify(Sz.subs(x , L)) , 0)
eq4 =  sym.Eq(sym.simplify(M.subs(x , L)) , 0)
display(eq1, eq2, eq3, eq4)

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

$\displaystyle C_{4} EI \alpha^{2} - \frac{EI \rho \left(- C_{2} - C_{3} \alpha\right)}{L} = 0$

$\displaystyle - C_{2} EI \alpha^{2} = 0$

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

# schrijf stelsel in matrixvorm voor vinden van de determinant
A, b = sym.linear_eq_to_matrix((eq1 , eq2 , eq3 , eq4) ,
               (C1  , C2  , C3  , C4))
display(A)
display(b)
# ... SOME WORK TO SORT OUT ...
# nog uitzoeken hoe de vergelijkingen kunnen worden opgeschoond ....
char_eq = sym.simplify( A.det()*L/(alpha**5*EI**3))
display(char_eq)

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

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

$\displaystyle - L \alpha \sin{\left(L \alpha \right)} + \rho \cos{\left(L \alpha \right)}$

Lval    = 1
rhoval  = 2  # try 1 .. 4
# bepaal de eerste vier wortels
# plot karakteristieke vergelijking:
%matplotlib widget
# plot oplossing voor (alpha L) :
AA = char_eq.subs(L,Lval).subs(rho,rhoval),(alpha,0,10)
# display(AA)
plot(AA,title='Karakteristieke vergelijking',ylim=[-1,1],xlabel=r'$\alpha$L',ylabel=r' kar_verg($\alpha$L)',adaptive=False);

Figure

Kniklasten

# bepaal de nulpunten van de karakteristieke vergelijking (=kniklasten)
print("determinant:");
display(char_eq)
f = sym.lambdify(alpha, char_eq.subs(L,Lval).subs(rho,rhoval))
# numerieke oplossing voor 4 startwaarden:
# display(f)
sol = so.fsolve(f,[1.7,4.5,7.3,10])
x1 = sol[0]
x2 = sol[1]
x3 = sol[2]
x4 = sol[3]
print('\N{greek small letter alpha}L1, \N{greek small letter alpha}L2, \N{greek small letter alpha}L3, \N{greek small letter alpha}L4')
print(x1, x2, x3, x4)
# kniklasten
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 - L \alpha \sin{\left(L \alpha \right)} + \rho \cos{\left(L \alpha \right)}$

αL1, αL2, αL3, αL4
1.0768739969518588 3.6435971661455393 6.578333731545392 9.629560343296307
Fk1, Fk2, Fk3, Fk4
1.15965760531107*EI/L**2 13.2758003091438*EI/L**2 43.2744746835879*EI/L**2 92.7284324051849*EI/L**2

Knikvormen

# vind de verhoudingen tussen C1, C2 en C4 t.o.v. C3 (sinus-term)
sol = sym.solve((eq1,eq2,eq3) ,(C1  , C2  , C4))
# display(sol)
# substitueer het resultaat in w en schaal C3 met 1/2
MyMode=w.subs(sol).subs(C3,1/2)
# substitueer de gevonden alpha*L per knikvorm in w1 .. w4:
w1=MyMode.subs(alpha,x1/L).subs(L,1).subs(rho,rhoval),(x,0,1)
w2=MyMode.subs(alpha,x2/L).subs(L,1).subs(rho,rhoval),(x,0,1)
w3=MyMode.subs(alpha,x3/L).subs(L,1).subs(rho,rhoval),(x,0,1)
w4=MyMode.subs(alpha,x4/L).subs(L,1).subs(rho,rhoval),(x,0,1)

# plot de knikvormen in 1 grafiek:
atitle = "Knikvormen 1 t/m 4 voor rho=" + str(rhoval)
p=plot(w1,w2,w3,w4,title=atitle,xlabel=r'x/L',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

Vergelijk oplossing met standaardformule

\begin{equation} \frac{1}{F_k}=\frac{1}{F_k star} + \frac{1}{F_k Euler} \end{equation}\begin{equation} met: \quad \quad F_k star=\frac{k}{L} \quad en \quad F_k Euler=\frac{\pi^2EI}{4L^2} \end{equation}\begin{equation} en \quad \quad \rho=\frac{kL}{EI} \end{equation}

# kniklast op basis van ingenieursaanpak:
Feuler = sym.pi**2*EI/(4*L**2)
Frigid = k/L
Fk     = sym.simplify(1/( (1/Frigid)+(1/Feuler)))
print('Fk - star:')
display(Frigid)
print('Fk - Euler:')
display(Feuler)
print('Fk:')
display(Fk)

Fk - star:

$\displaystyle \frac{EI \rho}{L^{2}}$

Fk - Euler:

$\displaystyle \frac{\pi^{2} EI}{4 L^{2}}$

Fk:

$\displaystyle \frac{\pi^{2} EI \rho}{L^{2} \cdot \left(4 \rho + \pi^{2}\right)}$

print("Benaderingsoplossing voor kniklast Fk voor rho=" + str(rhoval))
display(Fk.subs(rho,rhoval))
(Fk.subs(rho,rhoval)).evalf()

Benaderingsoplossing voor kniklast Fk voor rho=2

$\displaystyle \frac{2 \pi^{2} EI}{L^{2} \cdot \left(8 + \pi^{2}\right)}$

$\displaystyle \frac{1.10462483439059 EI}{L^{2}}$

FF = (1/(4+(sym.pi**2/rhoval)))*(sym.pi**2*EI/L**2)
display(FF)
print(FF.evalf())

$\displaystyle \frac{\pi^{2} EI}{L^{2} \cdot \left(4 + \frac{\pi^{2}}{2}\right)}$

1.10462483439059*EI/L**2