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

                        

Nichtlineare Schwingungsdynamik

Übung - Aufgabe 8

Bestimmen Sie für die Differentialgleichung

\begin{equation*} \ddot{q} + \omega_0^2 (q+\beta^2q^2) = P \cos(\Omega t) \end{equation*}

1) eine harmonische Lösung unter Nutzung der Störungsrechnung. Wie verhalten sich die Amplituden im Resonanzfall bei Abbruch nach dem ersten Glied bzw. nach dem zweiten Glied?

2) eine subharmonische Lösung mit der Periode $T = 4\pi/\Omega$

Aufgabenteil 1)

Überführen der DGL in die Form

\begin{align*} \eta^2 x^{\prime \prime} + x + \mu (x^2 - P_0 \cos(\tau)) = 0 \end{align*}

Potenzreihen:

\begin{align*} x(\tau) &= x_0(\tau) + \mu x_1(\tau) + \mu^2 x_2(\tau) \\ \Omega(\tau) &= \Omega_0(\tau) + \mu \Omega_1(\tau) + \mu^2 \Omega_2(\tau) \\ \text{bzw. }\; \eta(\tau) &= 1 + \mu \eta_1(\tau) + \mu^2 \eta_2(\tau) \\ \end{align*}

Aus den Lösungen der Koeffizientenvergleiche in $\mu^0, \mu^1, \mu^2$ unter der Forderung nach periodischen Lösungen und einer Anfangsbedingung mit $A_1 = 0$ folgt

\begin{align*} \eta = 1 - \mu \frac{P_0}{2 A_0} - \mu^2 \left( \frac{5}{12} A_0^2 + \frac{P_0^2}{8A_0^2} \right) + \mathcal{O}(\mu^3) \end{align*}
In [3]:
import matplotlib.pyplot as plt
import numpy as np
A = np.hstack((np.linspace(-1.0, -0.001, 100),  np.NaN,  np.linspace(0.001, 1.0, 100)))
A_abs = np.abs(A)

mu = 0.2
P0 = 2E-4

# Abbruch nach dem 1. Glied
eta1 = 1.0-mu*P0/2/A
# Abbruch nach dem 2. Glied
eta2 = 1.0-mu*P0/2/A-mu**2*(5/12*A**2+P0**2/8/A**2)

plt.plot(eta1, A_abs)
plt.plot(eta2, A_abs)

plt.xlabel("η")
plt.ylabel("|A|")
plt.legend(('Abbruch nach dem 1. Glied', 'Abbruch nach dem 2. Glied'))
plt.xlim(0.98, 1.02)
plt.show()

Aufgabenteil 2)

In [4]:
from sympy import *
init_printing()
In [5]:
## variables
t = symbols('t', real=True)
## functions
q = Function('q')(t)
## constants
om, Om, beta, P = symbols('omega_0 Omega beta P', real=True)
q0, A, B = symbols('q_0 A B', real=True)
In [6]:
dgl = Eq(Derivative(q, t, t) +  om**2*(q+beta**2*q**2) - P*cos(Om*t),0)
dgl
Out[6]:
$$- P \cos{\left (\Omega t \right )} + \omega_{0}^{2} \left(\beta^{2} q^{2}{\left (t \right )} + q{\left (t \right )}\right) + \frac{d^{2}}{d t^{2}} q{\left (t \right )} = 0$$

Ansatz: $q(t) = q_0 + A\cos(\frac{\Omega}{2} t) + B\sin(\frac{\Omega}{2}t)$

einsetzen:

In [7]:
dgl = dgl.replace(q,q0+A*cos(Om/2*t)+B*sin(Om/2*t))
dgl
Out[7]:
$$- P \cos{\left (\Omega t \right )} + \omega_{0}^{2} \left(A \cos{\left (\frac{\Omega t}{2} \right )} + B \sin{\left (\frac{\Omega t}{2} \right )} + \beta^{2} \left(A \cos{\left (\frac{\Omega t}{2} \right )} + B \sin{\left (\frac{\Omega t}{2} \right )} + q_{0}\right)^{2} + q_{0}\right) + \frac{\partial^{2}}{\partial t^{2}} \left(A \cos{\left (\frac{\Omega t}{2} \right )} + B \sin{\left (\frac{\Omega t}{2} \right )} + q_{0}\right) = 0$$

Ableitung ausführen und Terme zusammenfassen:

In [8]:
dgl = dgl.expand().doit().simplify().trigsimp()
# sin^2(x) bzw. cos^2(x) ersetzen ...
dgl = dgl.rewrite(sin, exp).rewrite(cos, exp).expand().rewrite(exp, sin).simplify()
dgl
Out[8]:
$$\frac{A^{2} \beta^{2} \omega_{0}^{2} \cos{\left (\Omega t \right )}}{2} + \frac{A^{2} \beta^{2} \omega_{0}^{2}}{2} + A B \beta^{2} \omega_{0}^{2} \sin{\left (\Omega t \right )} - \frac{A \Omega^{2} \cos{\left (\frac{\Omega t}{2} \right )}}{4} + 2 A \beta^{2} \omega_{0}^{2} q_{0} \cos{\left (\frac{\Omega t}{2} \right )} + A \omega_{0}^{2} \cos{\left (\frac{\Omega t}{2} \right )} - \frac{B^{2} \beta^{2} \omega_{0}^{2} \cos{\left (\Omega t \right )}}{2} + \frac{B^{2} \beta^{2} \omega_{0}^{2}}{2} - \frac{B \Omega^{2} \sin{\left (\frac{\Omega t}{2} \right )}}{4} + 2 B \beta^{2} \omega_{0}^{2} q_{0} \sin{\left (\frac{\Omega t}{2} \right )} + B \omega_{0}^{2} \sin{\left (\frac{\Omega t}{2} \right )} - P \cos{\left (\Omega t \right )} + \beta^{2} \omega_{0}^{2} q_{0}^{2} + \omega_{0}^{2} q_{0} = 0$$

