from IPython.display import HTML
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
Nichtlineare Schwingungsdynamik¶
Aufgabe 12¶
Ein Turm wird durch eine Windböe zu Schwingungen angeregt. Die Windlast sei durch eine dreiecksförmige Streckenlast $p(x)$ abgebildet, welche mit der Zeitfunktion
$ \begin{align} f(t) = \frac{1}{2} \left( 1-\cos\left( \frac{2\,\pi}{T_b} t \right) \right) \end{align} $
moduliert wird.
Lösung¶
Die Differentialgleichung der fremderregten Kontinuumsschwingung lautet
$ \begin{align} (EI\, w^{\prime\prime})^{\prime\prime} + \varrho A \ddot{w} = p(x)\,f(t) \, . \end{align} $
Der Separationsansatz $w(x,t) = \sum_{k=1}^\infty \varphi_k(x) \, q_k(t)$ liefert die Eigenformen in der Form:
from sympy import *
init_printing()
## variables
t = symbols('t', real=True)
x = symbols('x', real=True)
k = symbols('k', integer=True, positive=True)
## constants
EI, rho, A, l, lam, p0, Tb = symbols('EI, rho, A, l, lambda, p_0, T_b', real=True)
## functions
phi = Function('varphi')(lam,x)
p = Function('p')(x)
phi = (cos(lam*x/l)-cosh(lam*x/l)-(cos(lam)+cosh(lam))/(sin(lam)+sinh(lam))*(sin(lam*x/l)-sinh(lam*x/l)))
Eq(Function('varphi')(lam,x), phi)
%%html
<style>
table {float:left}
</style>
<!--use html to flush the tables to the left-->
Aus den Randbedingungen
$ \begin{align} w(x=0,t) &= 0 \\ w^\prime(x=0,t) &= 0 \\ M(x=l,t) &= 0 \\ Q(x=l,t) &= 0 \end{align} $
ergeben sich die Eigenwerte $\lambda_k$ :
| k | $\lambda_k$ |
|---|---|
| 1 | 1.8751 |
| 2 | 4.6941 |
| 3 | 7.8548 |
| >3 | $(2k-1)\pi/2$ |
Wird der Ansatz in die DGL eingesetzt, diese mit der $i$-ten Eigenform $\varphi_i$ multipliziert und beide Seiten über die Balkenlänge integriert, folgt
$ \begin{align} \sum_{k=1}^\infty \left( \int_0^l (EI\, \varphi_k^{\prime\prime})^{\prime\prime} \varphi_i \,\text{d}x \, q_k {\, +} \int_0^l \varrho A \varphi_k \varphi_i \,\text{d}x \, \ddot{q}_k \right) = \int_0^l p(x) \, \varphi_i \,\text{d}x \, f(t) \quad , \end{align} $
wobei alle Terme für $i\neq k$ aufgrund der Orthogonalitätsbeziehungen entfallen und somit
$ \begin{align} s_k q_k + m_k \ddot{q}_k & = \widetilde{p}_k \, f(t) \end{align} $
mit
$ \begin{align} s_k = \int_0^l (EI\varphi^{\prime \prime}_k)^{\prime \prime} \varphi_k \,\text{d}x \end{align} $
sk = EI*integrate(diff(phi,x,x,x,x)*phi,(x,0,l))
sk = sk.trigsimp()
Eq(Function('s_k')(lam,l,EI),sk)
$ \begin{align} m_k = \int_0^l \varrho \, A \, \varphi_k^2 \,\text{d}x \end{align} $
mk = rho*A*integrate(phi*phi,(x,0,l))
mk = mk.trigsimp()
Eq(Function('m_k')(lam,l,A,rho),mk)
p = p0*x/l
Eq(Function('p')(x), p)
f = (1-cos(2*pi*t/Tb))/2
Eq(Function('f')(t), f)
$ \begin{align} \widetilde{p}_k = \int_0^l p(x) \, \varphi_k \,\text{d}x \end{align} $
pk = integrate(p*phi,(x,0,l))
pk = pk.trigsimp()
Eq(Function('\widetilde{p}_k')(l,p0,lam), pk)
Die entsprechenden Werte ergeben sich für die ersten Eigenformen zu
k_max = 10
# save everything in arrays, too
lam_a = []
sk_a = []
mk_a = []
pk_a = []
str_out = ' k | $\lambda$ | $\sqrt{\\rho A / EI} \omega_k = \lambda^2$ | $s_k\,l^3\,/\,(EI)$ | $m_k$ / ($\\rho A l$) | $p_k$ / $(l p_0)$ \r\n'
str_out = str_out + '----|-------------|--------------------------------------------|-----------------------|---------------------------|-------------\r\n'
for i in range(k_max):
k = i+1
if k == 1:
lam_ = 1.8751
elif k == 2:
lam_ = 4.6941
elif k == 3:
lam_ = 7.8548
else:
lam_ = N((2*k-1)*pi/2.0)
sk_ = N(sk.subs(lam,lam_)*l**3/EI)
mk_ = N(mk.subs(lam,lam_)/rho/A/l)
pk_ = N(pk.subs(lam,lam_)/l/p0)
lam_a.append(lam_)
sk_a.append(sk_)
mk_a.append(mk_)
pk_a.append(pk_)
# k lam om sk mk pk
str_out = str_out + " {:2d} | {:8.4f} | {:15.4f} | {:14.6e} | {:17.4f} | {:12.8f}\r\n".format( k,lam_, lam_**2,sk_,mk_,pk_ )
# display the output as markdown table directly
from IPython.display import Markdown, display
display(Markdown(str_out))
und für jede Eigenform kann der "1-Massen-Schwinger"
$ \begin{align} s_k q_k + m_k \ddot{q}_k & = \widetilde{p}_k \, f(t) \end{align} $
dgl1 = Eq(symbols('s_k')*Function('q_k')(t)+symbols('m_k')*diff(Function('q_k')(t),t,t), symbols('\widetilde{p}_k')*Function('f')(t))
dgl1
dgl1 = dgl1.subs(Function('f')(t),(1-cos(2*pi/Tb*t))/2)
dgl1
gelöst werden. Unter der Annahme, dass sich die Böe periodisch wiederholt, lässt sich der eingeschwungene Zustand (partikuläre Lösung) berechnen und darstellen. Ansatz nach Art der rechten Seite:
$q_k(t) = C_1 + C_2 \, cos(2\,\pi\, t/T_b)$
C1, C2 = symbols('C_1, C_2', real=True)
dgl1 = dgl1.subs(Function('q_k')(t),C1+C2*cos(2*pi*t/Tb))
dgl1 = dgl1.doit()
dgl1 = dgl1.doit().expand()
dgl1
Der Koeffizientenvergleich liefert:
kvgl_1 = Eq(dgl1.lhs.as_independent(t)[0],dgl1.rhs.as_independent(t)[0])
C1_sol = solveset(kvgl_1,C1)
C1_sol = C1_sol.args[0]
kvgl_1, Eq(C1,C1_sol)
kvgl_2 = Eq(dgl1.lhs.coeff(cos(2*pi*t/Tb)),dgl1.rhs.coeff(cos(2*pi*t/Tb)))
C2_sol = solveset(kvgl_2,C2)
C2_sol = C2_sol.args[0]
kvgl_2, Eq(C2,C2_sol)
q_sol = C1_sol+C2_sol*cos(2*pi*t/Tb)
Eq(Function('q_k')(t),q_sol)
Man erkennt, dass der Anteil der ersten Eigenform dominant ist.
%matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np
import ipywidgets as widgets
# init figure
fig = plt.figure(figsize=(9.5, 8))
# set appropriate numerical values for the problem
Tb_ = 10.0
p0_ = 0.4
l_ = 100
EI_ = 1.24646*10**11
rhoA_ = 1500
# evaluate arrays of the generalised parameters numerically
sk_a_ = np.dot(sk_a,EI_/l_**3)
mk_a_ = np.dot(mk_a,rhoA_*l_)
pk_a_ = np.dot(pk_a,l_*p0_)
# scaling factors for plotting the eigenmodes
skal = [1E4, 1E6, 1E7]
# RHS function p(x)
ax1 = fig.add_subplot(2, 3, 1)
plt.xlabel("p (x)")
plt.ylabel("x")
plt.xlim(-0.5, p0_+0.5)
plt.ylim(0.0, l_*1.3)
plt.yticks([1],("l"))
plt.xticks([0, p0_],("0", "p0"))
plt.plot((p0_, p0_), (0, l_),color="black",linewidth=4.0)
plt.plot((0, 0), (0, l_),color="blue")
x_ = np.linspace(0,l_,10)
lam_p = lambdify(x, p.subs(p0,p0_).subs(l,l_), modules=['numpy']) # sympy expression -> numpy function
p_ = lam_p(x_)
plt.plot(p_,x_, color="blue")
num_ar = 7 # number of arrows to draw
# draw the arrows:
for aa in range(num_ar):
ax1.arrow(0, (aa+1)/num_ar*l_, (aa+1)/num_ar*p0_, 0,
length_includes_head=True,
head_width=2,
head_length=0.05,
color="blue")
# RHS time function f(t)
ax2 = fig.add_subplot(2, 3, 2)
plt.xlabel("t")
plt.ylabel("f (t)")
plt.xlim(0.0, Tb_)
plt.ylim(0.0, 1.3)
plt.yticks([1])
plt.xticks([0, Tb_],("0", "Tb"))
t_ = np.linspace(0,Tb_,100)
lam_f = lambdify(t, f.subs(Tb,Tb_), modules=['numpy']) # sympy expression -> numpy function
f_ = lam_f(t_)
plt.plot(t_,f_)
line, = ax2.plot((0, 0), (0, l_+0.3))
# plot the eigenmodes
# return the eigenfunction w_k(x,t):
def wk(kk):
# replace symbols in w = φ(x)*q(t) by numerical values :
wk = (phi*q_sol)\
.subs(lam,lam_a[kk-1])\
.subs(l,l_)\
.subs(symbols('\widetilde{p}_k'),pk_a_[kk-1])\
.subs(symbols('s_k'),sk_a_[kk-1])\
.subs(symbols('m_k'),mk_a_[kk-1])\
.subs(Tb,Tb_)
return wk
# put recurring plot commands in a function:
def plotmode(kk,axk,skal_):
plt.xlabel("w_{} (x,t)".format(kk))
plt.ylabel("x")
plt.xlim(-1, 1)
plt.ylim(0.0, l_*1.3)
plt.yticks([1],("l"))
x_ = np.linspace(0,l_,100)
lam_wk = lambdify([x, t], wk(kk), modules=['numpy']) # sympy expression -> numpy function
wk_ = lam_wk(x_, 0.0)
linewk, = axk.plot(skal_*wk_,x_,label='scal = {:.0e}'.format(skal_), color="black")
plt.legend()
return linewk,lam_wk
# first eigenmode
kk = 1 # eigenmode index
ax4 = fig.add_subplot(2, 3, 4)
linew1,lam_w1 = plotmode(kk,ax4,skal[kk-1])
# second eigenmode
kk = 2 # eigenmode index
ax5 = fig.add_subplot(2, 3, 5)
linew2,lam_w2 = plotmode(kk,ax5,skal[kk-1])
# third eigenmode
kk = 3 # eigenmode index
ax6 = fig.add_subplot(2, 3, 6)
linew3,lam_w3 = plotmode(kk,ax6,skal[kk-1])
# superposition of first three eigenmodes (scaled uniformly)
ax3 = fig.add_subplot(2, 3, 3)
plt.xlabel("w (x,t)")
plt.ylabel("x")
plt.xlim(-1, 1)
plt.ylim(0.0, l_*1.3)
plt.yticks([1],("l"))
x_ = np.linspace(0,l_,100)
linew, = ax3.plot(skal[0]*(lam_w1(x_,0.0) + lam_w2(x_,0.0) + lam_w3(x_,0.0)),x_ ,
label='scal = {:.0e}'.format(skal[0]),
color="black")
plt.legend()
plt.show()
# widgets / Animation
def on_value_change(play):
print("t = {}/100 * T_b".format(play))
# update line plot x/y - data
line.set_xdata((Tb_*play/100, Tb_*play/100))
linew1.set_xdata(skal[0]*lam_w1(x_,Tb_*play/100))
linew2.set_xdata(skal[1]*lam_w2(x_,Tb_*play/100))
linew3.set_xdata(skal[2]*lam_w3(x_,Tb_*play/100))
linew.set_xdata(skal[0]*(lam_w1(x_,Tb_*play/100) + lam_w2(x_,Tb_*play/100) + lam_w3(x_,Tb_*play/100)))
fig.canvas.draw() # redraw the new lines
inw = widgets.interact(on_value_change,
play=widgets.Play(value=0,min=0,max=100,
step=1,description="Press play"));
