Gelöste Aufgaben/Kw28

Aufgabenstellung

Das System besteht aus einer Kugel (Radius r, Masse m₂) und einer Plattform (Masse m₁). Wie skizziert ist die Plattform mit einer Parallelführung aus zwei Euler-Bernoulli-Balken (Biegesteifigkeit EI) elastisch gelagert. Aus der statischen Referenzkonfiguration wird die Kugel aus der Höhe H über der Plattform losgelassen. Kugel und Plattform stoßen also aufeinander und führen dann Ihre eigene vertikale Bewegung durch - bis zur nächsten Kollision. Der Stoß zwischen Kugel und Oberfläche sei ideal-elastisch.

Lageplan

Für das skizzierte System modellieren Sie die Kugel als elastisch, die elastisch gelagerte Plattform als starr. Gesucht ist eine numerische Lösung als Anfangswertproblem und die nichtlinearen Schwingungen der beiden Systemteile.

Lösung mit Maxima

Header

Die Kugel können Sie sich im unteren Teil durch eine Feder der Steifigkeit k₂ ersetzt denken - das geht analog zu Beispiel Kw23.

Die Federkraft K ist also Null, solange die Kugel die Oberfläche nicht berührt und sie ist proportional zur Federkompression w, wenn sich Kugel und Oberfläche berühren.

Kugel-Modell: elast. Kontakt mit Einfederung w.

Die beiden Körper haben jeweils einen Freiheitsgrad in vertikale Richtung. Die Bewegungsgleichungen für die beiden Körper sind stückweise linear (Kontakt: K = k w / kein Kontakt K = 0). Den Kontakt erfassen wir durch eine Kennlinie, die wir zwischen den beiden linearen Bereichen ausrunden. Das Ausrunden macht die numerische Integration schneller.

Wir lösen das Anfangswertproblem zu der zugeordneten nichtlinearen Bewegungsgleichung. Die Nichtlinearität kommt hier aus der Kontaktbedingung zwischen der Kugel und der Plattform.

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

Declarations

Die Koordinaten der Verschiebung der beiden Massen nennen wir

${\underline {Q}}(t)=\left({\begin{array}{l}u_{1}(t)\\u_{2}(t)\end{array}}\right){\text{ und damit }}\delta {\underline {Q}}=\left({\begin{array}{l}\delta u_{1}\\\delta u_{2}\end{array}}\right)$.

Für die beiden elastischen Balken verwenden die Ersatzfeder-Steifigkeit

$\displaystyle k_{1}=\displaystyle 2\cdot {\frac {12EI}{\ell ^{3}}}$.

Und als Abkürzungen verwenden wir

${\begin{array}{ll}\alpha &=\displaystyle {\frac {m_{2}}{m_{1}}}\\\kappa &=\displaystyle {\frac {k_{2}}{k_{1}}}\end{array}}$.

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

/* declare variational variables - see 6.3 Identifiers */
declare("δΠ", alphabetic);
declare("δW", alphabetic);
declare("δu", alphabetic);
declare("δQ", alphabetic);

/* coordinates */
 Q[t]: [ u[1](t), u[2](t)];
δQ[t]: [δu[1]   ,δu[2]   ];

/* abbreviations */
abbrev : [f(u[1](t)-u[2](t)) = f(w),
          gamma              = m[2]/m[1],
          k[2]               = kappa*k[1],
          m[2]               = alpha*m[1]];

Contact Characteristic

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

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

Kontakt-Kennlinie

Und so sieht die Kennlinie dann aus:

Die Parabel zwischen den beiden linearen Kennlinien-Stücken macht die Kraft K stetig differentierbar in w. Das macht die numerische Integration schneller - und genauer.

/**** define nonlinear spring characteristic ***/
C(w) := if w < -1 then
             0
        elseif w < 1 then
             1/4*(w+1)^2
        else
             w;        
plot2d(C(w),[w,-2,2],
          [ylabel,"C/1->"], [xlabel,"w/ε->"],
          [legend, "restoring force"]);

Equilibrium Conditions

Koordinaten

Die Gleichgewichtsbedingung konstruieren wir mit dem Prinzip der virtuellen Verrückungen - wir brauchen also kein Freikörperbild.

Dazu verwenden wir die Verschiebungs-Koordinaten u₁, u₂ der beiden Körper wie skizziert.

Für die allgemeine Gleichgewichtsbedingung

${\begin{array}{ll}\delta W&=0\\&=\delta W^{a}-\delta \Pi \end{array}}$

setzen wir für die virtuelle Arbeit von äußeren, eingeprägte Kräften (inklusive der D'Alembert'sche Trägheitskraft) an:

