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

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
#       - installeer ipympl (b.v. met pip install ipywidgets==8.0.4)
#
# eenheden kN, m
#
# 2024 (c) Hans Welleman
# --------------------------------------------------------------------
w1, w2 = sym.symbols('w1 w2', cls=sym.Function)
C1, C2, C3, C4, C5, C6, C7, C8 = sym.symbols('C1 C2 C3 C4 C5 C6 C7 C8')
x    = sym.symbols('x', positive=True)
F    = sym.symbols('F', real=True, positive=True)
a, b, q, EI1, EI2 = sym.symbols('a b q EI1 EI2', real=True, positive=True)
# vaste waarden voor parameters:
a   = sym.nsimplify(10*(1/3))
b   = sym.nsimplify(10*(2/3))
EI1 = 5000
EI2 = EI1/2

A   = sym.symbols('A')
w   = A*sym.sin(sym.pi*x/(a+b))
# display(w)
EE1 = (1/2)*EI1*sym.diff(w,x,2)**2
EE2 = (1/2)*EI2*sym.diff(w,x,2)**2
Ev  = sym.integrate(EE1,(x,0,a)) + sym.integrate(EE2,(x,a,a+b))
Atot= F*sym.integrate((1/2)*sym.diff(w,x)**2,(x,0,a+b)) + \
      F*sym.integrate((1/2)*sym.diff(w,x)**2,(x,0,a))
eq  = sym.Eq(Ev,Atot)
sol = sym.solve(eq,F)
FR  = sol[0]
print("Kniklast volgens Rayleigh:")
display(FR)

Kniklast volgens Rayleigh:

# differentiaalvergelijkingen per veld:
DV1   = sym.Eq(EI1*w1(x).diff(x,4)+2*F*w1(x).diff(x,2),0)
DV2   = sym.Eq(EI2*w2(x).diff(x,4)+F*w2(x).diff(x,2),0)
w1,w2 = sym.dsolve([DV1,DV2],[w1(x),w2(x)])
w1 = w1.rhs
w2 = w2.rhs
# algemene oplossing per veld:
# display(w1,w2)

# afgeleide functies:
phi1   = -sym.diff(w1 , x)
kappa1 =  sym.diff(phi1 , x)
M1     = EI1 * kappa1
Sz1    =  sym.diff(M1 , x)-2*F*sym.diff(w1,x)

phi2   = -sym.diff(w2 , x)
kappa2 =  sym.diff(phi2 , x)
M2     = EI2 * kappa2
Sz2    =  sym.diff(M2 , x)-F*sym.diff(w2,x)

# rand- en overgangsvooraarden - stelsel vergelijkingen -
eq1  = sym.Eq(  w1.subs(x , 0) , 0)
eq2  = sym.Eq(  M1.subs(x , 0) , 0)

eq3  = sym.Eq(  w1.subs(x , a) - w2.subs(x , a)  ,0)
eq4  = sym.Eq(  M1.subs(x , a) - M2.subs(x , a)  ,0)
eq5  = sym.Eq(phi1.subs(x , a) - phi2.subs(x , a),0)
eq6  = sym.Eq( Sz1.subs(x , a) - Sz2.subs(x , a) ,0)

eq7 = sym.Eq(w2.subs(x , a+b) , 0)
eq8 = sym.Eq(M2.subs(x , a+b) , 0)
# display(eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8)

# schrijf vergelijkingen als stelsel en bepaal de determinant:
A, b = sym.linear_eq_to_matrix((eq1 , eq2 , eq3 , eq4 , eq5 , eq6 , eq7 , eq8) ,
               (C1  , C2  , C3  , C4  , C5  , C6  , C7  , C8))
display(A)
display(b)
kar_vgl = sym.simplify( A.det()/EI1**3 )
display(kar_vgl)

# zoek eerste vier wortels
# plot karakteristieke vergelijking:
%matplotlib widget
# plot oplossing kar_vgl:
AA = kar_vgl,(F,1,7000)
plot(AA,title='karakteristieke vergelijking',ylim=[-1,1],xlabel='F',ylabel=' kar_vgl(F)',adaptive=False);

# Bepaal de wortels van de karakteristieke vergelijking (=kniklast)
print("determinant:");
display(kar_vgl)
g = sym.lambdify(F, kar_vgl)
# numeriek oplossen met 4 startwaarden:
sol = so.fsolve(g,[200,800,3500,6500])
Fk1 = sol[0]
Fk2 = sol[1]
Fk3 = sol[2]
Fk4 = sol[3]
print('Fk1, Fk2, Fk3, Fk4')
print(Fk1, Fk2, Fk3, Fk4)

determinant:

Fk1, Fk2, Fk3, Fk4
189.87349873751302 1060.4892969033594 3761.492820704937 6371.916754767171

Voorbeeld : Knik, kolom met 2 puntlasten F¶

Rayleigh benadering¶

DV-oplossing¶

Kniklasten¶