Gelöste Aufgaben/Kw25

Aufgabenstellung

Ein Körper (Masse m) liegt auf einem Band, das über eine Rolle (Radius r) mit konstanter Drehgeschwindigkeit Ω angetrieben wird. Der Körper ist über eine Feder (Steifigkeit k) mit der Umgebung verbunden. Zwischen Körper und Band wirkt trockene Reibung (μ, μ₀).

Lageplan

Gesucht ist die numerische Lösung für den Klotz auf dem Band als Anfangswertproblem, untersucht werden selbsterregte Schwingungen.

Lösung mit Maxima

Header

"Stick-Slip" Schwingungen sind klassische selbsterregte Schwingungen. Selbsterrung heißt: das System lenkt den Energiefluss im System so, dass Schwingungen "aus sich heraus" angeregt werden. Charakteristisch für "Stick-Slip" Schwingungen ist, dass zwei Körper zeitweise aneinander haften, dann auseinander gerissen werden und dann aneinander reiben - bis sie wieder aneinander haften.

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

Declarations

Im Zentrum dieser Aufgabe steht die Kennlinie (engl.: Characteristic), die den Zusammenhang zwischen Tangentialkräft (Haftkraft, Reibkraft), Normalkraft und Relativgeschwindigkeit erfasst.

Dafür brauchen wir hier drei Parameter.

• a=μ / μ_0;
• b : erfasst, welchen Ausschnitt der Sinus-Funktion wir nutzen;
• v₀: eine (kleine) Kriechgeschwindigkeit.

/*********************************************************/
/* declarations */
assume(a>0,b>0);

/* parameter */
params: [a=1/2, b = 5/6, v = 1., V0=0.001];

Friction Characteristic

Trockene Reibung im Modell beschreibt man mit zwei System-Parametern:

• dem Reibungskoeffizient μ und
• dem Haftungskoeffizient μ₀.

Haften Körper und Unterlage aneinander (Relativgeschwindigkeit v_r ist "Null"), so ist die Haftkraft

$|H|\leq \mu _{0}\cdot N$

mit der Normalkraft N. Bewegen sich die beiden relativ zueinander, so schreibt man

$|R|=\mu \cdot N$.

Die Beträge bedeuten hier: die Kräfte sind gegen die Bewegung orientiert.

Für die numerische Lösung des Anfangswertproblems ist diese Bescheibung eine Katastrophe!

1. Wir haben eine Kennlinie mit Sprung zwischen μ und μ₀.
2. Das Vorzeichen des Funktionswerts für Haftung hängt davon ab, mit welcher Orientierung eine resultierende Kraft auf den Körper wirkt.

Der Sprung ist eine effektive Handbremse für jeden Algorithmus zur numerischen Integration der Bewegungsgleichung, der von einer stetigen und stetig differenzierbaren "Rechten Seite" der Bewegungsgleichung ausgehen (vgl. Runge-Kutta-Verfahren 4.ter Ordnung) - also für praktisch alle Algorithmen. Und dass die Haftkraft nicht eindeutig von den Zustandsgrößen des Systems abhängt, hebelt jedes Standard-Lösungsverfahren aus.

Kommerzielle Programme führen deshalb "Schaltstellen" ein, in denen der Lösungsalgorithmus an diesen Unstetigkeiten dann "etwas anderes" macht. Oft wird die Lösung des Anfangswertproblems dadurch aber sehr ineffizient. Hier schreiben wir deshalb eine Kennlinie an, bei der die Tangentialkraft (hier nennen wir sie Reibkraft) stetig von der Relativgeschwindigkeit v_r zwischen Körper und Band abhängt.

Reibkennlinie

Wir suchen eine Kennlinie wie hier im Bild blau und rot gekennzeichnet. Die klassische Formulierung ist in grau hinterlegt. Die Kennlinie unserer Formulierung nennen wir Reibkennlinie

$R(\mu ,\mu _{0},v_{r})$.

Das Haften ist hier nur ein Sonderfall, den die Kennlinine mit abbildet.

Die gesuchte Kennlinie

$R=\mu (v_{r})\cdot N$

