Gelöste Aufgaben/FEAE

Aufgabenstellung

Wir schauen uns die Lösungsanteile für die Bewegungsgleichung des Zwei-Masse-Schwinger an: welche Rolle spielen die homogene und partikulare Lösung?

Lageplan

Gesucht ist Gesamtlösung des Systems im Zeitbereich beim Loslassen der beiden Massen aus der Referenz-Konfiguration, bei der beide Federn entspannt sind. Wir stellen die Bewegungsgleichung mit dem Prinzip der virtuellen Verrückungen auf.

Wir arbeiten mit den System-Parametern

$\begin{array}{ccc}{k}_1& =& \frac{3}{2}{k}_0\\ k_2& =& k_0\\ m_1& =& m_0\\ m_2& =& m_0\\ \end{array}$

Lösung mit Maxima

Header

Der Lösungsweg zu dieser Aufgabe wird bereits in Integration der Differentialbeziehung (IVP) gezeigt, hier verfolgende wir einen etwas allgemeineren Lösungspfad.

Dabei arbeiten wir u.a. mit quadratischen Gleichungen. Für deren Lösung müssen wir Maxima wissen lassen, welche Parameter positiv sind, also

$k_0>0,m_0>0,T>0;$

Koordinaten der Bewegung.

Die Lage-Koordinaten des Systems fassen wir in

$\underset{\_}{Q}=\left(\begin{array}{l}{w}_1(t)\\ w_2(t)\end{array}\right)$

zusammen.

Für die dimensionslose Schreibweise brauchen wir später eine Bezugszeit t_Bez.

Dazu nehmen wir eine "Anleihe" bei einem ähnlichen System:

Für dieses System "kennen" wir die Eigenkreisfrequenz

${\omega }_0=\sqrt{\frac{k_0}{m_0}}\text{ und }{\omega }_0=\frac{2\pi }{T}$.

Aus diesem Modell "leihen" wir uns die Bezugszeit zu

$T=2\pi \sqrt{\frac{m_0}{k_0}}$.

und wählen

$t_{Bez}:=T$

Hier kennen wir ω₀.

/*******************************************************/
/* MAXIMA script                                       */
/* version: wxMaxima 16.04.2                           */
/* author: Andreas Baumgart                            */
/* last updated: 2018-01-07                            */
/* ref: FEM, PVMPE using two coordinates               */
/* description: solve by principle virt. Work          */
/*******************************************************/
/* declare variational variables - see 6.3 Identifiers */
declare("δW", alphabetic);
declare("δA", alphabetic);
declare("δΠ", alphabetic);
declare("δw", alphabetic);
declare("δQ", alphabetic);

/* assumptions */
assume(k[0]>0,m[0]>0,T>0);

/* abbreviations */
abbrev : [omega[0,0]=2*%pi/T, omega[0,0]^2= k[0]/m[0]];
abbrev : append(abbrev, solve(abbrev,[omega[0,0],T])[2]);

/* system parameters and dimensionless time tau */
params: [m[1]=m[0],m[2]=m[0],k[1]=3/2*k[0],k[2]=k[0],
         t = T*tau, w[s] = m[0]*g/k[0]];

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

Equilibrium Conditions

Die Gleichgewichtsbedingung mit dem Prinzip der virtuellen Verrückungen lautet

$\begin{array}{lcl}\delta W& =& \underset{:=\delta W^a}{\underbrace{-m_1\cdot {\ddot{w}}_1\cdot \delta w_1-m_2\cdot {\ddot{w}}_2\cdot \delta w_2+m_1\cdot g\cdot \delta w_1+m_2\cdot g\cdot \delta w_2}}\\ & & -\underset{:=\delta \mathrm{\Pi }}{\underbrace{k_1\cdot w_1\cdot \delta w_1+k_2\cdot (w_2-w_1)\cdot (\delta w_2-\delta w_1)}}\\ & \stackrel{!}{=}0\end{array}$.

Aufgelöst nach den Koeffizienten der virtuellen Verrückungen folgt die gewöhnliche, lineare Differentialgleichung:

