In [1]:
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>''')
Out[1]:

Lehrstuhl Technische Dynamik, Juniorprofessur Fluid-Struktur Kopplung

                        

Aufgabe 10

Berechnen Sie die Schwingungen im Antriebsstrang einer elektrischen Lokomotive, die sich mit konstanter Geschwindigkeit $v$ bewegt. Betrachten Sie die Schubstangen als Balken, welche unter Wirkung einer definierten Längskraft

$F_s = \frac{M}{r} \cos\left( \frac{v}{r_a} t\right)$

Querschwingungen ausführen können. Welche relevanten Schwungungsphänomene treten auf und wie können diese behandelt werden.

Erzwungene Longitudinalschwingungen

${F}_{s}=\frac{M}{r}\mathrm{cos}\left(\frac{v}{{r}_{a}}t\right)$, $\qquad$ $\Omega = \frac{v}{r_a}$

Separationsansatz: \begin{align} u\left(x,t\right)&= U(x) T(x) \\ u\left(x,t\right)&=\left[\stackrel{-}{A}\mathrm{cos}\left(\frac{\omega }{c}x\right)+\stackrel{-}{B}\mathrm{sin}\left(\frac{\omega }{c}x\right)\right] \cos(\omega t - \alpha) \end{align}

Wellengleichung: \begin{align} \frac{{\partial }^{2}u}{\partial {x}^{2}}=\frac{1}{{c}^{2}}\frac{{\partial }^{2}u}{\partial {t}^{2}} \end{align}

Randbedingungen: \begin{align} N\left(l,t\right)&=EAu^{\prime}\left(l,t\right)={F}_{s}\left(t\right) \\ N\left(0,t\right)&=EAu^{\prime}\left(0,t\right)={F}_{s}\left(t\right) \end{align}

Homogene Lösung: Beidseitig gelenkig gelagerter Stab: Forderung: $\sin\left( \frac{\omega_k}{c} l \right) = 0$ $\qquad$ mit $c = \sqrt{E/\varrho}$ $\qquad \leadsto$ $\omega_k = k \pi \sqrt{{E}/({\varrho l^2})}$

\begin{align*} u_h(x,t) &= {U}_{k}\left(x\right)\cdot \left({C}_{k}\mathrm{cos}\left({\omega }_{k}t\right)+{D}_{k}\mathrm{sin}\left({w}_{k}t\right)\right) \\ \text{bzw.} \\ u_h(x,t) &= {U}_{k}\left(x\right)\cdot {E}_{k}\cdot \mathrm{cos}\left({\omega }_{k}t-{\alpha }_{k}\right) \\ \text{mit} \\ U_k(x) &= \cos\left( \frac{k\pi x}{l} \right) \end{align*}

Resonanzanregung bei $\Omega = \omega_k = k \pi \sqrt{{E}/({\varrho l^2})}$

Partikuläre Lösung:

\begin{align*} {u}_{p}\left(x,t\right)={U}_{p}\left(x\right)\cdot \mathrm{cos}\left(\Omega t\right) \end{align*}

in Wellengleichung einsetzen:

\begin{align*} &{U}^{\prime\prime}_{p}\mathrm{cos}\left(\Omega t\right)=-\frac{1}{{c}^{2}}{U}_{p}\left(x\right){\Omega }^{2}\mathrm{cos}\left(\Omega t\right) \\ \leadsto \qquad &{U}^{\prime\prime}_{p}\left(x\right)+\frac{{\Omega }^{2}}{{c}^{2}}{U}_{p}\left(x\right)=0 \end{align*}

Lösung:

\begin{align*} {U}_{p}\left(x\right)&={B}_{{1}}\mathrm{cos}\left(\frac{\Omega }{c}x\right)+{B}_{2}\mathrm{sin}\left(\frac{\Omega }{c}x\right) \\ \leadsto \qquad u_p(x,t) &= \left[ {B}_{{1}}\mathrm{cos}\left(\frac{\Omega }{c}x\right)+{B}_{2}\mathrm{sin}\left(\frac{\Omega }{c}x\right) \right] \cos(\Omega t) \end{align*}

Auswertung der Randbedingungen

\begin{align*} u^\prime_p (x,t) &= \left[ -B_1 \frac{\Omega}{c} \sin\left(\frac{\Omega}{c} x\right) + B_2 \frac{\Omega}{c} \cos\left( \frac{\Omega}{c} x \right) \right] \cos(\Omega t) \end{align*}

1) $\qquad u^\prime_p (l,t) = \left[ -B_1 \frac{\Omega}{c} \sin\left(\frac{\Omega}{c} l\right) + B_2 \frac{\Omega}{c} \cos\left( \frac{\Omega}{c} l\right) \right] \cos(\Omega t) = \frac{F_s(t)}{EA} = \frac{M}{EAr} cos(\Omega t)$

2) $\qquad u^\prime_p (0,t) = B_2 \frac{\Omega}{c} \cos(\Omega t) = \frac{F_s(t)}{EA} = \frac{M}{EAr} cos(\Omega t)$

$\qquad\leadsto B_2 = \frac{M}{EAr} \frac{c}{\Omega} = \frac{M}{EAr} \frac{\sqrt{E/\varrho}\,r_a}{v}$

in 1)

\begin{align*} \,- B_1 \frac{\Omega}{c} \sin\left( \frac{\Omega}{c} l \right) + \frac{M}{EAr} \cos\left( \frac{\Omega}{c} l \right) = \frac{M}{EAr} \\ \,- B_1 \frac{\Omega}{c} \frac{EAr}{M} \sin\left( \frac{\Omega}{c} l \right) = 1- \cos\left( \frac{\Omega}{c} l \right) \\ B_1 = \frac{\cos\left( \frac{\Omega}{c} l \right) -1}{\sin\left( \frac{\Omega}{c} l \right)} \frac{M}{EAr} \frac{\sqrt{E/\varrho}\,r_a}{v} \end{align*}

Eingeschwungene Lösung:

\begin{equation} u_p(x,t) = \frac{M r_a \sqrt{E/\varrho}}{vEAr} \left[ \frac{\cos(\Omega l/c)-1}{\sin(\Omega l /c)} \cos\left( \frac{\Omega}{c} x \right) + \sin\left( \frac{\Omega}{c} x \right)\right] \cos(\Omega t) \end{equation}\begin{equation} u_p(x,t) = \frac{M r_a \sqrt{E/\varrho}}{vEAr} \left[ -\tan\left( \frac{\Omega l}{2c} \right)\cos\left( \frac{\Omega}{c} x \right) + \sin\left( \frac{\Omega}{c} x \right)\right] \cos(\Omega t) \end{equation}

Erzwungene Transversalschwingungen

Differentialgleichung (Balken schubstarr, Rotationsträgheit vernachlässigt, $EI$ = konst.) :

\begin{equation} \varrho A \ddot{w} = \frac{\partial^2M_b}{\partial x^2} \end{equation}

homogenes Problem: Bernoulli: $M_b(x) = -EI w^{\prime\prime}$

$\leadsto \varrho A \ddot{w} = -EI w^{\prime\prime\prime\prime}$

hier inhomogenes Problem:

$M_b(x) = -EI w^{\prime\prime} + F_s(t)\,w$

\begin{align*} \leadsto \varrho A \ddot{w} = -EI w^{\prime\prime\prime\prime} + F_s(t)\,w^{\prime\prime} \\ \leadsto EI w^{\prime\prime\prime\prime} - \frac{M}{r} \cos\left(\frac{v}{r_a} t\right) w^{\prime\prime} + \varrho A \ddot{w} = 0 \end{align*}
In [2]:
from sympy import *
init_printing()
## variables
t = symbols('t', real=True)
x = symbols('x', real=True)
## functions
w = Function('w')(x,t)
W = Function('W')(x)
T = Function('T')(t)
## constants
k, pi, l, v, ra, EI, rho, A, M, r = symbols('k, pi, l, v, r_a, EI, varrho, A, M, r', real=True)

dgl = Eq(EI*diff(w,x,x,x,x)-M/r*cos(v/ra*t)*diff(w,x,x)+rho*A*diff(w,t,t),0)
dgl
Out[2]:
$$A \varrho \frac{\partial^{2}}{\partial t^{2}} w{\left (x,t \right )} + EI \frac{\partial^{4}}{\partial x^{4}} w{\left (x,t \right )} - \frac{M \cos{\left (\frac{t v}{r_{a}} \right )} \frac{\partial^{2}}{\partial x^{2}} w{\left (x,t \right )}}{r} = 0$$

Randbedingungen:

\begin{align} w(x=0,t) &= 0 \\ w(x=l,t) &= 0 \\ M_b(x=0,t) &= 0 \leadsto w^{\prime \prime}(x=0,t)=0 \\ M_b(x=l,t) &= 0 \leadsto w^{\prime \prime}(x=l,t)=0 \\ \end{align}

Separationsansatz:

In [3]:
Eq(w,W*T)
Out[3]:
$$w{\left (x,t \right )} = T{\left (t \right )} W{\left (x \right )}$$

Eigenformen für zweifache Lagerung mit $k = 1,2,3,\ldots$

In [4]:
W
Eq(W,sin(k*pi*x/l))
Out[4]:
$$W{\left (x \right )} = \sin{\left (\frac{k \pi x}{l} \right )}$$
In [5]:
sol = Eq(w,W*T)
sol = Eq(sol.lhs,sol.rhs.replace(W,sin(k*pi*x/l)))
sol
Out[5]:
$$w{\left (x,t \right )} = T{\left (t \right )} \sin{\left (\frac{k \pi x}{l} \right )}$$
In [6]:
w_sol = sol.rhs
Eq(diff(w,t,t),diff(w_sol,t,t))
Out[6]:
$$\frac{\partial^{2}}{\partial t^{2}} w{\left (x,t \right )} = \sin{\left (\frac{k \pi x}{l} \right )} \frac{d^{2}}{d t^{2}} T{\left (t \right )}$$
In [7]:
Eq(diff(w,x,x),diff(w_sol,x,x))
Out[7]:
$$\frac{\partial^{2}}{\partial x^{2}} w{\left (x,t \right )} = - \frac{k^{2} \pi^{2} T{\left (t \right )} \sin{\left (\frac{k \pi x}{l} \right )}}{l^{2}}$$
In [8]:
Eq(diff(w,x,x,x,x),diff(w_sol,x,x,x,x))
Out[8]:
$$\frac{\partial^{4}}{\partial x^{4}} w{\left (x,t \right )} = \frac{k^{4} \pi^{4} T{\left (t \right )} \sin{\left (\frac{k \pi x}{l} \right )}}{l^{4}}$$

Einsetzen:

In [9]:
dgl = Eq(dgl.lhs.replace(w,w_sol),0)
dgl
Out[9]:
$$A \varrho \frac{\partial^{2}}{\partial t^{2}} T{\left (t \right )} \sin{\left (\frac{k \pi x}{l} \right )} + EI \frac{\partial^{4}}{\partial x^{4}} T{\left (t \right )} \sin{\left (\frac{k \pi x}{l} \right )} - \frac{M \cos{\left (\frac{t v}{r_{a}} \right )} \frac{\partial^{2}}{\partial x^{2}} T{\left (t \right )} \sin{\left (\frac{k \pi x}{l} \right )}}{r} = 0$$

$sin(k\pi x/l)$ ausklammern:

In [10]:
dgl = Eq(dgl.lhs.doit().expand().collect(sin(k*pi*x/l)),0)
dgl
Out[10]:
$$\left(A \varrho \frac{d^{2}}{d t^{2}} T{\left (t \right )} + \frac{EI k^{4} \pi^{4} T{\left (t \right )}}{l^{4}} + \frac{M k^{2} \pi^{2} T{\left (t \right )} \cos{\left (\frac{t v}{r_{a}} \right )}}{l^{2} r}\right) \sin{\left (\frac{k \pi x}{l} \right )} = 0$$

Es ist eine Lösung für alle $x$ gesucht. Daher muss der linke Klammerausdruck verschwinden:

In [11]:
dgl2 = Eq(dgl.lhs.coeff(sin(k*pi*x/l)).collect(T),0)
dgl2
Out[11]:
$$A \varrho \frac{d^{2}}{d t^{2}} T{\left (t \right )} + \left(\frac{EI k^{4} \pi^{4}}{l^{4}} + \frac{M k^{2} \pi^{2} \cos{\left (\frac{t v}{r_{a}} \right )}}{l^{2} r}\right) T{\left (t \right )} = 0$$

da $A\varrho \neq 0$, folgt

In [12]:
dgl3 = Eq((dgl2.lhs/A/rho).expand().collect(T),dgl2.rhs/A/rho)
dgl3
Out[12]:
$$\left(\frac{EI k^{4} \pi^{4}}{A l^{4} \varrho} + \frac{M k^{2} \pi^{2} \cos{\left (\frac{t v}{r_{a}} \right )}}{A l^{2} r \varrho}\right) T{\left (t \right )} + \frac{d^{2}}{d t^{2}} T{\left (t \right )} = 0$$

Einführen der ,,dimensionslosen Zeit''

\begin{align} \tau &= \frac{v}{2r_a}t \\ \text{mit}\qquad \partial T/\partial t &= \frac{v}{2r_a} \partial T/\partial \tau \end{align}

liefert

In [13]:
tau = symbols('tau', real=True)
T2 = Function('T')(tau)
dgl4 = Eq(dgl3.lhs.replace(diff(T,t,t),v**2/(2*ra)**2*diff(T2,tau,tau)).replace(T,T2).replace(v/ra*t,2*tau),0)
dgl4
Out[13]:
$$\left(\frac{EI k^{4} \pi^{4}}{A l^{4} \varrho} + \frac{M k^{2} \pi^{2} \cos{\left (2 \tau \right )}}{A l^{2} r \varrho}\right) T{\left (\tau \right )} + \frac{v^{2} \frac{d^{2}}{d \tau^{2}} T{\left (\tau \right )}}{4 r_{a}^{2}} = 0$$

bzw. mit $v^2/(4r_a^2) \neq 0$

In [14]:
dgl5 = Eq((dgl4.lhs*4*ra**2/v**2).expand().collect(T2),dgl4.rhs*4*ra**2/v**2)
dgl5
Out[14]:
$$\left(\frac{4 EI k^{4} \pi^{4} r_{a}^{2}}{A l^{4} v^{2} \varrho} + \frac{4 M k^{2} \pi^{2} r_{a}^{2} \cos{\left (2 \tau \right )}}{A l^{2} r v^{2} \varrho}\right) T{\left (\tau \right )} + \frac{d^{2}}{d \tau^{2}} T{\left (\tau \right )} = 0$$

Der Abgleich mit der Mathieu-DGL

\begin{align} \frac{\partial^2 \varphi}{\partial \tau^2} + (\lambda - 2\gamma \cos(2\tau))\varphi = 0 \end{align}

liefert

In [15]:
lam = symbols('lambda', real=True)
gam = symbols('gamma', real=True)
lameq = Eq(lam,dgl5.lhs.coeff(T2).as_independent(tau)[0])
gameq = Eq(-gam,dgl5.lhs.coeff(T2).as_independent(tau)[1].coeff(cos(2*tau))/2)
display(lameq)
display(gameq)
$$\lambda = \frac{4 EI k^{4} \pi^{4} r_{a}^{2}}{A l^{4} v^{2} \varrho}$$
$$- \gamma = \frac{2 M k^{2} \pi^{2} r_{a}^{2}}{A l^{2} r v^{2} \varrho}$$

Die möglichen Punkte auf der Strutt'schen Stabilitätskarte liegen auf einer Ursprungsgeraden mit der Steigung

$p = -\gamma / \lambda$

In [16]:
p = symbols('p', real=True)
lam_sol = solveset(lameq,lam).args[0]
gam_sol = solveset(gameq,gam).args[0]
Eq(p,-gam_sol/lam_sol)
Out[16]:
$$p = \frac{M l^{2}}{2 EI k^{2} \pi^{2} r}$$

Gut für große Stabilitätsbereiche:

1) großes $r$, $EI$

2) kleines $M$, $l$

In [17]:
import matplotlib.pyplot as plt
import numpy as np
import ipywidgets as widgets

# plot data
d = np.array([
                [0, 0, 1, 1, 4, 4, 9, 9, 16, 16, 25, 25, 36],
                [1, -0.4551386, -0.1102488, 1.8591081, 3.9170248, 4.371301, 9.0477393, 9.0783688, 16.0329701, 16.0338323, 25.0208408, 25.0208543, 36.0142899],
                [2, -1.5139569, -1.3906765, 2.3791999, 3.6722327, 5.1726651, 9.14062277, 9.3703225, 16.127688, 16.1412038, 25.083349, 25.0837778, 36.057207],
                [3, -2.8343919, -2.7853797, 2.5190391, 3.276922, 6.0451969, 9.2231328, 9.91155063, 16.2727012, 16.3387207, 25.1870798, 25.1902855, 36.1288712],
                [4, -4.2805188, -4.2591829, 2.3180082, 2.746881, 6.8290748, 9.2614461, 10.6710271, 16.4520353, 16.6468189, 25.3305449, 25.3437576, 36.22944114],
                [5, -5.800046, -5.7900806, 1.8581875, 2.0994604, 7.4491097, 9.2363277, 11.548832, 16.6482199, 17.0965817, 25.510816, 25.5499717, 36.3588668],
                [6, -7.3688308, -7.363911, 1.2142782, 1.3513812, 7.8700645, 9.1379058, 12.4656007, 16.8446016, 17.688783, 25.7234107, 25.817272, 36.5170667],
                [7, -8.9737425, -8.9712024, 0.4383491, 0.5175454, 8.0866231, 8.9623855, 13.3584213, 17.0266608, 18.4166087, 25.9624472, 26.1561202, 36.7035027],
                [8, -10.6067292, -10.6053681, -0.4359436, -0.3893618, 8.1152388, 8.7099144, 14.1818804, 17.1825278, 19.2527051, 26.2209995, 26.5777533, 36.9172131],
                [9, -12.2624142, -12.2616617, -1.3867016, -1.3588101, 7.9828432, 8.3831192, 14.9036797, 17.303011, 20.1609264, 26.4915472, 27.0918661, 37.156695],
                [10, -13.93698, -13.9365525, -2.3991424, -2.3821582, 7.7173698, 7.9860691, 15.5027844, 17.3813807, 21.1046337, 26.7664264, 27.7037687, 37.4198588],
                [12, -17.332066, -17.3319184, -4.5701329, -4.5635399, 6.8787369, 7.0005668, 16.3015349, 17.3952497, 22.9721275, 27.3000124, 29.208055, 38.0060087],
                [14, -20.7760553, -20.7760004, -6.8934005, -6.8907007, 5.7363123, 5.7926295, 16.5985405, 17.2071153, 24.6505951, 27.7697667, 31.0000508, 38.6484719],
                [16, -24.2586795, -24.2586578, -9.33523671, -9.3341097, 4.3712326, 4.3978962, 16.4868843, 16.8186837, 26.0086783, 28.136359, 32.9308951, 39.3150108],
                [18, -27.7728422, -27.7728332, -11.8732425, -11.8727265, 2.8330567, 2.8459917, 16.0619754, 16.2420804, 26.9877664, 28.3738582, 34.8530587, 39.9723511],
                [20, -31.3133901, -31.3133862, -14.4913014, -14.4910633, 1.15422829, 1.1607057, 15.3958109, 15.4939776, 27.5945782, 28.4682213, 36.6449897, 40.5896641],
                [24, -38.4589732, -38.4589724, -19.9225956, -19.9225403, -2.5397657, -2.5380779, 13.5228427, 13.5527965, 27.8854408, 28.2153594, 39.5125519, 41.6057099],
                [28, -45.6733696, -45.6733694, -25.5617471, -25.5617329, -6.588063, -6.587585, 11.1110798, 11.1206227, 27.2833082, 27.4057488, 41.2349503, 42.2248415],
                [32, -52.942223, -52.9422229, -31.3651544, -31.3651505, -10.9143534, -10.914209, 8.2914962, 8.2946721, 26.0624482, 26.1083526, 41.9535112, 42.3939428],
                [36, -60.2555679, -60.2555679, -37.3026391, -37.302638, -15.4667703, -15.4667243, 5.1456363, 5.1467375, 24.3785094, 24.3960665, 41.9266646, 42.1183561]
              ])
In [18]:
%matplotlib notebook
w1 = widgets.FloatSlider(
    value=45.0, 
    min=0.0, 
    max=90.0, 
    step=5,
    description='tan(-γ/λ) [°] :', 
    readout_format='.0f',
)

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
clr = ("white","grey")
unstable = ax.fill_betweenx(d[:,0], 40, d[:,12], color=clr[1], label="unstable")
stable   = ax.fill_betweenx(d[:,0], 40, 41, color=clr[0], label="stable")

for i in range(d.shape[1]-1, 0, -1):
    ax.fill_betweenx(d[:,0], -10, d[:,i], color=clr[i%2])
    plt.plot(d[:,i], d[:,0], color="black")
    
plt.plot((0, 0), (0, 20), color="black", linestyle='--')

plt.legend(handles=[stable, unstable])
plt.xlabel("λ")
plt.ylabel("-γ")
plt.xlim(-10.0, 36.0)
plt.ylim(0.0, 20.0)
plt.show()

line, = ax.plot((0, 40), (0, 40/2*np.cos(w1.value*np.pi/180)))

def update(w1):
    line.set_ydata((0, 40/2*np.cos(w1*np.pi/180)))
    fig.canvas.draw()
    
widgets.interact(update, w1=w1);