$\displaystyle \delta W^{a}=-\sum _{i=1}^{2}m_{i}\cdot \left({\ddot {u}}_{i}+g\right)\cdot \delta u_{i}$

sowie für die virtuelle Formänderungsenergie der beiden Ersatz-Federn:

$\delta \Pi =k_{1}u_{1}(t)\cdot \delta u_{1}+\underbrace {K(u_{1}(t)-u_{2}(t))} _{\displaystyle {\text{Kontakt-Kraft}}}\cdot (\delta u_{1}-\delta u_{2})$.

Sortieren nach den virtuellen Verrückungen und anschrieben der Gleichgewichtsbedingungen in Matrix-Form liefert

${\begin{pmatrix}{{m}_{1}}&0\\0&{{m}_{2}}\end{pmatrix}}\cdot {\begin{pmatrix}{{\ddot {u}}_{1}}\\{{\ddot {u}}_{2}}\end{pmatrix}}+{\begin{pmatrix}{{k}_{1}}&0\\0&0\end{pmatrix}}\cdot {\begin{pmatrix}{{u}_{1}}\\{{u}_{2}}\end{pmatrix}}=-g\cdot {\begin{pmatrix}{{m}_{1}}\\{{m}_{2}}\end{pmatrix}}+K\left(u_{1}-u_{2}\right)\cdot {\begin{pmatrix}-1\\+1\end{pmatrix}}$.

/********************************/
/* derive ODE in dim'less coordinates U using the PvV */

/* virtual work of system (equilibrium condition) */
PvV : [δW   = δW^a - δΠ,
       δW^a = - sum(m[i]* diff(u[i](t),t,2)*δu[i],i,1,2) - sum(m[i]*g*δu[i],i,1,2),
       δΠ   =   k[1]*u[1](t)*δu[1] + f(u[1](t)-u[2](t))*(δu[1]-δu[2])];
 
/* equlilbrium condition */
 
eom: subst(PvV[3],subst(PvV[2],subst(PvV[1],δW=0)));

/* linear parts first */
M : -funmake('matrix,makelist(makelist(coeff(coeff(lhs(eom),δQ[t][i]),diff(Q[t][j],t,2)),i,1,2),j,1,2));
K : -funmake('matrix,makelist(makelist(coeff(coeff(lhs(eom),δQ[t][i]),     Q[t][j]     ),i,1,2),j,1,2));
G :  funmake('matrix,makelist(coeff(
            subst(makelist(diff(Q[t][i],t,2)=0,i,1,2),
            subst(makelist(     Q[t][i]=0     ,i,1,2),[lhs(eom)])),δQ[t][j]), j,1,2));
/* what's left ?*/
F :  funmake('matrix,makelist(coeff(expand([lhs(eom)] +
                    (M.diff(Q[t],t,2)+K.Q[t]-G).transpose(δQ[t])),δQ[t][j]), j,1,2));

/* control print-out */
print(M,"∙",transpose(diff(Q[t],t,2)),"+",K,"∙",transpose(Q[t]),"=",G,"+",F)$

Dimensionless Equations of Motion

Mit einer passenden Bezugszeit und -länge machen wir die Bewegungsgleichungen dimensionslos.

Hier liefert die Periodendauer der Plattform - allein, ohne Kugel - die Bezugszeit T. Die statische Absenkung der Plattform - ohne Einfluss der Kugel - liefert die Bezugslänge L.

Aus der Bewegungsgleichung oben lesen wir für K(u₁-u₂) = 0 ab:

$T=\displaystyle {\frac {2\pi }{\omega _{0}}}{\text{ mit }}\omega _{0}^{2}=\displaystyle {\frac {k_{1}}{m_{1}}}$

und

$L=\displaystyle {\frac {m_{1}}{k_{1}}}g$.

Nun ist also

$\displaystyle \tau ={\frac {t}{T}}{\text{ mit }}(.)':={\frac {d}{d\tau }}(.)$

und

$\displaystyle U_{i}={\frac {u_{i}}{L}}{\text{ für }}i=1,2$

Die Kennlinie überführen wir entsprechend in

$K(w)=k_{2}\cdot L\cdot C\left(\displaystyle {\frac {W}{\epsilon }}\right){\text{ mit }}W(\tau )=U_{1}-U_{2}$.

Einsetzen liefert die zwei Differentialgleichungen

${\begin{array}{llll}\displaystyle {\frac {1}{4\cdot \pi ^{2}}}U''_{1}\left(\tau \right)&+{{U}_{1}}\left(\tau \right)&+1&=-\kappa \cdot C\left({\frac {W}{\epsilon }}\right)\\\displaystyle {\frac {\alpha }{4\cdot \pi ^{2}}}U''_{2}\left(\tau \right)&&+\alpha &=+\kappa \cdot C\left({\frac {W}{\epsilon }}\right)\end{array}}$.