$\left(\begin{array}{cc}{m}_0& 0\\ 0& m_0\end{array}\right)\cdot \left(\begin{array}{c}{\ddot{w}}_1\\ {\ddot{w}}_2\end{array}\right)+\left(\begin{array}{cc}\frac{5}{2}{k}_0& -k_0\\ -k_0& k_0\end{array}\right)\cdot \left(\begin{array}{c}{w}_1\\ w_2\end{array}\right)=\left(\begin{array}{c}{m}_0\cdot g\\ m_0\cdot g\end{array}\right)$.

In dieser Beziehung stehen jetzt nur noch m₀ und k₀ - die Indizes "0" können wir also ab hier weglassen. In Matrix-Schreibweise lautet die Bewegungsgleichung jetzt:

$\underset{\_}{\underset{\_}{M}}\cdot \underset{\_}{\ddot{Q}}+\underset{\_}{\underset{\_}{K}}\cdot \underset{\_}{Q}=\underset{\_}{G}$.

/* virtual of system (equilibrium condition) */
PvV : [δW = δW^a - δΠ,
       δW^a = - sum(m[i]* diff(w[i](t),t,2)*δw[i],i,1,2) + sum(m[i]*g*δw[i],i,1,2),
       δΠ   =   sum(k[i]*(w[i](t)-w[i-1](t))*(δw[i]-δw[i-1]),i,1,2)];

/* equlilbrium condition */

eom: subst(PvV[3],subst(PvV[2],subst(PvV[1],δW=0)));
eom: expand(subst([w[0](t)=0,δw[0]=0],eom));

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));

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

Solving

Die Lösung des Anfangswertproblems setzt sich aus zwei Teilen zusammen:

• der partikularen Lösung, die die Rechte Seite "G" erfüllt und
• der homogenen Lösung, die die Rechte Seite "0" erfüllt.

Die Gesamtlösung Q_t setzt sich nun - bei diesem linearen System - additiv aus partikularer Q_p und Q_h homogener Lösung zusammen:

${\underset{\_}{Q}}_t(t)={\underset{\_}{Q}}_p(t)+{\underset{\_}{Q}}_h(t)$

/**************** total solution */
Q[t] : Q[p] + Q[h]

Right-Hand-Side Approach

Die rechte Seite der Bewegungsgleichung G ist nicht zeitabhängig - sie ist statisch. Also ist auch die Lösung Q_p - die partikulare Lösung - statisch, wir suchen nach der Lösung des Gleichungssystems

$\underset{\_}{\underset{\_}{K}}\cdot {\underset{\_}{Q}}_p=\underset{\_}{G}$

und die ist

${\underset{\_}{Q}}_p=\left(\begin{array}{c}\frac{4mg}{3k}\\ \frac{7mg}{3k}\end{array}\right)$

Wir führen die Abkürzung

$w_s=\frac{mg}{k}$

ein, also ist

${\underset{\_}{Q}}_p=w_s\cdot \left(\begin{array}{c}\frac{4}{3}\\ \frac{7}{3}\end{array}\right)$.

/**************** particular solution */
Q[p] : ratsimp(subst(params,linsolve_by_lu(K,G)[1]));

Homogeneous Solution

Für die Lösung der homogenen Bewegungsgleichung

$\underset{\_}{\underset{\_}{M}}\cdot \underset{\_}{\ddot{Q}}+\underset{\_}{\underset{\_}{K}}\cdot \underset{\_}{Q}=\underset{\_}{0}$

setzen wir (vgl. Modalanalyse) an

${\underset{\_}{Q}}_h={\underset{\_}{\hat{Q}}}_h\cdot e^{\lambda \cdot t}$

und erhalten das Eigenwertproblem

$\underset{\underset{\_}{\underset{\_}{D}}(\lambda )}{\underbrace{({\lambda }^2\underset{\_}{\underset{\_}{M}}+\underset{\_}{\underset{\_}{K}})}}\cdot {\underset{\_}{\hat{Q}}}_h=\underset{\_}{0}$

mit

