Gelöste Aufgaben/Kw24: Unterschied zwischen den Versionen
Keine Bearbeitungszusammenfassung |
|||
Zeile 19: | Zeile 19: | ||
== Lösung mit Maxima == | == Lösung mit Maxima == | ||
<!--------------------------------------------------------------------------------> | |||
{{MyCodeBlock|title=Header | {{MyCodeBlock|title=Header | ||
|text= | |text= | ||
Gesucht ist die Lösung zu dem nichtlinearen Anfangswertproblem. | |||
|code= | |code= | ||
<syntaxhighlight lang="lisp" line start=1> | <syntaxhighlight lang="lisp" line start=1> | ||
/*********************************************************/ | |||
/* MAXIMA script */ | |||
/* version: wxMaxima 15.08.2 */ | |||
/* author: Andreas Baumgart */ | |||
/* last updated: 2017-09-30 */ | |||
/* ref: Kw24 (TM-C, Labor 6) */ | |||
/* description: finds the solution for */ | |||
/* the nonlinear IVP */ | |||
/*********************************************************/ | |||
</syntaxhighlight> | </syntaxhighlight> | ||
}} | }} | ||
= | |||
<!--------------------------------------------------------------------------------> | |||
{{MyCodeBlock|title=Declarations | |||
|text= | |||
Wir wählen die Längenverhältnisse | Wir wählen die Längenverhältnisse | ||
<math>\displaystyle {{\ell}_{2}}={{\ell}_{1}}, {{\ell}_{3}}=\sqrt{2}\cdot {{\ell}_{1}}</math> | ::<math>\displaystyle {{\ell}_{2}}={{\ell}_{1}}, {{\ell}_{3}}=\sqrt{2}\cdot {{\ell}_{1}}</math> | ||
so dass die Feder für φ''(t) = 45''° entspannt ist. | so dass die Feder für φ''(t) = 45''° entspannt ist. | ||
|code= | |code= | ||
<syntaxhighlight lang="lisp" line start=1> | <syntaxhighlight lang="lisp" line start=1> | ||
1 | /*********************************************************/ | ||
/* declare variational variables - see 6.3 Identifiers */ | |||
declare("δW", alphabetic); /* total virtual work */ | |||
declare( "φ", alphabetic); /* coordinate of rotatin */ | |||
declare("δφ", alphabetic); /* variation of coordinate phi*/ | |||
declare("δw", alphabetic); /* variation of coordinate w */ | |||
/* parameter */ | |||
assume(l[1]>0, k>0, J>0); | |||
params: [l[2] = l[1], l[3] = sqrt(2)*l[1]]; | |||
</syntaxhighlight> | </syntaxhighlight> | ||
}} | }} | ||
= | <!--------------------------------------------------------------------------------> | ||
{{MyCodeBlock|title=Kinematics | |||
|text= | |||
Der Mechanismus hat genau einen Freiheitsgrad, den wir entweder durch die Federlängung ''w'' oder den Drehwinkel ''φ'' erfassen. Den Zusammen zwischen den beiden stellen wir über | Der Mechanismus hat genau einen Freiheitsgrad, den wir entweder durch die Federlängung ''w'' oder den Drehwinkel ''φ'' erfassen. Den Zusammen zwischen den beiden stellen wir über | ||
<math>{{\ell}_{1}^{2}}\cdot {{\mathrm{sin}\left( \mathrm{\varphi}\left( t\right) \right) }^{2}}+{{\left( {{\ell}_{3}}-{{\ell}_{1}}\cdot \mathrm{cos}\left( \mathrm{\varphi}\left( t\right) \right) \right) }^{2}}={{\left( {{\ell}_{2}}+w\right) }^{2}}</math> | ::<math>{{\ell}_{1}^{2}}\cdot {{\mathrm{sin}\left( \mathrm{\varphi}\left( t\right) \right) }^{2}}+{{\left( {{\ell}_{3}}-{{\ell}_{1}}\cdot \mathrm{cos}\left( \mathrm{\varphi}\left( t\right) \right) \right) }^{2}}={{\left( {{\ell}_{2}}+w\right) }^{2}}</math> | ||
her. Nun steht uns frei, ob wir diese Beziehung nach ''w'' oder ''φ'' auflösen. Wir wählen ''φ'' als Minimal-Koordinate und erhalten | her. Nun steht uns frei, ob wir diese Beziehung nach ''w'' oder ''φ'' auflösen. Wir wählen ''φ'' als Minimal-Koordinate und erhalten | ||
<math>w = \displaystyle \frac{{{\ell}_{1}}\cdot \sqrt{12-{{2}^{\frac{7}{2}}}\cdot \cos(\varphi) }-2\cdot {{\ell}_{1}}}{2}</math> | ::<math>w = \displaystyle \frac{{{\ell}_{1}}\cdot \sqrt{12-{{2}^{\frac{7}{2}}}\cdot \cos(\varphi) }-2\cdot {{\ell}_{1}}}{2}</math> | ||
Nach den Regeln der [[Sources/Lexikon/Variationsmethoden|Variationsmethoden]] erhalten wir daraus: | Nach den Regeln der [[Sources/Lexikon/Variationsmethoden|Variationsmethoden]] erhalten wir daraus: | ||
<math>\delta w = \displaystyle \frac{{2^{\frac{3}{2}}}\cdot {{\ell}_{1}}\cdot \sin( \varphi ) \cdot \delta \varphi}{\sqrt{12-{{2}^{\frac{7}{2}}}\cdot \cos(\varphi)}}</math>. | ::<math>\delta w = \displaystyle \frac{{2^{\frac{3}{2}}}\cdot {{\ell}_{1}}\cdot \sin( \varphi ) \cdot \delta \varphi}{\sqrt{12-{{2}^{\frac{7}{2}}}\cdot \cos(\varphi)}}</math>. | ||
|code= | |code= | ||
<syntaxhighlight lang="lisp" line start=1> | <syntaxhighlight lang="lisp" line start=1> | ||
1+1 | /*********************************************************/ | ||
/**** kinematics ***/ | |||
/* coordinates */ | |||
coords : [φ, δφ]; | |||
kin : (l[1]*sin(φ(t)))^2+(l[3]-l[1]*cos(φ(t)))^2=(l[2]+w)^2; | |||
kin : trigsimp(ratsimp(solve(subst(params, kin),w)[2])); | |||
kin : [kin, | |||
δw = subst([epsilon=0], | |||
diff(subst([φ(t)=φ(t)+epsilon*δφ], | |||
subst(kin,w)),epsilon))]; | |||
</syntaxhighlight> | </syntaxhighlight> | ||
}} | }} | ||
= | |||
<!--------------------------------------------------------------------------------> | |||
{{MyCodeBlock|title=Equilibrium Conditions | |||
|text= | |||
Mit dem Prinzip der virtuellen Verrückungen schreiben wir | Mit dem Prinzip der virtuellen Verrückungen schreiben wir | ||
<math>\displaystyle \delta W=-\left( \frac{{{d}^{2}}}{d\,{{t}^{2}}}\cdot \mathrm{\varphi}\left( t\right) \right) \cdot J\cdot \delta \varphi-k\cdot w\cdot \delta w</math> | ::<math>\displaystyle \delta W=-\left( \frac{{{d}^{2}}}{d\,{{t}^{2}}}\cdot \mathrm{\varphi}\left( t\right) \right) \cdot J\cdot \delta \varphi-k\cdot w\cdot \delta w</math> | ||
und erhalten die Bewegungsgleichung | und erhalten die Bewegungsgleichung | ||
<math>\displaystyle J \cdot \ddot{\varphi} + \sqrt{2}\cdot l_1^2 \cdot k \cdot \sin\varphi - \frac{{{2}^{\frac{3}{2}}}\cdot l_1^2 \cdot k \cdot \sin\varphi}{\sqrt{12-{{2}^{\frac{7}{2}}}\cdot \cos\varphi}} = 0</math> | ::<math>\displaystyle J \cdot \ddot{\varphi} + \sqrt{2}\cdot l_1^2 \cdot k \cdot \sin\varphi - \frac{{{2}^{\frac{3}{2}}}\cdot l_1^2 \cdot k \cdot \sin\varphi}{\sqrt{12-{{2}^{\frac{7}{2}}}\cdot \cos\varphi}} = 0</math> | ||
|code= | |code= | ||
<syntaxhighlight lang="lisp" line start=1> | <syntaxhighlight lang="lisp" line start=1> | ||
1 | /*********************************************************/ | ||
/* equilibrium conditions */ | |||
PvV : δW = -J*'diff(φ(t),t,2)*δφ - k*w*δw; | |||
/* dgl for φ ist the coefficient of δφ in δW */ | |||
dgl[1] : expand(coeff(subst([δW = 0], subst(kin, -PvV)),δφ)); | |||
</syntaxhighlight> | </syntaxhighlight> | ||
}} | }} | ||
= | |||
<!--------------------------------------------------------------------------------> | |||
{{MyCodeBlock|title=Equilibrium Points | |||
|text= | |||
Die Bewegungsgleichung machen wir dimensionslos, indem wir charakteristische Größen aus der linearisierten Bewegungsgleichung verwenden. | Die Bewegungsgleichung machen wir dimensionslos, indem wir charakteristische Größen aus der linearisierten Bewegungsgleichung verwenden. | ||
Die drei Gleichgewichts-Lagen bzgl. derer wir die Bewegungsgleichung linearisieren können, sind | Die drei Gleichgewichts-Lagen bzgl. derer wir die Bewegungsgleichung linearisieren können, sind | ||
<math>\begin{array}{ll}\varphi_{0,1} &= \displaystyle -\frac{\pi}{4}\\\varphi_{0,2} &= 0\\\varphi_{0,3} &= \displaystyle +\frac{\pi}{4}\end{array}</math>. | ::<math>\begin{array}{ll}\varphi_{0,1} &= \displaystyle -\frac{\pi}{4}\\\varphi_{0,2} &= 0\\\varphi_{0,3} &= \displaystyle +\frac{\pi}{4}\end{array}</math>. | ||
Eine Taylor-Reihenentwicklung des Rückstellmoments bzgl. φ''<sub>0,3</sub>'' ,also | Eine Taylor-Reihenentwicklung des Rückstellmoments bzgl. φ''<sub>0,3</sub>'' ,also | ||
<math>\varphi = \varphi_{0,3} + \phi</math> | ::<math>\varphi = \varphi_{0,3} + \phi</math> | ||
liefert die linearisierte Bewegungsgleichung | liefert die linearisierte Bewegungsgleichung | ||
<math>\displaystyle\ddot{\phi}+\omega_0^2\cdot\phi = 0 \text{ mit } \omega_0^2 = \frac{{{\ell}_{1}^{2}}\cdot k}{J}</math>. | ::<math>\displaystyle\ddot{\phi}+\omega_0^2\cdot\phi = 0 \text{ mit } \omega_0^2 = \frac{{{\ell}_{1}^{2}}\cdot k}{J}</math>. | ||
Dabei ist ''ω<sub>0</sub>'' die Eigenkreisfreuquenz der linearisierten Bewegungsgleichung. | Dabei ist ''ω<sub>0</sub>'' die Eigenkreisfreuquenz der linearisierten Bewegungsgleichung. | ||
|code= | |code= | ||
<syntaxhighlight lang="lisp" line start=1> | <syntaxhighlight lang="lisp" line start=1> | ||
1 | /*********************************************************/ | ||
/* Zeros of restoring function */ | |||
equilib : solve(subst(['diff(φ(t),t,2)=0],dgl[1]),φ(t)); | |||
/*********************************************************/ | |||
/* derive the eigenfrequency */ | |||
/* of linearized ODE at stable point of equilob. */ | |||
LinODE : taylor (subst([φ(t)=psi],rhs(dgl[1]/J)), psi, %pi/4,1); | |||
ref : omega[0]^2 = coeff(expand(LinODE),psi); | |||
ref : [t[B] = 1/(subst(solve(ref,omega[0])[2],omega[0])/(2*%pi))]; | |||
</syntaxhighlight> | </syntaxhighlight> | ||
}} | }} | ||
= | |||
<!-------------------------------------------------------------------------------->{{MyCodeBlock|title=Make Equations of Motion Dimensionless | |||
|text= | |||
Eine Bezugszeit für unsere nichtlineare Bewegungsgleichung kommt nun aus | Eine Bezugszeit für unsere nichtlineare Bewegungsgleichung kommt nun aus | ||
<math>\displaystyle \omega_0 = \frac{2\pi}{t_{Bez}}</math> | ::<math>\displaystyle \omega_0 = \frac{2\pi}{t_{Bez}}</math> | ||
Damit lautet die dimensionslose Bewegungsgleichung | Damit lautet die dimensionslose Bewegungsgleichung | ||
<math>\displaystyle \frac{d^2}{d\tau^2}\Phi + \frac{{{2}^{\frac{5}{2}}}\cdot {{\pi }^{2}}\cdot \sqrt{12-{{2}^{\frac{7}{2}}}\cdot \mathrm{cos}\left( \Phi\right) }\cdot \mathrm{sin}\left( \Phi\right) -{{2}^{\frac{7}{2}}}\cdot {{\pi }^{2}}\cdot \mathrm{sin}\left( \Phi\right) }{\sqrt{12-{{2}^{\frac{7}{2}}}\cdot \mathrm{cos}\left( \Phi\right) }} = 0 \text{ mit } \tau = t/t_{Bez}</math> | ::<math>\displaystyle \frac{d^2}{d\tau^2}\Phi + \frac{{{2}^{\frac{5}{2}}}\cdot {{\pi }^{2}}\cdot \sqrt{12-{{2}^{\frac{7}{2}}}\cdot \mathrm{cos}\left( \Phi\right) }\cdot \mathrm{sin}\left( \Phi\right) -{{2}^{\frac{7}{2}}}\cdot {{\pi }^{2}}\cdot \mathrm{sin}\left( \Phi\right) }{\sqrt{12-{{2}^{\frac{7}{2}}}\cdot \mathrm{cos}\left( \Phi\right) }} = 0 \text{ mit } \tau = t/t_{Bez}</math> | ||
Nullstellen des Rückstellmoments zeigen die Gleichgewichtslagen des Systems an. Die Steigung dieses Rückstellmoments in den Gleichgewichtslagen deutet dabei auf die Stabilität der linearisierten Bewegung (oder Bewegung im Kleinen) hin. So ist die Lösung ''Φ'' = 0 eine instabile Gleichgewichtslage, die beiden anderen sind stabil. | [[Datei:Kw24-11.png|mini|Rückstellmoment der Feder|alternativtext=]]Den nichtlinearen (zweiten) Summand können wir als dimensionsloses "Rückstellmoment auffassen und über dem Drehwinkel auftragen: | ||
Nullstellen des Rückstellmoments zeigen die Gleichgewichtslagen des Systems an. Die Steigung dieses Rückstellmoments in den Gleichgewichtslagen deutet dabei auf die Stabilität der linearisierten Bewegung (oder Bewegung im Kleinen) hin. So ist die Lösung ''Φ'' = 0 eine instabile Gleichgewichtslage, die beiden anderen sind stabil. | |||
|code= | |code= | ||
<syntaxhighlight lang="lisp" line start=1> | <syntaxhighlight lang="lisp" line start=1> | ||
1 | /*********************************************************/ | ||
/* make dgl dimensionless employing the eigenfrequency */ | |||
dgl[2] : subst([diff(φ(t),t,2) = 'diff(Phi,t,2)/t[B]^2, | |||
φ(t) = Phi ] , dgl[1]); | |||
dgl[2] : solve(subst(ref,dgl[2]), 'diff(Phi,t,2)); | |||
/* equilibrium angles für Φ */ | |||
plot2d(-rhs(dgl[2][1]), [Phi,-1,1], | |||
[legend, "Rückstell'moment'"], | |||
[xlabel, "Φ/1->"], [ylabel, "M/1->"]); | |||
</syntaxhighlight> | </syntaxhighlight> | ||
}} | }} | ||
<math>\ddot{\Phi} = f(\Phi)</math> | <!--------------------------------------------------------------------------------> | ||
{{MyCodeBlock|title=Numerical Solution | |||
|text=Für die numerische Lösung müssen wir die Bewegungsgleichung | |||
::<math>\ddot{\Phi} = f(\Phi)</math> | |||
umschreiben auf eine Differentialgleichung erster Ordnung: | umschreiben auf eine Differentialgleichung erster Ordnung: | ||
<math>\left(\begin{array}{c}\Phi\\\Omega\end{array}\right) = \left(\begin{array}{c}\Omega\\f(\Phi)\end{array}\right)</math> | ::<math>\left(\begin{array}{c}\Phi\\\Omega\end{array}\right) = \left(\begin{array}{c}\Omega\\f(\Phi)\end{array}\right)</math> | ||
Diese lösen wir mit einer fertigen Routine - hier einem Runge-Kutta-Verfahren - als Anfangswertproblem. Anders als bei linearisierten Bewegungsgleichungen hängt die Charakteristik der Lösung von den Anfangsbedingungen ab. | Diese lösen wir mit einer fertigen Routine - hier einem Runge-Kutta-Verfahren - als Anfangswertproblem. Anders als bei linearisierten Bewegungsgleichungen hängt die Charakteristik der Lösung von den Anfangsbedingungen ab. | ||
|code= | |code= | ||
<syntaxhighlight lang="lisp" line start=1> | <syntaxhighlight lang="lisp" line start=1> | ||
1 | /*********************************************************/ | ||
/* numerical solution of IVP */ | |||
fpprintprec:3$ | |||
lhs(dgl[2][1]) = float(expand(rhs(dgl[2][1]))); | |||
times : subst([t0 = 0, tmax = 6, dt = 0.01], | |||
[t, t0, tmax, dt]); | |||
dgl1stOrder : [Omega,float(expand(rhs(dgl[2][1])))]; | |||
stateVabs : [Phi,Omega]; | |||
initiVals : [[1,0],[-1,0], | |||
[0,0.001],[0,3], | |||
[0,8], [0,10]]; | |||
/* solution of IVP (ivs) */ | |||
for i: 1 thru length(initiVals) do | |||
ivs[i] : rk(dgl1stOrder, stateVabs, initiVals[i], times)$ | |||
</syntaxhighlight> | </syntaxhighlight> | ||
}} | }} | ||
Zeile 153: | Zeile 215: | ||
Das schauen wir uns am Verlauf von Lösungen im Zeitbereich an: | Das schauen wir uns am Verlauf von Lösungen im Zeitbereich an: | ||
<table class="wikitable" style="background-color:white; float: left; margin-right:14px;"> | |||
<tr><td>a)[[Datei:K24-21.png|rahmenlos]] | |||
</td><td>[[Datei:Kw24-24.png|rahmenlos]] | |||
b) | |||
</td></tr><tr><td>[[Datei:Kw24-22.png|rahmenlos]] | |||
c) | |||
</td><td>[[Datei:Kw24-25.png|rahmenlos]] | |||
d) | |||
</td></tr><tr><td>[[Datei:Kw24-23.png|rahmenlos]] | |||
e) | |||
</td><td>[[Datei:Kw24-26.png|rahmenlos]] | |||
f) | |||
</td></tr></table> | |||
Lösung f) fällt aus dem Raster - hier gibt es keine periodische Lösung mehr: Bei einer anderen Skalierung der Φ-Achse erkennen wir, dass der Stab AC um den Punkt A rotiert:[[Datei:Kw24-30.png|mini|Lösung f) im Zeitbereich.|alternativtext=|ohne]] | Lösung f) fällt aus dem Raster - hier gibt es keine periodische Lösung mehr: Bei einer anderen Skalierung der Φ-Achse erkennen wir, dass der Stab AC um den Punkt A rotiert:[[Datei:Kw24-30.png|mini|Lösung f) im Zeitbereich.|alternativtext=|ohne]] | ||
Zeile 164: | Zeile 238: | ||
Analog sehen Sie in der Lösung c) (grün) eine Separatrix zwischen der periodischen Lösung der kleinen Schwingung um die Gleichgewichtslage und solchen periodischen Lösungen, bei denen der Stab die Horizontale durchbricht: [[Datei:Kw24-32.png|mini|Phasendiagramme|alternativtext=|ohne]] | Analog sehen Sie in der Lösung c) (grün) eine Separatrix zwischen der periodischen Lösung der kleinen Schwingung um die Gleichgewichtslage und solchen periodischen Lösungen, bei denen der Stab die Horizontale durchbricht: [[Datei:Kw24-32.png|mini|Phasendiagramme|alternativtext=|ohne]] | ||
<!--------------------------------------------------------------------------------> | <!--------------------------------------------------------------------------------> | ||
{{MyCodeBlock|title=Post-Processing | {{MyCodeBlock|title=Post-Processing | ||
Zeile 170: | Zeile 243: | ||
|code= | |code= | ||
<syntaxhighlight lang="lisp" line start=1> | <syntaxhighlight lang="lisp" line start=1> | ||
1 | |||
/*********************************************************/ | |||
/* plot results */ | |||
for i: 1 thru length(initiVals) do | |||
(T[i] : makelist(ivs[i][j][1],j,1,length(ivs[i])), | |||
P[i] : makelist(ivs[i][j][2],j,1,length(ivs[i])), | |||
O[i] : makelist(ivs[i][j][3],j,1,length(ivs[i])), | |||
plot2d([discrete, T[i], P[i]], [y,-%pi,%pi], | |||
[title, sconcat("start @: ",string(initiVals[i]))], | |||
[xlabel,"τ/1->"], [ylabel,"Φ/1->"])); | |||
/* phase plot */ | |||
curves : makelist([discrete,P[i],O[i]],i,1,length(initiVals))$ | |||
plot2d(curves, [legend, false], [x,-2,10], | |||
[ylabel,"Ω/1->"], [xlabel,"Φ/1->"]); | |||
plot2d(curves, [legend, false], [x,-2,2], [y,-5,5], | |||
[ylabel,"Ω/1->"], [xlabel,"Φ/1->"]); | |||
plot2d(curves, [legend, false], [x,-2,2], [y,-5,5], [svg_file, "Kw24-20.svg"], | |||
[ylabel,"Ω/1->"], [xlabel,"Φ/1->"]); | |||
</syntaxhighlight> | </syntaxhighlight> | ||
}} | }} | ||
<hr/> | <hr/> |
Version vom 29. März 2021, 06:02 Uhr
Aufgabenstellung
Ein Mechanismus besteht aus einem starren Stab AC (Länge ℓ1, Massenmoment J) und einer Feder BC (Länge ℓ2 im entspannten Zustand, Steifigkeit k).
Die Körper sind drehbar gelagert, der Abstand AB ist ℓ3.
Hier geht es um die verschiedenen Lösungstypen des Systems - in Abhängigkeit von den Anfangsbedingungen. Gesucht sind Lösungen als Anfangswertproblem, nichtlineare Schwingungen und Aussagen über die Stabilität der Bewegungen.
Lösung mit Maxima
Header
Gesucht ist die Lösung zu dem nichtlinearen Anfangswertproblem.
/*********************************************************/
/* MAXIMA script */
/* version: wxMaxima 15.08.2 */
/* author: Andreas Baumgart */
/* last updated: 2017-09-30 */
/* ref: Kw24 (TM-C, Labor 6) */
/* description: finds the solution for */
/* the nonlinear IVP */
/*********************************************************/
Declarations
Wir wählen die Längenverhältnisse
so dass die Feder für φ(t) = 45° entspannt ist.
/*********************************************************/
/* declare variational variables - see 6.3 Identifiers */
declare("δW", alphabetic); /* total virtual work */
declare( "φ", alphabetic); /* coordinate of rotatin */
declare("δφ", alphabetic); /* variation of coordinate phi*/
declare("δw", alphabetic); /* variation of coordinate w */
/* parameter */
assume(l[1]>0, k>0, J>0);
params: [l[2] = l[1], l[3] = sqrt(2)*l[1]];
Kinematics
Der Mechanismus hat genau einen Freiheitsgrad, den wir entweder durch die Federlängung w oder den Drehwinkel φ erfassen. Den Zusammen zwischen den beiden stellen wir über
her. Nun steht uns frei, ob wir diese Beziehung nach w oder φ auflösen. Wir wählen φ als Minimal-Koordinate und erhalten
Nach den Regeln der Variationsmethoden erhalten wir daraus:
- .
/*********************************************************/
/**** kinematics ***/
/* coordinates */
coords : [φ, δφ];
kin : (l[1]*sin(φ(t)))^2+(l[3]-l[1]*cos(φ(t)))^2=(l[2]+w)^2;
kin : trigsimp(ratsimp(solve(subst(params, kin),w)[2]));
kin : [kin,
δw = subst([epsilon=0],
diff(subst([φ(t)=φ(t)+epsilon*δφ],
subst(kin,w)),epsilon))];
Equilibrium Conditions
Mit dem Prinzip der virtuellen Verrückungen schreiben wir
und erhalten die Bewegungsgleichung
/*********************************************************/
/* equilibrium conditions */
PvV : δW = -J*'diff(φ(t),t,2)*δφ - k*w*δw;
/* dgl for φ ist the coefficient of δφ in δW */
dgl[1] : expand(coeff(subst([δW = 0], subst(kin, -PvV)),δφ));
Equilibrium Points
Die Bewegungsgleichung machen wir dimensionslos, indem wir charakteristische Größen aus der linearisierten Bewegungsgleichung verwenden.
Die drei Gleichgewichts-Lagen bzgl. derer wir die Bewegungsgleichung linearisieren können, sind
- .
Eine Taylor-Reihenentwicklung des Rückstellmoments bzgl. φ0,3 ,also
liefert die linearisierte Bewegungsgleichung
- .
Dabei ist ω0 die Eigenkreisfreuquenz der linearisierten Bewegungsgleichung.
/*********************************************************/
/* Zeros of restoring function */
equilib : solve(subst(['diff(φ(t),t,2)=0],dgl[1]),φ(t));
/*********************************************************/
/* derive the eigenfrequency */
/* of linearized ODE at stable point of equilob. */
LinODE : taylor (subst([φ(t)=psi],rhs(dgl[1]/J)), psi, %pi/4,1);
ref : omega[0]^2 = coeff(expand(LinODE),psi);
ref : [t[B] = 1/(subst(solve(ref,omega[0])[2],omega[0])/(2*%pi))];
Make Equations of Motion Dimensionless
Eine Bezugszeit für unsere nichtlineare Bewegungsgleichung kommt nun aus
Damit lautet die dimensionslose Bewegungsgleichung
Den nichtlinearen (zweiten) Summand können wir als dimensionsloses "Rückstellmoment auffassen und über dem Drehwinkel auftragen:
Nullstellen des Rückstellmoments zeigen die Gleichgewichtslagen des Systems an. Die Steigung dieses Rückstellmoments in den Gleichgewichtslagen deutet dabei auf die Stabilität der linearisierten Bewegung (oder Bewegung im Kleinen) hin. So ist die Lösung Φ = 0 eine instabile Gleichgewichtslage, die beiden anderen sind stabil.
/*********************************************************/
/* make dgl dimensionless employing the eigenfrequency */
dgl[2] : subst([diff(φ(t),t,2) = 'diff(Phi,t,2)/t[B]^2,
φ(t) = Phi ] , dgl[1]);
dgl[2] : solve(subst(ref,dgl[2]), 'diff(Phi,t,2));
/* equilibrium angles für Φ */
plot2d(-rhs(dgl[2][1]), [Phi,-1,1],
[legend, "Rückstell'moment'"],
[xlabel, "Φ/1->"], [ylabel, "M/1->"]);
Numerical Solution
Für die numerische Lösung müssen wir die Bewegungsgleichung
umschreiben auf eine Differentialgleichung erster Ordnung:
Diese lösen wir mit einer fertigen Routine - hier einem Runge-Kutta-Verfahren - als Anfangswertproblem. Anders als bei linearisierten Bewegungsgleichungen hängt die Charakteristik der Lösung von den Anfangsbedingungen ab.
/*********************************************************/
/* numerical solution of IVP */
fpprintprec:3$
lhs(dgl[2][1]) = float(expand(rhs(dgl[2][1])));
times : subst([t0 = 0, tmax = 6, dt = 0.01],
[t, t0, tmax, dt]);
dgl1stOrder : [Omega,float(expand(rhs(dgl[2][1])))];
stateVabs : [Phi,Omega];
initiVals : [[1,0],[-1,0],
[0,0.001],[0,3],
[0,8], [0,10]];
/* solution of IVP (ivs) */
for i: 1 thru length(initiVals) do
ivs[i] : rk(dgl1stOrder, stateVabs, initiVals[i], times)$
tmp
Das schauen wir uns am Verlauf von Lösungen im Zeitbereich an:
a) |
b) |
c) |
d) |
e) |
f) |
Lösung f) fällt aus dem Raster - hier gibt es keine periodische Lösung mehr: Bei einer anderen Skalierung der Φ-Achse erkennen wir, dass der Stab AC um den Punkt A rotiert:
Oft erkennt man das Charakteristische einer Lösung im Phasendiagramm. Hier wird die Geschwindigkeit Ω über dem Drehwinkel Φ aufgetragen:
Periodische Lösungen - die oft gesucht sind - erkennt man an einem geschlossenen Verlauf. Die hellblaue Lösung gehört zu den Anfangswerten von f), bei dem der Drehwinkel immer größer wird - sich also keine periodische Lösung ergibt.
Wir sehen, dass die Lösung e) (schwarz) eine "gerade noch" periodische Lösung ist - etwas mehr Anfangsgeschwindigkeit im System und die Lösung sieht wie in f) aus. Diese Lösung nennt man auch Separatrix - sie separiert zwei Charakteristiken von Lösungen.
Analog sehen Sie in der Lösung c) (grün) eine Separatrix zwischen der periodischen Lösung der kleinen Schwingung um die Gleichgewichtslage und solchen periodischen Lösungen, bei denen der Stab die Horizontale durchbricht:
Post-Processing
Text
/*********************************************************/
/* plot results */
for i: 1 thru length(initiVals) do
(T[i] : makelist(ivs[i][j][1],j,1,length(ivs[i])),
P[i] : makelist(ivs[i][j][2],j,1,length(ivs[i])),
O[i] : makelist(ivs[i][j][3],j,1,length(ivs[i])),
plot2d([discrete, T[i], P[i]], [y,-%pi,%pi],
[title, sconcat("start @: ",string(initiVals[i]))],
[xlabel,"τ/1->"], [ylabel,"Φ/1->"]));
/* phase plot */
curves : makelist([discrete,P[i],O[i]],i,1,length(initiVals))$
plot2d(curves, [legend, false], [x,-2,10],
[ylabel,"Ω/1->"], [xlabel,"Φ/1->"]);
plot2d(curves, [legend, false], [x,-2,2], [y,-5,5],
[ylabel,"Ω/1->"], [xlabel,"Φ/1->"]);
plot2d(curves, [legend, false], [x,-2,2], [y,-5,5], [svg_file, "Kw24-20.svg"],
[ylabel,"Ω/1->"], [xlabel,"Φ/1->"]);
Links
- ...
Literature
- ...