Gelöste Aufgaben/Kw26
Aufgabenstellung
Ein Aggregat der Masse 2m ist durch zwei Gummifedern wie skizziert gelagert und wird senkrecht durch seine Gewichtskraft und eine äußere Last F=2 H sin(Ω t) belastet.
Für dieses System mit einem nichtlineares Federelement soll die Lösung als Anfangswertproblem gefunden und mit der linearen Näherungslösung verglichen werden.
Gesucht sind
- die Lösung der zugeordneten linearisierten Bewegungsgleichung,
- die numerische Lösung des nichtlinearen Problems.
Für die Gummifeder gibt der Lieferant folgende eine Kennlinie an:
Gegeben:
- u0,
- m, g, H, Ω'
Lösung mit Maxima
Header
Gummi-Elemente zeigen eine deutliche Nichtlinearität bei großen Auslenkungen. Wie man diese Nichtlinearität bei dynamischer Belastung berücksichtigt, zeige ich Ihnen hier.
/*********************************************************/
/* MAXIMA script */
/* version: wxMaxima 15.08.2 */
/* author: Andreas Baumgart */
/* last updated: 2018-11-27 */
/* ref: Kw26 (TM-C, Labor 6) */
/* description: finds the solution for */
/* the nonlinear IVP */
/*********************************************************/
Declarations
Wir berücksichtigen nur die vertikale Bewegung des Aggregats, dafür nehmen wir die Koordinate u:
Als System-Parameter setzten wir an:
(Kraft-Amplitude als Vielfaches der Gewichtskraft) und
(Erreger-Periodendauer ist ein Vielfaches der Schwingungsdauer des zugeordneten linearen Systems).
/*********************************************************/
/* declarations */
assume(g>0, u[0]>0);
/* parameter */
params: [F(t) = 2*gamma*m*g*sin(Omega*t),
Omega=2*%pi/T[Omega],
b = 2*D*m*omega[0],
omega[0] = 2*%pi/T[l],
T[Omega]=eta*T[l]];
Spring Characteristic
Die Gummi-Element Federkennlinie approximieren wir über die punktsymmetrische Funktion
Die beiden Konstanten c1 und c2 erhalten wir aus den Bedingungen
- ,
die wir aus der Kennlinie ablesen zu
So sieht die Kennlinie dann aus:
/*********************************************************/
/* define friction characteristics (Reib-Kennlinie) */
curve : c[1]*u[t]/u[0] + c[3]*(u[t]/u[0])^3;
sol: solve([subst([u[t]=u[0]],curve) = 1,
subst([u[t]=2*u[0]],curve) = 3],[c[1],c[3]])[1];
curve : subst(sol,curve);
plot2d(subst([u[t]=mu*u[0]],curve), [mu,-1,2],
[legend,"characteristic"],
[style, [lines,3,1]], [xlabel, "ut/u0 →"],
[ylabel, "K/mg →"])$
Equilibrium Conditions
Die Bewegungsgleichungen des Systems bekommen wir direkt aus dem Kräftegleichgewicht am Körper (mit dem Prinzip von d'Alembert) zu
Damit die "homogene Lösung" der Bewegung (der Einschwingvorgang in die periodische Lösung) nach endlicher Zeit weg ist, mogeln wir noch etwas Dämpfung hinzu:
- .
Wir brauchen auch die linearisierte Bewegungsgleichung. Dazu entwickeln wir die Federkennlinie um die stationäre Ruhelage u = u0 des Aggregats in eine Taylor-Reihe bis zum linearen Glied, setzen also
und finden die linearisierte Form der Bewegungsgleichung zu
- .
/*********************************************************/
/* define ODE in dim'less coordinate U */
eom : 2*m*'diff(u[t],t,2) + 2*b*'diff(u[t],t,1) + 2*m*g*curve = 2*m*g+F(t);
/* static solution (F(t)=0, diff(u[t],t,2)=0) */
stationary: solve(subst([F(t)=0,
'diff(u[t],t,1)=0,
'diff(u[t],t,2)=0,
u[t]=u[S]], eom),u[S])[3];
linear : subst([u[t]=u[0]],
subst(['diff(u[t],t,2)='diff(u[l],t,2),
'diff(u[t],t,1)='diff(u[l],t,1)], lhs(eom))) +
subst([u[t]=u[0]],diff(lhs(eom), u[t]))*u[l] = rhs(eom);
linear : linear - 2*m*g;
nonlin : expand(solve(eom, 'diff(u[t],t,2)));
Dimensionless Representation
Dimensionslos machen wir die Bewegungsgleichung mit der Referenz-Auslenkung u0 und der Periodendauer T0 der zugeordneten linearen Bewegungsgleichung, also
- und
Umschreiben der linearen und nichtlinearen Bewegungsgleichungen in dimensionslose Form liefert mit
- die nichtlineare Form
und - die linearisierte Form
.
/********************************/
/* redefine ODE in dim'less coordinate U */
dimless: ['diff(u[l],t,2)=u[0]*'diff(mu[l],tau,2)/T[l]^2,
'diff(u[t],t,2)=u[0]*'diff(mu[n],tau,2)/T[l]^2,
'diff(u[l],t,1)=u[0]*'diff(mu[l],tau,1)/T[l],
'diff(u[t],t,1)=u[0]*'diff(mu[n],tau,1)/T[l],
u[l]=mu[l]*u[0],
u[t]=mu[n]*u[0],
t=tau*T[l]];
linear : subst(dimless,subst(params,linear));
linear : linear/coeff(lhs(linear),'diff(mu[l],tau,2));
linear : expand(linear);
dimless : append(dimless,
[solve(coeff(lhs(linear),mu[l]) = 4*%pi^2,T[l])[2] ]);
linear: subst(dimless,linear);
nonlin : expand(subst(dimless,subst(params,3*%pi^2*nonlin/g)));
Solving
Teil A: Partikulare Lösung der linearen Bewegungsgleichung
Von der linearen Bewegungsgleichung wollen wir nur die partikulare Lösung - also die Lösung nach Abklingen des Einschwingvorgangs (dafür brauchten wir die Dämpfung!). Wir wählen als Ansatz
(Ansatz vom Typ der rechten Seite)
und erhalten für die Koeffizienten
- .
Die Lösung plotten wir unten.
Teil B: Numerische Lösung des Anfangswertproblems
Für die numerische Lösung brauchen die Bewegungsgleichung als Differentialgleichung erster Ordnung als
- .
Die lösen wir nun als Anfangswertproblem
mit
- .
Als Anfangsbedingungen wählen wir die stationäre Ruhelage des Aggregats, also
- .
Und wir führen die numerische Integration mit drei Parameter-Sätzen durch:
- Überkritische Anregung:
Hier ist η = 1/2, also erfolgt die Anregung mit einer Periodendauer, die halb so groß ist, wie die Periodendauer des linearisierten Systems - oder: die Erregerfrequenz ist doppelt so hoch wie die Eigenfreuqenz. - Anregung in der Nähe der Resonanzfrequenz:
Mit η = 1.1 regen wir in der Nähe der Erregerfrequenz des linearen Systems an - etwas darunter, die Phase des Systems ist also fast noch konphas mit der Erregung. - Unterkritische Anregung:
Hier ist Erregerfrequenz nur hald so groß, wie die Eigenfrequenz des linearen Systems, die Phase des Systems ist also konphas zur Erregung.
/***********************************************/
/* inhomogenous part of linearized solution */
ansatz: [mu[l] = U[c]*cos((2*%pi*tau)/eta)+U[s]*sin((2*%pi*tau)/eta)];
linSol: subst(ansatz,linear);
linSol: expand(ev(linSol,nouns));
linSol: solve([coeff(linSol,cos((2*%pi*tau)/eta)),
coeff(linSol,sin((2*%pi*tau)/eta))],[U[c],U[s]])[1];
numpar: [[eta=0.5,gamma=1, D=0.1],
[eta=1.1,gamma=1, D=0.1],
[eta=2.0,gamma=1, D=0.1]];
/* plot results, time domain */
for i: 1 thru length(numpar) do
plot2d(subst(numpar[i],subst(linSol, subst(ansatz,mu[l]))),[tau,0,2],
[xlabel, "τ ->"], [ylabel, "μ ->"], [legend, sconcat("@: ",string(numpar[i]))]);
/********************************/
/* numerical solution of IVP */
times : subst([t0 = 0, tmax = 20, dt = 0.001],
[t, t0, tmax, dt]);
dgl1stOrder : [V,
subst(solve(nonlin, 'diff(mu[n],tau,2)),'diff(mu[n],tau,2))];
dgl1stOrder : subst(['diff(mu[n],tau,1)=V,mu[n]=U, tau=t],dgl1stOrder);
stateVabs : [U,V];
initiVals : [1,0];
/* solution of IVP (ivs) */
for i: 1 thru length(numpar) do
ivs[i] : rk(subst(numpar[i], dgl1stOrder), stateVabs, initiVals, times)$
Post-Processing
Für die drei Anfangsbedingungen erhalten wir die Auslenkungen U(τ):
# | Auslenkung | Geschwindigkeit |
---|---|---|
1) | ||
2) | ||
3) |
Für das nichtlineare System können wir uns den Einschwingvorgang im Phasendiagramm ansehen:
Und stellen fest:
- nur für die überkritische Lösung ist die Auslenkungen des Systems so klein, dass eine Linearisierung gerechtfertigt ist.
Für die Parameter-Sätze 3 und besonders 2 entstehen Schwingungen mit hohen Amplituden, die wir in jedem Fall als nichtlineares Problem lösen müssen.
Und im eingeschwungenen Zustand können wir die Lösung für lineares und nichtlineares System einfach vergleichen. Die Lösungen des linearen Systems sind immer Ellipsen, die Lösungen des nichtlinearen Systems sind deutlich "kunstvoller":
/********************************/
/* plot results, time domain */
for i: 1 thru length(numpar) 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]],
[legend, sconcat("@: ",string(numpar[i]))],
[xlabel,"τ/1->"], [ylabel,"μ/1->"]));
/* phase plot */
curves : makelist([discrete,uD[i],vD[i]],i,1,length(numpar))$
plot2d(curves, [legend, false], [x,-2.2,3.2],
[ylabel,"V/1->"], [xlabel,"U/1->"]);
/********************************/
/* phase plot - linear and nonlinear */
curves : makelist(subst(numpar[i],subst([tau=t],subst(linSol,subst(ansatz,[parametric, 1+mu[l], 'diff(mu[l],tau),[t,0,2]])))),i,1,length(numpar));
curves : ev(curves,nouns);
/* get last 2 "seconds" of numerical solution */
for i: 1 thru length(numpar) do
(uDL[i] : makelist(uD[i][j],j,length(uD[i])-round(2/times[4]),length(uD[i])),
vDL[i] : makelist(vD[i][j],j,length(vD[i])-round(2/times[4]),length(vD[i])));
curves : append(curves,makelist([discrete,uDL[i],vDL[i]],i,1,length(numpar)))$
plot2d(curves, [legend, false], [x,-2.2,3.2],
[style, [lines,2,1], [lines,2,2], [lines,2,3], [lines,3,1], [lines,3,2], [lines,3,3]],
[ylabel,"V/1->"], [xlabel,"U/1->"]);
🧨 Unzulässige Belastung von Gummifedern: |
Das Phasendiagramm zeigt: die Gummifeder wird für die Parameter-Sets 2 und 3 deutlich auf Zug belastet - das ist nicht zulässig! |
Links
- ...
Literature
- ...