muss punkt-symmetrisch sein (die Reibkraft ändert ihr Vorzeichen mit der Orientierung der Relativgeschwindigkeit). Welche Funktionen wir verwenden ist egal, hier probieren wir Sinus und Exponential-Funktionen aus. Wie diese Funktioen aneinader-gestückelt werden, steht im Lexikon unter "Reibkennlinie".

Wir arbeiten mit

$\displaystyle \mu =\left\{{\begin{array}{lllll}\left(a-\sin \left(\pi \cdot {{b}^{2}}\right)\right)\cdot {{e}^{\displaystyle -{\frac {\pi \cdot b\cdot \mathrm {cos} \left(\pi \cdot {{b}^{2}}\right)\cdot \left(b+\upsilon \right)}{\mathrm {sin} \left(\pi \cdot {{b}^{2}}\right)-a}}}}-a&{\text{für}}&\upsilon &<-b\\\sin(\pi \;b\;\upsilon )&{\text{für}}-b<&\upsilon &<+b\\a-\left(a-\sin \left(\pi \cdot {{b}^{2}}\right)\right)\cdot {{e}^{\displaystyle -{\frac {\pi \cdot b\cdot \cos \left(\pi \cdot {{b}^{2}}\right)\cdot \left(b-\upsilon \right)}{\sin \left(\pi \cdot {{b}^{2}}\right)-a}}}}&{\text{für}}+b<&\upsilon &\end{array}}\right.$.

/*********************************************************/
/*  define friction characteristics (Reib-Kennlinie) */
ch : [sin(b*%pi*upsilon),0,0];

ch[2]: E*%e^(kappa*(upsilon+b))-a;
ch[2] : subst(solve([subst([upsilon=-b],ch[2]=ch[1]),subst([upsilon=-b],diff(ch[2]=ch[1],upsilon))],[E,\kappa])[1],ch[2]);
ch[3] : subst([upsilon=-upsilon],-ch[2]);
plot2d([[parametric, t, subst(params,subst(t,upsilon,ch[1])),subst(params,[t,-b, b])],
        [parametric, t, subst(params,subst(t,upsilon,ch[2])),subst(params,[t,-3,-b])],
        [parametric, t, subst(params,subst(t,upsilon,ch[3])),subst(params,[t,+b,+3])]],[legend,"sec. 0","sec -1","sec +1"],
        [y,-1,1],
        [style, [lines,3,1], [lines,3,2], [lines,3,2]], [xlabel, "Geschwindigkeit v/V0 →"], [ylabel, "Reib-Koeffizient μ/1 →"])$

/**** define nonlinear friction characteristic ***/
fric(upsilon, a,b) := 
  if upsilon < -b then
     (a-sin(%pi*b^2))*%e^(-(%pi*b*cos(%pi*b^2)*(b+upsilon))/(sin(%pi*b^2)-a))-a       /*ch[2]*/
  elseif upsilon > +b then
     a-(a-sin(%pi*b^2))*%e^(-(%pi*b*cos(%pi*b^2)*(b-upsilon))/(sin(%pi*b^2)-a))       /*ch[3]*/
  else
     sin(%pi*b*upsilon)                                                               /*ch[1]*/
     ;

Equilibrium Conditions

Freikörperbild

Die Bewegungsgleichungen des Systems bekommen wir direkt aus dem Kräftegleichgewicht am Körper (mit dem Prinzip von d'Alembert) zu

$m\;{\ddot {u}}+k\;u-R(\mu ,\mu _{0},v_{r})=0$.

Dimensionslos machen wir die Bewegungsgleichung mit

• der Referenzzeit T aus
  ${\begin{array}{ll}\omega _{0}&=2\pi /T\\&={\sqrt {k/m}}\end{array}}$
  mit der Eigenkreisfrequenz
  $\omega _{0}$
  des frei schwingenden Körpers (R=0) und
• der Auslenkung des Körpers an der Feder im Erd-Schwerefeld
  $\displaystyle {\hat {u}}={\frac {m\;g}{k}}$

Wir erhalten mit

$u={\hat {u}}\cdot U{\text{ und }}t=T\cdot \tau$

als Bewegungsgleichung

$\displaystyle {\frac {{{\mu }_{0}}\cdot g\cdot m}{4\cdot {{\pi }^{2}}}}\cdot {\frac {{d}^{2}U}{d\,{{\tau }^{2}}}}+{{\mu }_{0}}\cdot g\cdot m\cdot U-g\cdot m\cdot \mu (v_{r})=0$

bzw.

$\displaystyle {\frac {{d}^{2}U}{d\,{{\tau }^{2}}}}+4\pi ^{2}\cdot U-4\pi ^{2}\cdot {\frac {\mu (v_{r})}{\mu _{0}}}=0$

/*********************************************************/
/* define ODE in dim'less coordinate U */
equil : m*'diff(u,t,2) + k*u - R = 0;

fstOrd : ['diff(U,tau,2)='diff(V,tau,1)];
ref: [solve((2*%pi/T)^2 = k/m,T)[2]];

equil : subst([R=m*g*mu[0]*fric((v-V)/V0,a,b), 'diff((mu[0]*m*g/k*U),tau,2)=mu[0]*m*g/k*'diff(U,tau,2)],subst(u=mu[0]*m*g/k*U, subst(ref,subst(['diff(u,t,2)= 'diff(u,tau,2)/T^2],equil))));
equil : [subst(fstOrd,equil), 'diff(U,tau,1)=V];

dydt : solve(equil,['diff(V,tau,1),'diff(U,tau,1)])[1];

Solving

Die Differentialgleichung erster Ordnung

${\begin{array}{ll}\displaystyle {\frac {d\,U}{d\,\tau }}&=V,\\\displaystyle {\frac {d\,V}{d\,\tau }}&=4\cdot {{\pi }^{2}}\cdot \left(\mu \left(v_{r}\right)-U\right)\end{array}}$

lösen wir nun als Anfangswertproblem

${\underline {\dot {q}}}={\underline {f}}({\underline {q}})$

mit

${\underline {q}}=\left({\begin{array}{c}U\\V\end{array}}\right)$.

Als Anfangsbedingungen wählen wir diese vier:

${\underline {q}}_{0,1}=\left({\begin{array}{l}1\\0\end{array}}\right)$,

${\underline {q}}_{0,2}=\left({\begin{array}{l}-1\\+0\end{array}}\right)$,

${\underline {q}}_{0,3}=\left({\begin{array}{l}0.\\0.001\end{array}}\right)$,

${\underline {q}}_{0,4}=\left({\begin{array}{l}0\\3\end{array}}\right)$.

/********************************/
/* numerical solution of IVP */
		
times : subst([t0 = 0, tmax = 4, dt = 0.001],
                    [t, t0, tmax, dt]);
dgl1stOrder : [rhs(dydt[2]),float(subst(params,rhs(dydt[1])))];
stateVabs : [U,V];
initiVals : [[1,0],[-1,0],
             [0,0.001],[0,3]];
/* solution of IVP (ivs) */
for i: 1 thru length(initiVals) do
   ivs[i] : rk(dgl1stOrder, stateVabs, initiVals[i], times)$

Post-Processing

Für die vier Anfangsbedingungen finden wir diese Zeitverläufe der Auslenkung U(τ):

a)	b)
c)	d)