/* dimensionless equations of motion */
dimless : ['diff(u[1](t),t,2) = L*'diff(U[1](tau),tau,2)/T^2,
           'diff(u[2](t),t,2) = L*'diff(U[2](tau),tau,2)/T^2,
                 u[1](t)      = L*U[1](tau),
                 u[2](t)      = L*U[2](tau),
                 L            = m[1]*g/k[1],
                 T^2          = (2*%pi/omega[1])^2,
                 omega[1]^2   = k[1]/m[1]];

eom : makelist((M.diff(Q[t],t,2)+K.Q[t]-(G+F))[i][1]=0,i,1,2);
eom : expand(subst(dimless,subst(abbrev,eom)));
eom : subst([f(w) = subst(dimless,subst(abbrev,L*k[2]*C(W/epsilon)))],eom);
eom : expand(solve(eom, ['diff(U[1](tau),tau,2),'diff(U[2](tau),tau,2)])[1]);

params: [kappa = 5, alpha = 1/2, W=(U[1](tau)-U[2](tau)), epsilon=1/100];
eom : expand(subst(params,eom));

Solving

Für die numerische Lösung schreiben wir die zwei Bewegungsgleichung zweiter Ordnung um als vier Bewegungsgleichungen erster Ordnung zu

$\displaystyle {\frac {d}{d\,\tau }}\left({\begin{array}{c}U_{1}\\U_{2}\\V_{1}\\V_{2}\end{array}}\right)={\underline {f}}(U_{1},U_{2}){\text{ mit }}V_{i}={\frac {dU_{i}}{d\,\tau }}$.

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

${\underline {f}}=\left({\begin{array}{c}1\\1\\\displaystyle 4\pi ^{2}(-U_{1}-1-\kappa \cdot C\cdot \left({\frac {W}{\epsilon }}\right))\\\displaystyle 4\pi ^{2}(-1+{\frac {\kappa \cdot C}{\alpha }}\cdot \left({\frac {W}{\epsilon }}\right))\end{array}}\right)$

Die dimensionslosen Parameter wählen wir zu

$\displaystyle \kappa =10,\alpha ={\frac {1}{2}},\epsilon ={\frac {1}{100}}$

und lösen die Bewegungsgleichungen numerisch - hier mit dem Runge-Kutta-Verfahren 4.ter Ordnung.

/********************************/
/* numerical solution of IVP */
times : subst([t0 = 0, tmax = 20, dt = 0.001],
                    [t, t0, tmax, dt]);
dgl1stOrder : subst(params,[v1,v2,rhs(eom[1]),rhs(eom[2])]);
dgl1stOrder : subst([U[1](tau) = h1,U[2](tau)=h2],dgl1stOrder);
stateVabs : [h1,h2,v1,v2];
initiVals : [-1,20,0,0];
 
ivs : rk(dgl1stOrder, stateVabs, initiVals, times)$

Post-Processing

Als Lösung tragen wir die Zeitverläufe der Höhen U₁, U₂ über der Zeit τ

Zeitverlauf der Lösung

und U₁, U₂  jeweils über V₁, V₂  im Phasendiagramm auf:

Phasendiagramm

/********************************/
/* solution samples */
Ti :  makelist(ivs[j][1],j,1,length(ivs))$
Hi : [makelist(ivs[j][2],j,1,length(ivs)),
     makelist(ivs[j][3],j,1,length(ivs))]$
Vi : [makelist(ivs[j][4],j,1,length(ivs)),
     makelist(ivs[j][5],j,1,length(ivs))]$

/* time-domaine plot */
plot2d([[discrete, Ti, Hi[1]],[discrete, Ti, Hi[2]]],
     [legend, "Plattform", "Kugel"], 
     [title, sconcat("start @: ",string(initiVals))],
                  [xlabel,"τ/1->"], [ylabel,"U/1->"]);
 
/* phase plot */
plot2d([[discrete,Hi[1],Vi[1]],[discrete,Hi[2],Vi[2]]],
       [legend, "Plattform", "Kugel"], 
                   [ylabel,"V/1->"], [xlabel,"H/1->"]);

✔ Wie genau ist die Lösung?:

Überlegen Sie, wie Sie bei diesem Ergebnis die Genauigkeit prüfen können. Das ist hier besonders einfach, weil es um ein konservatives System geht - also eins, bei dem keine Energie hinzugefügt oder dissipiert wird.

---

Links

• ...

Literature

• ...