Gelöste Aufgaben/Kw23

Aufgabenstellung

Eine Kugel (Masse m, Radius r) wird aus der Höhe H im Erdschwerefeld losgelassen, kommt auf den Boden auf und „springt“ dann wie ein Flummi auf und ab. Der Stoß zwischen Kugel und Oberfläche sei ideal-elastisch. Gefragt ist eine numerische Lösung des Problems als Anfangswertproblem.

Dabei modellieren Sie die Kugel als elastisch, die Unterlage als starr. Die Kugel können Sie sich so wie unten skizziert im unteren Teil durch eine Feder ersetzt denken.

Lageplan

Eine elastische Kugel wird im Erdschwerefeld losgelassen und „springt“ wie ein Flummi auf und ab. Gesucht ist die numerische Lösung als Anfangswertproblem.

Lösung mit Maxima

Die Schwierigkeit kommt aus der Modellierung der Kontaktkraft zwischen Boden und Kugel. Oft reicht es, ein phänomenologisches Modell zu implementieren - also die Kontaktkraft als reine "Federkraft" zu interpretieren.

Kugel-Modell: elastischer Kontakt mit Einfederung w.

Die Federkraft K ist also Null, solange die Kugel die Oberfläche nicht berührt und sie sei hier - eine weitere drastische Vereinfachung - proportional zur Federkompression w, wenn sich Kugel und Oberfläche berühren:

Header

Wir lösen das Anfangswertproblem zu einer nichtlinearen Bewegungsgleichung. Die Nichtlinearität kommt aus der Kontaktbedingung der Kugel mit der Oberfläche.

Die Bewegungsgleichung ist stückweise linear, hier erfassen wir sie durch eine Kennlinie, die wir zwischen den linearen Bereich ausrunden. Das Ausrunden macht die numerische Integration schneller.

/*********************************************************/
/* MAXIMA script                                         */
/* version: wxMaxima 15.08.2                             */
/* author: Andreas Baumgart                              */
/* last updated: 2017-09-24                              */
/* ref: Kw23 (TM-C, Labor 5)                             */
/* description: finds the solution for                   */
/*              the nonlinear IVP                        */
/*********************************************************/

Declarations

Wir verwenden die Parameter

${\begin{array}{ll}k&=\displaystyle \kappa {\frac {m\;g}{r}}\\t_{B}&={\sqrt {\displaystyle {\frac {r}{g}}}}\\\kappa &=100\end{array}}$

Kontakt-Kennlinie

Die Kennlinie für den Kontakt definieren wir stückweise zu

$K(w,\epsilon )=\left\{{\begin{array}{ccl}0&{\text{für}}&w<-\epsilon \\\displaystyle k\cdot {\frac {\epsilon }{4}}\left(1+{\frac {w}{\epsilon }}\right)^{2}&{\text{für}}&-\epsilon <w<+\epsilon \\k\cdot w&{\text{für}}&\epsilon <w\\\end{array}}\right.$

Und so sieht sie dann aus:

/*********************************************************/
/* declarations */
assume(g>0, r>0);

/* parameter */
params: [k = kappa*m*g/r,  /*spring stiffness*/
         t[B] = sqrt(2*r/g), /*reference time*/
         kappa = 100];     /*dim'less spring stiffness*/

/**** define nonlinear spring characteristic ***/

K(w,epsilon) := if w <- epsilon then
                   0
                elseif w< epsilon then
                   (w/epsilon+1)^2*epsilon/4
                else
                w;         
plot2d(K(w,0.5),[w,-1,1], 
          [ylabel,"K/(k*W)->"], [xlabel,"w/W->"],
          [legend, "restoring force, ε=0.5"]);

Equilibrium Conditions

Freikörperbild

Die Gleichgewichtsbedingung lesen wir mit dem Prinzip von d'Alembert aus dem Freikörperbild ab:

$-m\cdot {\ddot {h}}-m\cdot g+K(w)=0$

Mit der dimensionslosen Zeit und Höhe

$\tau =t/t_{B}$,

${\tilde {h}}=h/r$

sowie

$K=k\cdot r\cdot {\tilde {K}}$

wird daraus

${\frac {\displaystyle m\cdot r}{\displaystyle {{t}_{B}^{2}}}}\left(\displaystyle {\frac {{d}^{2}}{d\,{{\tau }^{2}}}}\cdot {\tilde {h}}\right)-k\cdot r\cdot {\tilde {K}}\left(-{\tilde {h}},{\frac {1}{100}}\right)+g\cdot m=0$.

Weil wir ab hier mit den dimensionslosen Koordinaten weiterarbeiten, lassen wir die Tilde über dem h gleich wieder weg - so wird's übersichtlicher.

/********************************/
/* define ODE in dim'less coordinate h */
dgl : m*r/t[B]^2*'diff(h,t,2)-k*r*K(-h,1/100)+m*g = 0;
dgl : expand(solve(dgl, 'diff(h,t,2)))[1];
dgl : expand(subst(params,dgl));

Solving

Die Bewegungsgleichung zweiter Ordnung schrieben wir nun um als zwei Bewegungsgleichungen erster Ordnung zu

$\displaystyle {\frac {d}{d\,\tau }}\left({\begin{array}{c}h\\v\end{array}}\right)={\underline {f}}(h)$.

Die rechte Seite der gewöhnlichen Differentialgleichung ("Right-hand-side of ODE") lautet damit

${\begin{pmatrix}v\\100.0\cdot \left({\text{ if }}-1.0\cdot h<-0.01{\text{ then }}0.0{\text{ elseif }}-1.0\cdot h<0.01{\text{ then }}0.5\cdot {{h}^{2}}-0.01\cdot h+5.0\cdot {{10}^{-5}}\mathrm {else} -1.0\cdot h\right)-1.0\end{pmatrix}}$

Und diese können wir - hier mit dem Runge-Kutta-Verfahren 4.ter Ordnung - lösen.

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

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

Post-Processing

Wir tragen die Höhe h über der Zeit und h über v im Phasendiagramm auf:

Height h:

Lösung im Zeitbereich

Phase Diagram:

Lösung im Phasenraum

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

/* phase plot */
plot2d([discrete,H,V], [legend, false], [x,-1,10],
                   [ylabel,"V/1->"], [xlabel,"H/1->"]);

---

Links

• ...

Literature

• ...