Unabhängig von den Anfangsbedingungen läuft das System sofort in die "Stick-Slip" Schwingung hinein.

Phasen Diagramme:

a)

b)

Hier erkennt man gut, wie der Körper am Band haften bleibt - für V=1 sind Band und Körper-Geschwindigkeit gleich.

/********************************/
/* plot results, time domain  */
for i: 1 thru length(initiVals) do
	(T[i] : makelist(ivs[i][j][1],j,1,length(ivs[i])),
	 uD[i] : makelist(ivs[i][j][2],j,1,length(ivs[i])),
	 vD[i] : makelist(ivs[i][j][3],j,1,length(ivs[i])),
	 plot2d([discrete, T[i], uD[i]],
	      [title, sconcat("start @: ",string(initiVals[i]))],
	                   [xlabel,"τ/1->"], [ylabel,"U/1->"]));

/* phase plot */
curves : makelist([discrete,uD[i],vD[i]],i,1,length(initiVals))$
plot2d(curves, [legend, false], [x,-1.2,1.2],
                   [ylabel,"V/1->"], [xlabel,"U/1->"]);
plot2d(curves, [legend, false], [x,-0.1,0.8], [y,-1,2],
                   [ylabel,"V/1->"], [xlabel,"U/1->"]);

---

Links

• ...

Literature

• ...