Koeffizientenvergleich in den Skalaren (unabhängig von $t$):

In [9]:
KV_skal = Eq(dgl.lhs.as_independent(t)[0],dgl.rhs.as_independent(t)[0])
KV_skal
Out[9]:
$$\frac{A^{2} \beta^{2} \omega_{0}^{2}}{2} + \frac{B^{2} \beta^{2} \omega_{0}^{2}}{2} + \beta^{2} \omega_{0}^{2} q_{0}^{2} + \omega_{0}^{2} q_{0} = 0$$

Koeffizientenvergleich in $\cos(\frac{\Omega}{2}t)$:

In [10]:
KV_cos1_2 = Eq(dgl.lhs.coeff(cos(Om/2*t)),dgl.rhs.coeff(cos(Om/2*t)))
KV_cos1_2
Out[10]:
$$- \frac{A \Omega^{2}}{4} + 2 A \beta^{2} \omega_{0}^{2} q_{0} + A \omega_{0}^{2} = 0$$

Koeffizientenvergleich in $\cos(\Omega t)$:

In [11]:
KV_cos = Eq(dgl.lhs.coeff(cos(Om*t)),dgl.rhs.coeff(cos(Om*t)))
KV_cos
Out[11]:
$$\frac{A^{2} \beta^{2} \omega_{0}^{2}}{2} - \frac{B^{2} \beta^{2} \omega_{0}^{2}}{2} - P = 0$$

Koeffizientenvergleich in $sin(\frac{\Omega}{2}t)$:

In [12]:
KV_sin1_2 = Eq(dgl.lhs.coeff(sin(Om/2*t)),dgl.rhs.coeff(sin(Om/2*t)))
KV_sin1_2
Out[12]:
$$- \frac{B \Omega^{2}}{4} + 2 B \beta^{2} \omega_{0}^{2} q_{0} + B \omega_{0}^{2} = 0$$

Koeffizientenvergleich in $sin(\Omega t)$:

In [13]:
KV_sin = Eq(dgl.lhs.coeff(sin(Om*t)),dgl.rhs.coeff(sin(Om*t)))
KV_sin
Out[13]:
$$A B \beta^{2} \omega_{0}^{2} = 0$$

Aus dem Koeffizientenvergleich in $sin(\Omega t)$ folgt $A = 0 \lor B=0$. In Kombination mit dem Koeffizientenvergleich in $\cos(\Omega t)$ mit $\beta^2 > 0$, $\omega_0^2>0$, $A^2>0$, $B^2>0$ folgt $B = 0$.

Die gesuchte Lösung $A(P,\Omega)$ aus dem Koeffizientenvergleich in $\cos(\Omega t)$ :

In [14]:
Ares = solve(KV_cos.replace(B,0),A)
Ares[0] = Ares[0].replace(om,Om/2)
Ares[1] = Ares[1].replace(om,Om/2)
Eq(A,Ares[1].replace(om,Om/2))
Out[14]:
$$A = \frac{2 \sqrt{2} \sqrt{P}}{\Omega \beta}$$

$q_0$ kann anschließend aus dem skalaren Koeffizientenvergleich ermittelt werden:

In [15]:
q0res = solve(KV_skal.replace(B,0).replace(A,Ares[1]),q0)
q0res
Out[15]:
$$\left [ \frac{- \Omega + \sqrt{\Omega^{2} - 16 P \beta^{2}}}{2 \Omega \beta^{2}}, \quad - \frac{\Omega + \sqrt{\Omega^{2} - 16 P \beta^{2}}}{2 \Omega \beta^{2}}\right ]$$

$q_0 = -\frac{1}{2\beta^2} \pm \sqrt{\frac{1}{4\beta^4}-\frac{4P}{\Omega^2\beta^4}}$

Zuletzt folgt aus dem Koeffizientenvergleich in $\cos(\Omega/2 t)$:

In [16]:
KV_cos1_2_1 = KV_cos1_2.replace(A,Ares[1].replace(om,Om/2)).replace(q0,q0res[0]).replace(om,Om/2).expand().doit().simplify()
KV_cos1_2_1
Out[16]:
$$\frac{\sqrt{2} \sqrt{P} \left(- \Omega + \sqrt{\Omega^{2} - 16 P \beta^{2}}\right)}{2 \beta} = 0$$

Mit $P \neq 0$, $\beta \neq 0$ folgt eine implizite Darstellung für $\Omega$.

In [17]:
KV_cos1_2_2 = Eq(KV_cos1_2_1.lhs*2*beta/sqrt(2)/sqrt(P)+Om,KV_cos1_2_1.rhs*2*beta/sqrt(2)/sqrt(P)+Om)
KV_cos1_2_2
Out[17]:
$$\sqrt{\Omega^{2} - 16 P \beta^{2}} = \Omega$$

Die Lösung existiert nur für $16 P \beta^2 < \Omega^2$ bzw. $4 P \beta^2 < \omega_0^2$