Gelöste Aufgaben/SKEB

Aus numpedia
Zur Navigation springen Zur Suche springen


Aufgabenstellung

Die Bewegung des Balkens wird durch das Zusammenspiel von elastischen Verformungen und Trägheitskräften bestimmt. Man nennt das "Schwingungen von Kontinua" - diese untersuchen wir hier. Der zentrale Aufgabe besteht in der Berechnung der homogenen Lösung - und der Anpassung der Lösungsanteile an die Anfangsbedingungen.


Lageplan

Gesucht ist analytische Lösung für Schwingungen des Euler-Bernoulli-Balkens beim Loslassen aus der enspannten Ruhelage.


Lösung mit Maxima

Header

Für die mathematische Behandlung - insbesondere der Auflösung quadratischer Gleichungen - setzen wir in Maxima voraus, dass

.

/*******************************************************/
/* MAXIMA script                                       */
/* version: wxMaxima 16.04.2                           */
/* author: Andreas Baumgart                            */
/* last updated: 2018-01-15                            */
/* ref: oscillations of continuous media               */
/* description: analytic solutioun for EBB             */
/*******************************************************/

/* settings */
workdir: "C:\\Users\\X...X\\plots\\";

/* assumptions */
assume(rho>0, E>0,I>0,A>0,l[0]>0,alpha>0,gamma>0);

/* abbreviations used later .... */
abbrev : [alpha = (sqrt(omega[0])*A^(1/4)*rho^(1/4)*l[0])/(E^(1/4)*I^(1/4)),
          gamma = ((cosh(alpha)+cos(alpha)))/(sinh(alpha)+sin(alpha))];
abbrev : append(abbrev, solve(abbrev,[sqrt(omega[0]),cosh(alpha)])[1]);
abbrev : append(abbrev, [abbrev[3]^2,abbrev[3]^3]);




Equations of Motion

In der Gleichgewichtsbeziehung für den Euler-Bernoulli-Balken setzen wir als eingeprägte, äußere Streckenlast q(x,t) die D'Alembert'sche Trägheitskraft und die Gewichtskraft an, also

Die Lösung der linearen, partielle Bewegungsgleichung

setzt sich dann aus zwei Lösungsanteilen, der partikularen Lösung wp und der homogenen Lösung wh, zusammen; wir schreiben

.

/* equation of motion               */
eom: rho*A*diff(w(x,t),t,2) + E*I*diff(w(x,t),x,4) = rho*A*g;




Particular Solution

Die partikulare Lösung wp erfüllt die "rechte Seite" der Bewegungsgleichung, also ϱ A⋅g:

.

Die rechte Seite ist zeit-unveränderlich - so auch die partikulare Lösung.

Wir integrieren die Bewegungsgleichung vier Mal und erhalten

.

Die vier Integrationskonstanten Ci müssen wir nun an die Randbedingungen

anpassen, wir erhalten mit dem linearen Gleichungssystem

die partikulare Lösung

.
Statische Auslenkung
Die maximale Auslenkung - am rechten Rand - nutzen wir als Bezugslänge

Und so sieht wp aus:


/***************************************************************************/
/* p a r t i c u l a r   s o l u t i o n                                   */
/***************************************************************************/
peom : subst(['diff(w(x,t),t,2)=0],eom);
w[p] : integrate(integrate(integrate(integrate(peom,x),x),x),x);
w[p] : subst([%c1=C[3],%c2=C[2],%c3=C[1],%c4=C[0]],w[p]);

BouCon : [subst([x= 0  ],subst([      w(x,t)     =0],     w[p]     )),
          subst([x= 0  ],subst(['diff(w(x,t),x,1)=0],diff(w[p],x,1))),
          subst([x=l[0]],subst(['diff(w(x,t),x,2)=0],diff(w[p],x,2))),
          subst([x=l[0]],subst(['diff(w(x,t),x,3)=0],diff(w[p],x,3)))];
w[p] : subst(solve(BouCon,makelist(C[i],i,0,3))[1],w[p]);
w[p] : ratsimp(lhs(w[p]/(E*I)) = subst([x=l[0]*xi], rhs(w[p])/(E*I)));

params: [w[s] = subst([xi=1],subst(w[p],w(x,t))), omega[0]=2*%pi/T];

plot2d(ratsimp(subst(params,subst(w[p],w(x,t))/w[s])),[xi,0,1], 
               [legend, "w[p]/w[s]"],
               [style, [lines,3,1]], [xlabel, "ξ →"], [ylabel, "← w"],
               [gnuplot_preamble, "set yrange [1.1:-0.1] reverse"])$




Homogeneous Solution

Die homogene Bewegungsgleichung

lässt Lösungen vom Typ

zu. Wir bekommen zu jedem ω0 vier κ

,

also


Die komplexen Exponenten deuten schon darauf hin: im Allgemeinen sind die Koeffizienten Wi ebenfalls komplex. Für dieses Beispiel macht es deshalb viel Sinn, diese allgemeine Lösung der homogenen Differentialgleichung etwas anders hinzuschreiben. Wir nutzen dafür:

  • die Hyperbelfunktionen
  • die Euler'sche Formel

Dann lautet die allgemeine Lösung  

,

bei der wir davon ausgehen dürfen, dass alle Wi reellwertig sind.

Wie bei der partikularen Lösung wp sind es die vier Integrationskonstanten Wi, die wir nun an die Randbedingungen

anpassen müssen. Wir erhalten das lineare Gleichungssystem

,

vereinfacht und in Matrix-Schreibweise

.

Für dieses lineare, homogenen Gleichungen gibt es nur dann nicht-triviale Lösungen, wenn die Matrix einen Rangabfall >= 1 hat, also

.

Wir finden die charakteristische Gleichung

.

Beim Plotten der Funktion sehen wir bereits, dass es mehr als eine Nullstelle gibt, hier bei α≈+/- 2 und α≈+/- 5.

Charakteristische Gleichung

Die charakteristische Gleichung können wir umformen zu

und beide Anteile auftragen. So wird klar, dass die beiden Funktionen sich wegen der Periodizität der cos-Funktion unendlich oft schneiden - es also unendlich viele Nullstellen gibt:

Nullstellensuche (Ausschnitt)
Nullstellensuche

Analytisch finden wir keine Lösungen dieser charakteristischen Gleichung - wir müssen sie numerisch bestimmen – hier mit der Maxima-Funktion "find_root". Die ersten (positiven) Lösungen sind

.

Sie zu finden ist einfach, weil die αi für i>1 ungefähr π auseinander liegen - da können wir keine Nullstelle verpassen.

In diesen αi stecken über  auch der Periodendauern, z.B. für die erste Mode

.

Wie bei einem normalen Eigenwertproblem berechnen wir nun die Eigenvektoren W aus der Beziehung

.

Aus unendlich vielen Lösungen αi folgen nun also auch unendlich viele Wj, die wir zu

bestimmen. Die homogene Lösung der Bewegungsgleichung ist also

.

Denken Sie daran: die Eigenkreisfrequenz ω0j ist dabei wieder eine Funktion von αi - siehe oben. Die Ci sind wieder Konstanten, die wir brauchen, um die Lösung an die Anfangsbedingungen anzupassen. Zunächst schauen wir uns die Lösung jeweils zu einem αi an. Diese Funktionen heißen Modalformen ϕ und deren Schwingungen können wir plotten:

ModeModalform ϕj ModeModalform ϕj
#1
#2
#3
#4
#5
#6
Modalformen

Diese Moden ϕj können wir auch einzeln übereinander darstellen:


/***************************************************************************/
/* h o m o g e n e o u s   s o l u t i o n                                 */
/***************************************************************************/
ansatz: w(x,t) = W[i] * %e^(kappa*x) * %e^(%i*omega[0]*t);
heom : subst(ansatz,subst([g=0],eom));
heom : ratsimp(subst(ansatz,ev(heom,nouns)/w(x,t)));
sol[1] : subst(abbrev,solve(heom,kappa));
/* ansatz: w(x,t) = sum(subst(sol[1][i],subst(ansatz,w(x,t))),i,1,length(sol[1])); */
/* this is possible - but results in complex W[i] - too demanding */
ansatz : w(x,t) = subst(sol[1][1], W[1]*%e^(kappa*x)*%e^(%i*omega[0]*t))
                 +subst(sol[1][2], W[2]*%e^(kappa*x)*%e^(%i*omega[0]*t))
                 +subst(sol[1][3], W[3]*%e^(kappa*x)*%e^(%i*omega[0]*t))
                 +subst(sol[1][4], W[4]*%e^(kappa*x)*%e^(%i*omega[0]*t));
/* .. but we can rephrase w(x,t) with
             https://de.wikipedia.org/wiki/Sinus_hyperbolicus_und_Kosinus_hyperbolicus

      to have only real-values W[i]:            */
ansatz : w(x,t) =( W[1]*sinh(kappa*x)
                  +W[2]*cosh(kappa*x)
                  +W[3]*sin (kappa*x)
                  +W[4]*cos (kappa*x))*%e^(%i*omega[0]*t);

/* test:subst(ansatz,heom);ratsimp(ev(%,nouns));subst([x=xi*l[0]],%);subst(abbrev[1],%);*/

/* boundary conditions */
BouCon : [subst([x=  0 ],      subst(ansatz,w(x,t))     )/%e^(%i*omega[0]*t)=0,
          subst([x=  0 ], diff(subst(ansatz,w(x,t)),x  ))/%e^(%i*omega[0]*t)=0,
          subst([x=l[0]], diff(subst(ansatz,w(x,t)),x,2))/%e^(%i*omega[0]*t)=0,
          subst([x=l[0]], diff(subst(ansatz,w(x,t)),x,3))/%e^(%i*omega[0]*t)=0];
D : submatrix(augcoefmatrix(expand(BouCon),makelist(W[i],i,1,4)),5);
/* what happens to the six alpha=0 ???? */
D : diag_matrix (1, 1/kappa, 1/kappa^2,1/kappa^3).subst(abbrev,D);
D : subst(solve(sol[1][4],kappa),D);
CharEqua : expand(determinant(D))=0;
CharEqua : ratsimp(subst([sin(alpha)^2=1-cos(alpha)^2,sinh(alpha)^2=cosh(alpha)^2-1], CharEqua/2));

/* numerically find zeros for omega[0] (alpha) */
plot2d(lhs(CharEqua),[alpha,-6,6], [y,-6,6], [ylabel,"lhs(CharEqua)→"], [xlabel,"α→"], [style, [lines,3]])$
plot2d([cos(alpha),-1/cosh(alpha)],[alpha,0,60], [ylabel,"f(α)→"], [xlabel,"α→"], [style, [lines,3]])$
guessed: [2,5,8,11,14,17];
zeros : makelist(alpha = find_root(CharEqua/alpha^6, alpha, guessed[i]-1/2, guessed[i]+1/2),i,1,length(guessed));
sol[2] : makelist(omega[0] = subst(zeros[i],subst(abbrev[3],sqrt(omega[0]))^2),i,1,length(zeros));
params : append(params, [T[ref] = subst(zeros[1],subst(abbrev,2*%pi/omega[0]))]);

/* use upper triangular matrix instead of D */
[P,L,U] : get_lu_factors(lu_factor(D))$

/* get eigenvectors */
EigVec : subst(solve(makelist(sum(U[j][i]*W[i],i,1,4)=0,j,1,3),[W[1],W[2],W[3]])[1],
                                   matrix([W[1]],[W[2]],[W[3]],[W[4]]));
EigVec : subst(abbrev,subst([W[4]=1],EigVec));
Gammas : makelist(gamma = subst(zeros[i],subst(abbrev[2],gamma)),i,1,length(zeros));

modeshape : subst(abbrev[2],subst([x=xi*l[0]],subst(sol[1][4],subst(makelist(W[i]=EigVec[i][1],i,1,4),
                       subst(params,subst([t=T[ref]*tau],subst(abbrev,subst(ansatz,w(x,t)))))))));
modes : makelist(subst(zeros[i],modeshape),i,1,length(zeros));

/* plot mode snapshots for animation -> convert with imageMagick */
fpprintprec: 3$
titles : makelist(simplode(["ω/ω1 = ",subst(sol[2][1],subst(sol[2][i],omega[0])/omega[0])]),i,1,6);
NSteps : 19$
for mode: 1 thru 6 do
   for j:0 thru NSteps do
      (tstep : simplode([if j<10 then "0" else "",j]),
       plot2d(subst([tau=j/(NSteps+1)],subst([tau=0],modes[mode])*sin(2*%pi*tau)),[xi,0,1],
                        [title, titles[mode]], [legend, simplode(["t/T=",float(j/NSteps+0.01)])],
                        [style, [lines,3,1]], [xlabel, "ξ →"], [ylabel, "← w "],
                        [png_file,simplode([workdir,"mode-",mode,"-step-",tstep,".png"])],
                        [gnuplot_preamble, "set yrange [2.1:-2.1] reverse"]))$

/* plot all modeshapes (1...6)*/
legends : append([legend],makelist(simplode(["mode ",i]), i,1,length(zeros)));
plot2d(subst([tau=0],modes), [xi,0,1], legends,
             [xlabel, "ξ →"], [ylabel, "← w "],
             [gnuplot_preamble, "set yrange [2.1:-2.1] reverse"])$




Adapt to Initial Condition

Zum Zeitpunkt t=0 soll

gelten. Das kriegen wir hin, indem wir die (unendlich vielen) Ci der homogenen Lösung passend wählend. Das klingt komplizierter als es ist - die Mathematik bereitet uns den Weg.

Die Beziehung oben multiplizieren wir dazu mit einer der Modalformen ϕk und bilden das Integral über die Balken-Lange , also

Wir projizieren also die Modalformen auf sich selbst und auf die partikulare Lösung. Und nutzen aus, dass diese Faltungs-Integrale der Modalformen eine besondere Eigenschaft haben:

Man sagt: die Funktionen sind zueinander verallgemeinert orthogonal. Damit berechnen wir

,

Nur die erste Modalform liefert einen signifikanten Beitrag bei der Anpassung an die Anfangsbedingungen, wir nehmen bei der Auftragung der Gesamtlösung wt(x,t) trotzdem die ersten sechs Modalformen mit:

Animation: Loslassen des Systems aus der unverformten Rugelage.

In der Animation sehen Sie das Loslassen des Stabes aus der unverformten, geraden Referenzlage für genau eine Periodenlänge T1 der ersten Modalform. Allerdings ist die Schwingung, die entsteht, nicht T1-periodisch - die anderen Modalformen leisten einen (kleinen) Beitrag, dessen Schwingungsperiode kein Vielfaches der Grundschwingung ist. Das sehen sie an dem kleinen "Ruck" in der Darstellung der Lösung, wenn nach der Dauer T1 der Grundschwingung wieder die Anfangslage eingeblendet wird.


/***************************************************************************/
/* t o t a l   ( c o m p o s e d )   s o l u t i o n                       */
/***************************************************************************/

w[h] : sum(C[i]*modes[i],i,1,length(zeros));

/* adapt to inital condition w[h]+w[p]=0  */
/* Faltungsintegrale */
FI : [integrate(expand(subst([tau=0],modeshape^2)),xi,0,1),
      integrate(expand(subst([tau=0],modeshape*ratsimp(subst(params,subst(w[p],w(x,t))/w[s])))),xi,0,1)];
/* coefficients for mode shapes */
IVs : makelist(C[i] = -subst(zeros[i],FI[2]/FI[1]),i,1,length(zeros));

w[h] : subst(IVs,realpart(w[h]));

/* produce slides for animation */
fpprintprec: 3$
NSteps : 39$
for j:0 thru NSteps-1 do
   (tstep : simplode([if j<10 then "0" else "",j]),
    plot2d(subst([tau=j/NSteps],w[h])+ratsimp(subst(params,subst(w[p],w(x,t))/w[s])),[xi,0,1], 
               [legend, simplode(["t/T=",float(j/NSteps+0.01)])],
               [style, [lines,3,1]], [xlabel, "ξ →"], [ylabel, "← w"],
               [png_file,simplode([workdir,"initial-value-solution-step-",tstep,".png"])],
               [gnuplot_preamble, "set yrange [2.5:-0.5] reverse"]))$




🧨 Timoshenko-Balken:
Wieviele Moden man berücksichtigen muss, hängt von vielen praktischen Faktoren - z.B. der modalen Dämpfung - ab. Vergessen sollte man nicht, dass je nach Stab-Länge, der Einfluss der Schubverformung bei den höheren Moden zunimmt. Dann muss man ggf. das Modell wechseln ...
Saiteninstruente:
Bei Saiteninstrumenten möchte man Biegeschwingungen der Saiten verhindern. Darum sind besonders die Saiten, die eine hohe Massebelegung brauchen, um in niedrigen Frequenzen zu schwingen, besondere Konstruktionen. Dabei wird oft ein Draht um die eigentliche Saite gewickelt. Können Sie aus diesen Untersuchungen erkennen, was das Ziel dieser Konstruktion ist?

Links

  • ...

Literature

  • ...