Gelöste Aufgaben/FEC1/FEC1 - Teil III: Simulieren

⬅ Zurück zu Teil II

Computerprogramm schreiben

Um unsere Bewegungsgleichungen

${\underline {\underline {M}}}\;{\underline {\ddot {Q}}}+{\underline {\underline {G}}}\;{\underline {\dot {Q}}}+\left({\underline {\underline {K}}}_{C}\cos(2\Omega t)+{\underline {\underline {K}}}_{S}\sin(2\Omega t)+{\underline {\underline {K}}}_{0}\right)\;{\underline {Q}}={\underline {P}}$

zu lösen, war der Computer und Maxima bis hierher eine gute Hilfe - wir hätten aber auch alle Gleichungen per Hand aufstellen können. Ab hier geht das nicht mehr. Für die Lösung der Bewegungsgleichungen brauchen wir Computer. Für einfache Probleme wie dieses reichen Interpreter-Programme wie Matlab©, wenn's komplizierter wird, müssten wir auf Kompiler wie C++ umsteigen. Für alle etablierten Programme gibt es Numerik-Bibliotheken, um die Probleme zu lösen, die wir jetzt angehen.

Bei allen Lösungsansätzen steht in der Bewegungsgleichung noch ein Parameter: die Dimensionslose Drehzahl η mit

$\displaystyle \Omega =\eta \cdot \omega _{0,0}{\text{ und der Eigenfrequenz }}\omega _{0,0}^{2}=3{\frac {EI}{\ell ^{3}m_{I}}}$

Partikulare Lösung für konstante Koeffizienten

Für die Bewegungsgleichung oben gibt es keine geschlossene partikulare Lösung - die periodischen Koeffizienten machen das unmöglich. Anschaulich berücksichtigen diese Koeffizienten die Drehwinkel-veränderliche Steifigkeit der Lager beim Drehen. 

Wir könnten von Δk = 0 ausgehen und die Lösung von

${\underline {\underline {K}}}_{0}\;{\underline {Q}}_{P}={\underline {P}}$

berechnen. Aber hier steht nichts, was wir nicht mit einer Handrechnung auch hinbekommen würden: die Auslenkung des Systems wächst quadratisch mit der Drehzahl η(Ω) - also der Zentrifugalkraft.

tmp

Homogene Lösung für konstante Koeffizienten

Auch für die homogene, lineare Differentialgleichung (für P=0) gibt es keine geschlossene analytische Lösung. Wir schauen und also auch hier erst einmal nach einer Lösung der homogenen Differentialgleichung mit konstanten Koeffizienten

${\underline {\underline {M}}}\;{\underline {\ddot {Q}}}+{\underline {\underline {G}}}\;{\underline {\dot {Q}}}+{\underline {\underline {K}}}_{0}\;{\underline {Q}}={\underline {0}}$

um und führen wir eine Modalanalyse durch. Damit wir mit Standard-Lösungsroutinen arbeiten können, transformieren wir die Bewegungsgleichung in eine Differentialgleichung erster Ordnung, hier

$\left({\begin{array}{c}{\underline {\dot {P}}}\\{\underline {\dot {Q}}}\end{array}}\right)=\left({\begin{array}{cc}-{\underline {\underline {M}}}^{-1}\cdot {\underline {\underline {G}}}&-{\underline {\underline {M}}}^{-1}\cdot {\underline {\underline {K}}}\\{\underline {\underline {1}}}&{\underline {\underline {0}}}\end{array}}\right)\cdot \left({\begin{array}{c}{\underline {P}}\\{\underline {Q}}\end{array}}\right)$.

Dabei ist P = dQ/dτ, also die Geschwindigkeit von Auslenkung und Kippung. Wir haben also nun eine neue Schreibweise für die homogene Differentialgleichung als

${\underline {\dot {X}}}={\underline {\underline {H}}}(\eta )\cdot {\underline {X}}$

mit

${\underline {X}}=\left({\begin{array}{c}{\underline {P}}\\{\underline {Q}}\end{array}}\right){\text{ und }}{\underline {P}}={\underline {\dot {Q}}}$.

Für dieses Standard-Problem können wir in Maxima eine Adaption der Lapack-Routinen nutzen.

/*******************************************************/
/******************* PART III **************************/
/*******************************************************/
/* rewrite as first-order differential equations       */
/* X' = H(τ)*X + R, X = (P,Q)^T, R = (0,P)^T           */

/* inverse of  M */
I: dgesv (mats[1], diagmatrix (8, 1))$
/* System Matrix */
H  : zeromatrix(2*8,2*8)$
/* G*M^(-1) */
GI : expand(mats[2].I)$

/* populate system matrix */
for row:1 thru 8 do
   for col:1 thru 8 do
      (H[  row,8+col] : + I     [row,col],
       H[8+row,  col] : -mats[3][row,col],
       H[8+row,8+col] : -GI     [row,col])$

/*******************************************************/
/* solve eigenvalue Problem                            */

evals: []$
for Eta:0 thru 30 step 0.1 do
    (/*print("eta = ", Eta),*/
     evals: append(evals, [append([Eta], dgeev (subst([eta = Eta],H))[1])]))$

etas : flatten(makelist(makelist(evals[i][1],j,1,length(evals[i])-1),i,1,length(evals)))$
reals: flatten(makelist(makelist(realpart(evals[i][j]),j,2,length(evals[i])),i,1,length(evals)))$
imags: flatten(makelist(makelist(imagpart(evals[i][j]),j,2,length(evals[i])),i,1,length(evals)))$
plot2d([[discrete, etas, reals],[discrete, etas, imags]], [y,-100,100], 
        [style,points,points], [point_type,bullet,plus],
        [color,red, blue], [legend, "real part", "imag part"],
        [xlabel, "η->"], [ylabel, "ω, γ ->"],
        [gnuplot_term, "png size 2000,2000 font 'Helvetica,30'"], [gnuplot_out_file, "C:\\Users\\abs384\\OneDrive\\X\\tmp\\plots\\eigenvalues.png"] )$