den Eigenwerten	$\lambda$ und
den Eigenvektoren	${\underset{\_}{\hat{Q}}}_h$

Die Berechnung der Eigenwerte λ wird einfacher mit

${\lambda }^2=-\mathrm{\Lambda }$

und wir transformieren parallel die Bewegungsgleichungen auf die dimensionslose Zeit

$\tau :=\frac{t}{t_{Bez}}$ mit $t_{Bez}=T.$,

also

$\underset{\_}{\ddot{Q}}=\frac{d^2\underset{\_}{Q}}{d{\tau }^2}\cdot \underset{\frac{1}{T^2}}{\underbrace{\frac{d^2\tau }{d{t}^2}}}$.

Wir finden zwei Eigenwerte aus der Bedingung det(D)=0:

$\begin{array}{lr}{\mathrm{\Lambda }}_1=& 2\cdot {\pi }^2\\ {\mathrm{\Lambda }}_2=& 12\cdot {\pi }^2\\ \end{array}$

und

$\begin{array}{lr}{\lambda }_1=& j\cdot \sqrt{2}\cdot \pi \\ {\lambda }_2=& j\cdot 2\sqrt{(}3)\cdot \pi \\ \end{array}\text{ mit }j:=\sqrt{-1}$.

Für die Matrix

$\underset{\_}{\underset{\_}{D}}({\lambda }_i)={\lambda }^2\cdot \frac{1}{T^2}\cdot \underset{\_}{\underset{\_}{M}}+\underset{\_}{\underset{\_}{K}}$

stellen wir fest:

Die Matrix hat

• einen Rang (rank) von 1 - es gibt eine linear abhängige Zeile im Gleichungssystem und - entsprechend -
• einen Rangabfall (nullity)  von 1.

Die Eigenvektoren spannen nun den Nullraum (nullspace) der Matrix auf. Normmert auf die Länge 1 sind das

• für λ₁:
  ${\underset{\_}{\hat{Q}}}_{h,1}=\frac{1}{\sqrt{5}}\cdot \left(\begin{array}{c}1\\ 2\end{array}\right)$ und
• für λ₂:
  ${\underset{\_}{\hat{Q}}}_{h,2}=\frac{1}{\sqrt{5}}\cdot \left(\begin{array}{c}2\\ -1\end{array}\right)$.

Die homogene Lösung der Bewegungsgleichung lautet damit

${\underset{\_}{Q}}_h(t)=c_1\cdot \frac{1}{\sqrt{5}}\cdot \left(\begin{array}{r}1\\ 2\end{array}\right)\cdot e^{j\cdot \sqrt{2}\pi \cdot \tau }+c_2\cdot \frac{1}{\sqrt{5}}\cdot \left(\begin{array}{r}2\\ -1\end{array}\right)\cdot e^{j\cdot 2\sqrt{3}\pi \cdot \tau }$

mit den Integrationskonstanten c₁ und c₂. Durch die komplexen Eigenwerte sind die beiden Integrationskonstanten nun auch komplexwertig, also

$\begin{array}{l}{c}_1=c_{R,1}+j\cdot c_{I,1}\\ c_2=c_{R,2}+j\cdot c_{I,2}\end{array}$.

Diese vier reell-wertigen Integrationskonstanten bekommen wir aus vier Anfangsbedingungen für die beiden Massen, hier

$\begin{array}{l}{w}_1(0)=0\\ w_2(0)=0\\ {\dot{w}}_1(0)=0\\ {\dot{w}}_2(0)=0\end{array}$.

Wir finden

$c_{R,1}=-\frac{6mg}{\sqrt{5}k},c_{I,1}=0,c_{R,2}=-\frac{mg}{3\sqrt{5}k},c_{I,2}=0$.

und damit

${\underset{\_}{Q}}_h(t)=15\cdot \left(\begin{array}{c}-2\cos (2\sqrt{3}\pi \cdot \tau )-18\cos (\sqrt{2}\pi \cdot \tau )\\ \cos (2\sqrt{3}\pi \cdot \tau )-36\cos (\sqrt{2}\pi \cdot \tau )\end{array}\right)$.

