Gelöste Aufgaben/Kv52

Aufgabenstellung

Lageplan

Ein starrer Stab AB (Massenmoment J_A, Länge ℓ₁) wird aus dem Winkel φ₀ im Erdschwerefeld losgelassen und stößt in C auf einen Anschlag. Der Stoß zwischen Stab und Oberfläche sei ideal-elastisch.

Gesucht ist die nichtlineare Bewegungsgleichung und die numerische Lösung als Anfangswertproblem. Dabei denken wir uns den Anschlag als elastische Feder, die nur für φ<-π/2 Kontakt zum Stab hat. Die Federkraft ist also Null, solange der Stab die Oberfläche in C nicht berührt und sie ist proportional zur Federkompression w, wenn sich Kugel und Oberfläche berühren.

Lösung mit Maxima

Header

Wir stellen die Bewegungsgleichungen des Systems als System con Differentialgleichungen erster Ordnung auf. Die Nichtlinearität kommt aus den großen Winkeln φ(t) und dem Kontakt mit der Wand.

Mit unterschiedlichen Steifigkeiten für den Kontakt testen wir die Möglichkeiten der numerischen Integration aus.

/*********************************************************/
/* MAXIMA script                                         */
/* version: wxMaxima 15.08.2                             */
/* author: Andreas Baumgart                              */
/* last updated: 2018-05-14                              */
/* ref: Kv52 (TM-C, Labor 5)                             */
/* description: finds the solution for                   */
/*              the nonlinear IVP                        */
/*********************************************************/

Declarations

Die System-Parameter sind

$\begin{array}{l}{\ell }_2=\frac{3}{4}{\ell }_1\\ J_A=\frac{1}{3}m{\ell }_1^2\end{array}$.

Zum Dimensionslos-Machen der Bewegungsgleichungen brauchen wir später eine Bezugszeit t_B, die wir mit Hilfe der Eigenfrequenz der zugeordneten linearen Systems so wählen:

$t_B=1/{\omega }_0$

(Achtung: das macht 2π-periodische Lösungen)

Die Eigenkreisfrequenz bei Schwingungen um φ =- π/2 herum ist

${\omega }_0^2=\frac{mg\ell }{2J_A}$

Außerdem wählen wir eine dimensionslose Federsteifigkeit κ, so dass

$k=\kappa \cdot \frac{mg}{{\ell }_2}$

Kennlinie

Für den nichtlinearen Kontakt wählen wir eine Kennlinie wie in Kw23 zu

$C(\tilde{u},\epsilon )=\{\begin{array}{ll}0& \text{ wenn }\tilde{u}<-\epsilon \\ \frac{1}{2}(\tilde{u}+\epsilon {)}^2& \text{ wenn }-\epsilon <\tilde{u}<\epsilon \text{ und }\\ \tilde{u}& \text{ sonst }\end{array}$

Und so wie rechts im Bild sieht sie dann aus:

/* declarations */
assume(g>0, l[1]>0, l[2]>0, m>0);

/* parameter */
params: [k    = kappa*m*g/l[2],   /*spring stiffness*/
         l[2] = 3/4*l[1],
         t[B] = 1/omega[0]     ,  /*reference time*/
         omega[0] = sqrt((m*g*l[1]/2)/J[A]),
         J[A] = m*l[1]^2/3,
         kappa= 1000];            /*dim'less spring stiffness*/

/**** define nonlinear spring characteristic ***/
C(u,epsilon) := if u <- epsilon then
                   0
                elseif u< epsilon then
                   1/2*(u+epsilon)^2
                else
                u;         
plot2d(K(u,0.5),[u,-1,1],
          [ylabel,"K/1->"], [xlabel,"u/1->"],
          [legend, "contact characteristic, ε=0.5"]);

Equilibrium Conditions

Die Gleichgewichtsbeziehungen erhalten wir aus einem Momentengleichgewicht um den Punkt A.

Das Freikörperbild

liefert die Bewegungsgleichung

$J_A\ddot{\phi }-K(\phi ){\ell }_2+mg\frac{{\ell }_1}{2}\cos (\phi )=0$

Die Kontaktkraft konstruieren wir mit Hilfe der nichtlinearen Kennlinie C und

$\tilde{u}=-\frac{\pi }{2}-\phi$

zu

$K(\phi )={\ell }_2\cdot C(-\pi /2-\phi ,\epsilon )$.

Mit der dimensionslosen Zeit

$\tau =t/t_B$

finden wir als Bewegungsgleichung des Systems

$\frac{d^2}{d{\tau }^2}\phi +\cos (\phi )-\frac{3}{2}\kappa \cdot C(-\pi /2-\phi ,\epsilon )=0$.

/********************************/
/* define ODE in dim'less coordinate h */
dgl : J[A]*'diff(phi,t,2)-k*l[2]^2*C(-%pi/2-phi,1/100)+m*g*l[1]/2*cos(phi) = 0;
dgl : subst(['diff(phi,t,2)='diff(phi,tau,2)/t[B]^2],dgl);
dgl : expand(solve(dgl, 'diff(phi,tau,2)))[1];
dgl : expand(subst(params,dgl));

Solving

Umschreiben in ein Differentialgleichungssystem erster Ordnung - wie wir es für die numerische Lösung brauchen - liefert

$\left(\begin{array}{c}\dot{\phi }\\ \dot{\omega }\end{array}\right)=\left(\begin{array}{c}\omega \\ \frac{3}{2}\kappa \cdot C(-\pi /2-\phi ,\epsilon )-\cos (\phi )\end{array}\right)$.

Mit den Anfangsbedingungen

$\begin{array}{l}\phi (0)=0\\ \omega (0)=0\end{array}$

und einer Lösungsroutine nach dem Runge-Kutta-Verfahren 4.ter Ordnung liefert das numerische Ergebnis.

/********************************/
/* numerical solution of IVP */
times : subst([t0 = 0, tmax = 10, dt = 0.01],
                    [t, t0, tmax, dt]);
dgl1stOrder : subst(params,[omega,float(expand(rhs(dgl)))]);
stateVabs : [phi,omega];
initiVals : [0,0];

ivs : rk(dgl1stOrder, stateVabs, initiVals, times)$

Post-Processing

Wir tragen hier die Ergebnisse zweier Simulationsrechnungen auf:

... in dieser Spalte den Referenzfall für $\kappa =100.$	... und hier für $\kappa =1000.$
zunächst im Zeitbereich
und dann im Phasenraum:

Im Referenzfall erhalten wir - wie erhofft - eine periodische Lösung: die Phasenkurve ist geschlossen. Schon für eine um den Faktur 10 größere Steifigkeit des Anschlags jedoch weicht die Lösung sichtbar von einer periodischen ab. Das klassische Problem von numerischen Integrationsroutinen wird hier sichtbar.

/********************************/
/* plot time functions */
T : makelist(ivs[j][1],j,1,length(ivs))$
P : makelist(ivs[j][2],j,1,length(ivs))$
O : makelist(ivs[j][3],j,1,length(ivs))$
plot2d([discrete, T, P],
[title, sconcat("start @: ",string(initiVals))],
[xlabel,"τ/1->"], [ylabel,"φ/1->"]);

/* phase plot */
plot2d([discrete,P,O], [legend, false], /*[x,-1,10],*/
[ylabel,"φ'/1->"], [xlabel,"φ/1->"]);

---

Links

• ...

Literature

• ...