/* direct approach */
load(eigen); eigensystem: uniteigenvectors(subst(params, invert(M).K));
/* -> not employed */

/**************** homogeneous solution */
/* solve eigenvalue problem */
Q[E] : matrix([q[1,i]],[q[2,i]]);
ansatz: Q[h] = Q[E]*%e^(%lambda[i]*tau);

D : -Lambda*(1/T)^2*M + K;
D : subst(params, D);
charPoly : determinant(D);

/* eigenvalues */
eigval : solve(subst(params,charPoly)=0,Lambda);
eigval : subst(abbrev, eigval);

/* save for later ....*/
eigvalList : makelist(subst(eigval[j],[%lambda[j]=%i*sqrt(Lambda)]),j,1,2);

/* eigenvectors */
prop: [rnk = rank(subst(eigval[1],subst(params,D))), nly = nullity(subst(eigval[1],D))];

eigvec : makelist(apply (addcol, args (nullspace(
                           subst(abbrev,subst(eigval[i],D/k[0]))))),i,1,2);

/* norm eigenvectors (to be of length "1" */
eigvec : makelist(eigvec[i]/sqrt(eigvec[i].eigvec[i]),i,1,2);

/* adapt to initial values w[1](0)=0, w[2](0)=0 */
Q[h] : sum(subst(eigval[j],subst([%lambda[j]=%i*sqrt(Lambda),i=j],
                           (c[R,j]+%i*c[I,j])*eigvec[j]*%e^(%lambda[j]*tau))),j,1,2);

/* update Q[t]*/
Q[t] : ''(Q[t]);

/* adapt to initial values */
intcons : [c[R,1],c[I,1],c[R,2],c[I,2]];
intcond : append(makelist(subst([tau=0],     Q[t]     )[i][1]=0,i,1,2),
                 makelist(subst([tau=0],diff(Q[t],tau))[i][1]=0,i,1,2));

/* integration constants: */
sol[1] : ratsimp(solve(realpart(intcond),intcons)[1]);

Q[t] : realpart(ratsimp(subst(params,subst(sol[1],Q[t]/w[s]))))

Post-Processing

Die Lösung für

• w₁ (blau) und
• w₂ (rot)

im Zeitbereich, angepasst an die Anfangsbedingungen, sieht nun so aus:

Ergebnis-Plots für w₁(t), w₂(t).

Interessant ist die separate Auftragung der Lösungsanteile nach partikularer Lösung (Index "p") und den beiden homogenen Lösungsanteilen (Index "h1", "h2"). Man erkennt, wie die Summe jeweils der blauen und der roten Anteile die Gesamtlösung ergibt.

Ergebnis-Plot mit den Lösungsanteilen.

/* plot results */
preamble: "set yrange [] reverse";
plot2d([Q[t][1][1],Q[t][2][1]],[tau,0,4], [legend, "w[1]", "w[2]"],
    [xlabel, "τ/1 →"], [ylabel, "← w/(m*g/k)"],
    [gnuplot_preamble, preamble]);

/* deconstruct solution into modes */
timeFcts : makelist(subst(eigvalList[i],%e^(%lambda[i]*tau)),i,1,2)
Q[l] : subst(params,append([Q[p]/w[s]],
            realpart(makelist(coeff(subst(sol[1],Q[h]),timeFcts[i])/w[s]*timeFcts[i],i,1,2))));

plot2d([Q[l][1][1][1],Q[l][1][2][1],
        Q[l][2][1][1],Q[l][2][2][1],
        Q[l][3][1][1],Q[l][3][2][1]],
        [tau,0,4],
        [legend, "w[p,1]", "w[p,2]", "w[h1,1]", "w[h1,2]", "w[h2,1]", "w[h2,2]"],
        [color, blue, red, blue, red, blue, red],
        [style, [lines,3], [lines,3], [lines,1], [lines,1], [lines,1], [lines,1]],
        [xlabel, "τ/1 →"], [ylabel, "← w/(m*g/k)"],
        [gnuplot_preamble, preamble]);

Links

• ...

Literature

• ...