Gelöste Aufgaben/FEAG

Aufgabenstellung

Analog zu FEAF untersuchen wir hier die Schwingungen eines Kontinuums beim Loslassen aus der entspannten Rugelage. Hier nicht mit einem Dehnstab, sondern einem Euler-Bernoulli-Balken.

Lageplan

Gesucht ist die Schwingung eines Euler-Bernoulli-Balkens beim Loslassen aus der Ruhelage. Wir gehen nach dem Standardrezept der Finite Elemente Methode vor, arbeiten also mit dem Prinzip der virtuellen Arbeit.

Lösung mit Matlab

Interessant ist hier, dass - im Gegensatz zu Stablängsschwingungen - die Eigenfrequenz nicht ein gerades Vielfaches der untersten Eigenfrequenz ist. Falls Sie ein Saiteninstrument spielen, verstehen Sie sofort, warum das wichtig ist.

Maxima können wir hier nicht gut gebrauchen: die Gleichungen werden zu umfangreich. Wir arbeiten also mehr mit numerischen Verfahren, da ist Matlab geeigneter. Allerdings können wir Matlab-Inhalte nicht gut auf dieser Seite unterbringen - deshalb gibt es dafür die Seite FEAG-Matlab, die der gleichen Struktur folgt.

Header

Im Programm arbeiten wir mit einer dimensionslosen Formulierung - wir brauchen dafür eine Bezugszeit t_Bez und eine Bezugslänge l_Bez.

Bezugsgrößen wählen

Dazu nehmen wir eine "Anleihe" bei der analytischen Lösung des Schwingungsproblems (vgl. Aufgabe SKEB):

Analytische Lösung: homohener Lösungsanteil

Analytische Lösung: partikularer Lösungsanteil

Die homogene Bewegungsgleichung des Euler-Bernoulli-Balkens $\varrho \,A\cdot {\ddot {w}}+E\,I\cdot {w}^{IV}=0$ hat Lösungen vom Typ $\displaystyle w(x,t)=\sum _{i=1}^{\infty }W_{i}\cdot e^{\displaystyle \kappa _{i}\cdot x}\cdot e^{\displaystyle \omega _{0,i}\cdot t}$ Wir bekommen zu jedem ω0 vier κ ${\begin{array}{l}\displaystyle \kappa _{1}=+{\frac {j\,\alpha }{\ell }},\\\displaystyle \kappa _{2}=-{\frac {\alpha }{\ell }},\\\displaystyle \kappa _{3}=-{\frac {j\,\alpha }{\ell }},\\\displaystyle \kappa _{4}=+{\frac {\alpha }{\ell }}\end{array}}\;\;{\text{ mit }}\;\;\displaystyle \alpha ^{4}={\frac {{{\ell }^{4}}\,A\,\rho }{E\,I}}\cdot {{\omega }_{0}^{2}}$ Das α ist eine praktische Abkürzung, hinter der wir ω0 verstecken. Anders als beim Dehnstab (SKER) finden wir hier keine analytische Beziehung, sondern nur die numerischen Beziehungen, für die unsere Randbedingungen erfüllt sind: ${\begin{array}{c}\alpha _{1}=1.875\ldots ,\\\alpha _{2}=4.694\ldots ,\\\alpha _{3}=7.854\ldots ,\\\vdots \end{array}}$ Die langsamste Eigenmode gehört zu α1 mit der Periodendauer $\displaystyle T^{*}=0.5688\ldots \pi \ell ^{2}{\sqrt {\frac {\varrho \,A}{E\,I}}}$. Also wählen wir $t_{Bez}:=T^{*}$.

Die zugeordnete inhomogene Bewegungsgleichung des Euler-Bernoulli-Balkens $E\,I\cdot {w}^{IV}=\varrho \,A\cdot g$ hat die partikulare Lösung $\displaystyle E\,I\,{w}_{p}\left(x\right)=A\,g\,\varrho \,\ell ^{4}\cdot \left({\frac {{\xi }^{4}}{24}}-{\frac {{\xi }^{3}}{6}}+{\frac {{\xi }^{2}}{4}}\right)$ Die statische Auslenkung am unteren Ende ist demnach $\displaystyle {{w}_{s}}={\frac {{{\ell }^{4}}\,A\,g\,\rho }{8\,E\,I}}$. Wir wählen $\ell _{Bez}:=w_{s}$.

System-Parameter des FEM-Modells

Für die Diskretisierung wählen wir als Anzahl der Finiten Elemente (Number Of Elements):

$\displaystyle I=4\;\;\;({\text{ im Matlab-Code: }}NOE)$

und damit je Element

$\displaystyle \ell _{i}={\frac {\ell _{0}}{I}}$.

Wir wählen die kubische Ansatzfunktionen aus dem Abschnitt Finite Elemente Methode je Element, also

${\begin{array}{crc}\phi _{1}=&\left(\xi -1\right)^{2}\cdot \left(1+2\cdot \xi \right)\\\phi _{2}=&\left(\xi -1\right)^{2}\cdot \xi \cdot \ell _{i}\\\phi _{3}=&\left(\xi ^{2}\cdot \left(3-2\cdot \xi \right)\right)\\\phi _{4}=&\left(\left(\xi -1\right)\cdot \xi ^{2}\cdot \ell _{i}\right)\end{array}}$

Für die numerisch Implementierung stört in dieser Darstellung das Element-spezifische ℓ_i - die Elementlänge. Wir behelfen uns, indem wir die Ansatzfunktionen schreiben als

$\displaystyle \underbrace {\left({\begin{array}{c}\phi _{1}\\\phi _{2}\\\phi _{3}\\\phi _{4}\end{array}}\right)} _{\displaystyle =:{\underline {\phi }}}=\underbrace {\left({\begin{array}{cccc}1&0&0&0\\0&\ell _{i}&0&0\\0&0&1&0\\0&0&0&\ell _{i}\\\end{array}}\right)} _{\displaystyle =:{\underline {\underline {d}}}}\cdot \underbrace {\left({\begin{array}{c}\left(\xi -1\right)^{2}\cdot \left(1+2\cdot \xi \right)\\\left(\xi -1\right)^{2}\cdot \xi \\\left(\xi ^{2}\cdot \left(3-2\cdot \xi \right)\right)\\\left(\left(\xi -1\right)\cdot \xi ^{2}\right)\end{array}}\right)} _{\displaystyle =:{\underline {\tilde {\phi }}}}$,

wobei wir die Ansatz-Polynome in ${\underline {\tilde {\phi }}}$ verpacken.

Das schaut unnötig komplex aus - allerdings stecken wir die Diagonal-Matrizen d (für jedes Element eine) in eine Matlab-Variable, so dass sie dort nicht weiter auffällt.

Die abhängigen Koordinaten des FEM-Modells sind dann zunächst - bis zum Einarbeiten der Randbedingungen - diese:

${\underline {Q}}_{t}=\left({\begin{array}{c}W_{0}(t)\\\Phi _{0}(t)\\W_{1}(t)\\\vdots \\\Phi _{I}(t)\\W_{I}(t)\\\end{array}}\right)$.

Equilibrium Conditions

Die Gleichgewichtsbedingung mit dem Prinzip der virtuellen Verrückungen setzen sich additiv aus

$\displaystyle \delta W=\sum _{i=1}^{I}\left(\delta W_{i}^{a}-\delta \Pi _{i}\right)+\delta W_{R}^{a}$,

zusammen - also den virtuellen Arbeiten je Element zuzüglich von virtuellen Arbeiten am Rand. In unserem Beispiel haben wir allerdings keine eingeprägten, äußeren Lasten am Rand, also ist

$\delta W_{R}^{a}=0$.

Je Element erfassen wir die virtuellen Arbeiten durch Element-Matrizen. So ist z.B. für die virtuelle Formänderungsenergie (vgl. Finite Elemente Methode)

$\delta \Pi _{i}=\left(\delta W_{i-1},\delta \Phi _{i-1},\delta W_{i},\delta \Phi _{i}\right)\cdot {\underline {\underline {K}}}_{i}\cdot \left({\begin{array}{c}W_{i-1}(t)\\\Phi _{i-1}(t)\\W_{i}(t)\\\Phi _{i}(t)\end{array}}\right)$

mit der Element-Steifigkeitsmatrix ${\underline {\underline {K}}}_{i}$.

Aus der FEM-Formulierung für den Euler-Bernoulli-Balken könnten wir die Element-Matrizen der virtuellen Arbeit der D'Alembert'schen Trägheitskraft sowie der Virtuelle Formänderungsenergie herauskopieren. Auch die "rechte Seite" mit der virtuelle Arbeit der Gewichtskraft aus

$\displaystyle \delta W_{i}^{a,g}=\int _{\ell _{i}}\varrho \;A\;g\cdot \delta w(x)\;dx_{i}$,

könnten wir als Element-Lastmatrix (Spaltenmatrix) abspeichern.

$\displaystyle {\underline {G}}_{i}={\frac {\varrho \,A\,g\,\ell _{i}}{2}}\cdot \left({\begin{array}{c}+1\\\displaystyle +{\frac {\ell _{i}}{6}}\\+1\\\displaystyle -{\frac {\ell _{i}}{6}}\end{array}}\right)$.

Dass es auch anders geht,sehen Sie in der Definition der Matlab-Klasse FEM_Element_Modell, bei der die Element direkt aus den Trial-Functions - hier den ${\underline {\tilde {\phi }}}$ -

hergeleitet. werden:

                 :
          K(row,col) = diff(...
                            polyval(...
                                 polyint(...
                                     conv(polyder(self.phi(row,:)),polyder(self.phi(col,:)))),...
                                     [0,1]));
                 :

Komponieren der System-Matrizen

Die System-Matrizen des Gesamt–Systems komponieren wir nun durch Hinzuaddieren der Anteile je Element.

Beispiel: die Gesamt-Steifigkeitsmatrix für I=2:

$\displaystyle {\underline {\underline {K}}}={\frac {E\,I}{\ell _{i}^{3}}}\cdot \left({\begin{array}{ccccc}{\color {red}{+12}}&{\color {red}{+6}}&{\color {red}{-12}}&{\color {red}{+6}}&0&0\\{\color {red}{+6\ell _{i}}}&{\color {red}{+4\ell _{i}^{2}}}&{\color {red}{-6\ell _{i}}}&{\color {red}{+2\ell _{i}^{2}}}&0&0\\{\color {red}{-12}}&{\color {red}{-6}}&{\color {red}{+12}}{\color {green}{+12}}&{\color {red}{-6}}{\color {green}{+6}}&{\color {green}{-12}}&{\color {green}{+6}}\\{\color {red}{+6\ell _{i}}}&{\color {red}{2\ell _{i}^{2}}}&{\color {red}{-6\ell _{i}}}{\color {green}{+6\ell _{i}}}&{\color {red}{+2\ell _{i}^{2}}}{\color {green}{+2\ell _{i}^{2}}}&{\color {green}{-6\ell _{i}}}&{\color {green}{+2\ell _{i}^{2}}}\\0&0&{\color {green}{-12}}&{\color {green}{-6}}&{\color {green}{+12}}&{\color {green}{-6}}\\0&0&{\color {green}{+6\ell _{i}}}&{\color {green}{2\ell _{i}^{2}}}&{\color {green}{-6\ell _{i}}}&{\color {green}{+2\ell _{i}^{2}}}\end{array}}\right)$

Die Beiträge der zwei Elemente sind hier in rot bzw. grün eingefärbt.

Dieses "Hinzuaddieren" passiert in Matlab hier:

                 :
             e = 0;% counter for all elements of all sections
             for section = 1: length(element.noe)
                            :
                 for i = 1:element.noe(section)
                     e = e+1;
                            :
                     j = 2*(e-1)+1;
                            :
                     K(j:j+3,j:j+3) = K(j:j+3,j:j+3)+b*  d(:,:,e)*element.K*d(:,:,e);
                            :
                 end    
                 :

In Matrix-Schreibweise lautet die Bewegungsgleichung des Gesamt-Systems jetzt:

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

Einarbeiten der Randbedingungen

Randbedingungen einarbeiten.

Die Randbedingungen arbeiten wir ein, indem wir zeilenweise (für δW und δΦ) und spaltenweise (für W und Φ) streichen.

Hier sind es

• die erste Zeile/Spalte für W₀ (blau) und
• die zweite Zeile / Spalte für Φ₀ (grün).

Solving

Die Lösung des Anfangs- und Randwertproblems ist in diesem Lösungsansatz mit Matlab um die Klasse "FEM_Container" herum aufgebaut - in ihr sind alle Parameter und Lösungsprozesse beschrieben.

Hier heißt die Instanz des Modells

• mathModel

In diesem Container sind alle Parameter und Zustandsgrößen des Modells bereits dimensionslos gemacht - und zwar mit den oben genannten Bezugsgrößen, die in der Excel-Eingabedatei im Blatt "0) Bezugsgrößen" definiert sind.

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:

${\underline {Q}}_{t}(t)={\underline {Q}}_{p}(t)+{\underline {Q}}_{h}(t)$

Particular Solution

Particulare Lösung.

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

${\underline {\underline {K}}}\cdot {\underline {Q}}_{p}={\underline {G}}$

und die ist - mit der Normierung durch w_s  -

 >> Qp
 
 Qp =
 
          0
          0
     0.1055
     0.0042
     0.3542
     0.0064
     0.6680
     0.0072
     1.0000
     0.0073

Diese Lösung tragen wir - elementweise - auf.

Homogene Lösung

Zur Lösung der homogenen Bewegungsgleichung

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

setzen wir an

$\displaystyle {\underline {Q}}_{h}={\underline {\hat {Q}}}_{h}\cdot e^{\lambda \cdot t}$

und erhalten das Eigenwertproblem

$\underbrace {\left(\lambda ^{2}{\underline {\underline {M}}}+{\underline {\underline {K}}}\right)} _{=:\displaystyle {\underline {\underline {D}}}}\cdot {\underline {\hat {Q}}}_{h}={\underline {0}}$

mit

den Eigenwerten	$\lambda$ und
den Egenvektoren	$\displaystyle {\underline {\hat {Q}}}_{h}$.

Für die Berechnung der Eigenwerte λ müssen wir die Abkürzung

$\lambda ^{2}=\Lambda$

einführen. Damit arbeitet die Matlab-Routine "eig()", die in "mathModel.eigensystem()" implementiert ist.

Acht Eigenwerte kommen aus der Bedingung $\det({\underline {\underline {D}}})=0$, diese $\Lambda _{i}$ werden auf der Spur der Matrix D abgelegt:

 >> D
 
 D =
 
    1.0e+06 *
 
    -2.9006         0         0         0         0         0         0         0
          0   -1.0774         0         0         0         0         0         0
          0         0   -0.4287         0         0         0         0         0
          0         0         0   -0.1662         0         0         0         0
          0         0         0         0   -0.0480         0         0         0
          0         0         0         0         0   -0.0123         0         0
          0         0         0         0         0         0   -0.0016         0
          0         0         0         0         0         0         0   -0.0000

Die zugehörigen dimensionslosen Periodendauern sind

${\begin{array}{ll}T_{1}=&0.0037\\T_{2}=&0.0061\\T_{3}=&0.0096\\T_{4}=&0.0154\\T_{5}=&0.0287\\T_{6}=&0.0566\\T_{7}=&0.1594\\T_{8}=&1.0000\end{array}}$.

Und offensichtlich fällt die längste Periodendauer T₈ mit der analytisch berechneten untersten Schwingunsperiode T^* zusammen.

Die zugehörigen Eigenvektoren stehen in V:

 >> V
 
 V =
 
    -0.0579    0.2951   -0.0687    0.3043    0.4561    0.5309   -0.3197   -0.0780
    -0.0129    0.0882   -0.0805   -0.0578   -0.0115    0.0064   -0.0096   -0.0032
    -0.0884   -0.1020    0.4873   -0.0417   -0.4707    0.0160   -0.5469   -0.2721
    -0.0364    0.1230    0.0134    0.0659   -0.0007   -0.0222    0.0019   -0.0051
    -0.1652   -0.4100   -0.2420   -0.1979    0.4058   -0.4268   -0.1035   -0.5270
    -0.0763    0.0330    0.0704   -0.0601    0.0098    0.0106    0.0163   -0.0059
    -0.9540   -0.8313   -0.8232    0.9212   -0.6357    0.7308    0.7663   -0.8013
    -0.2099   -0.1390   -0.0999    0.0809   -0.0389    0.0315    0.0200   -0.0060

Die Eigenvektoren in V spannen nun den Nullraum (nullspace) der Matrix auf, es ist

${\underline {\underline {V}}}=\left(\;{\underline {\hat {Q}}}_{h,1}\;,\;{\underline {\hat {Q}}}_{h,2}\;,\;{\underline {\hat {Q}}}_{h,3}\;,\;{\underline {\hat {Q}}}_{h,4}\;,\;{\underline {\hat {Q}}}_{h,5}\;,\;{\underline {\hat {Q}}}_{h,6}\;,\;{\underline {\hat {Q}}}_{h,7}\;,\;{\underline {\hat {Q}}}_{h,8}\;\right)$

Die homogene Lösung der Bewegungsgleichung lautet damit

$\displaystyle {\underline {Q}}_{h}(t)=\sum _{i=1}^{8}C_{i}\cdot {\underline {\hat {Q}}}_{h,i}\cdot e^{\displaystyle j\cdot \omega _{0,i}\cdot \tau }$

mit den acht Integrationskonstanten C_i . Durch die komplexen Eigenwerte sind die Integrationskonstanten nun (eigentlich) auch komplexwertig. Darum kommen in diesem Fall herum, weil die Anfangsgeschwindigkeit des Balkens beim Loslassen Null ist - wir also nur cos-Anteile berücksichtigen müssen. Und die gehören wiederum zum Realteil der Exponential-Funktion.

Die C_i sind nun die Konstanten, die wir brauchen, um die Lösung an Anfangsbedingungen anzupassen.

Vorher schauen wir uns die Lösung jeweils zu einer Eigenfrequenz ω_i an. Diese Funktionen heißen Modalformen ϕ(x) und deren Schwingungen können wir plotten:

Mode	Modalform ϕj(x)	Mode	Modalform ϕj(x)
#8 $\displaystyle \omega _{0,8}=1.000\cdot {\frac {2\pi }{T_{ref}}}$	Mode #1	#7 $\displaystyle \omega _{0,7}=6.2742\cdot {\frac {2\pi }{T_{ref}}}$	Mode #2
#6 $\displaystyle \omega _{0,6}=17.6833\cdot {\frac {2\pi }{T_{ref}}}$	Mode #3	#5 $\displaystyle \omega _{0,5}=34.8854\cdot {\frac {2\pi }{T_{ref}}}$	Mode #4
#4 $\displaystyle \omega _{0,4}=64.8852\cdot {\frac {2\pi }{T_{ref}}}$	Mode #5	#3 $\displaystyle \omega _{0,3}=104.2059\cdot {\frac {2\pi }{T_{ref}}}$	Mode #6

Modalformen

Alle acht Moden ϕ_j können wir auch zum Zeitpunkt τ=0 übereinander darstellen:

✔ Zur Euler-Bernoulli-Hypothese:

Praktisch haben die höheren Moden kaum Relevanz - meist klingen sie durch Dämpfung schnell ab. Man sieht allerdings bereits an der Schwingungs-Form, dass hier die Länge zwischen zwei Knoten nicht mehr sehr klein ist im Vergleich zur Höhe. Und das führt dazu, dass Schubverformungen eine Rolle spielen. Wir müssen zum Timoshenko-Balken wechseln ....

Adapt to Initial Conditions

Die acht reell-wertigen Integrationskonstanten C_i bestimmen wir aus den Anfangsbedingungen für den Balken, nämlich

${\underline {Q}}_{t}(0)={\underline {0}}{\text{ und }}{\underline {\dot {Q}}}_{t}(0)={\underline {0}}$ ,

also

${\begin{array}{llcll}W_{1}(0)&=0&\;\;\;&{\dot {W}}_{1}(0)&=0\\\Phi _{1}(0)&=0&&{\dot {\Phi }}_{1}(0)&=0\\W_{2}(0)&=0&&{\dot {W}}_{2}(0)&=0\\\Phi _{2}(0)&=0&&{\dot {\Phi }}_{2}(0)&=0\\W_{3}(0)&=0&&{\dot {W}}_{4}(0)&=0\\\Phi _{3}(0)&=0&&{\dot {\Phi }}_{3}(0)&=0\\W_{4}(0)&=0&&{\dot {W}}_{4}(0)&=0\\\Phi _{4}(0)&=0&&{\dot {\Phi }}_{4}(0)&=0\\\end{array}}$.

Wenn wir die C_i komplexwertig denken, so stehen hier 16 Anfangsbedingungen für acht komplexe C_i .

Denken wir die C_i gleich reelwertig, so können wir die acht Anfangsbedingungen in der rechten Spalte weglassen - sie sind dann implizit mit erfüllt - für unser spezielles Anfangswertproblem.

Die Anfangsbedingung für die Gesamtlösung lautet nun:

${\underline {Q}}_{p}+{\underline {Q}}_{h}|_{\tau =0}={\underline {0}}$ ,

also

${\underline {\underline {V}}}\cdot \underbrace {\left({\begin{array}{c}C_{1}\\C_{2}\\C_{3}\\C_{4}\\C_{5}\\C_{6}\\C_{7}\\C_{8}\end{array}}\right)} _{\displaystyle =:{\underline {C}}}=-{\underline {Q}}_{p}$

Wir finden

${\underline {C}}=\left({\begin{array}{l}+6.87897955904839E^{-07}\\-4.64941679695326E^{-06}\\+1.66295172107590E^{-05}\\+4.37871009978632E^{-05}\\+0.000288375976158946\\+0.00145539720136673\\-0.0186576428599944\\-1.26472054045511\end{array}}\right)$ .

und damit

$\displaystyle {\underline {Q}}_{t}(t)={\underline {Q}}_{p}+\sum _{i=1}^{8}C_{i}\cdot {\underline {\hat {Q}}}_{h,i}\cos {\omega _{0,i}\cdot \tau }$ .

Offensichtlich ist der Beitrag der Mode 8 (also der langsamsten Mode) bei weitem am größten Anteil an Q_h. Das können wir auch in der Animation der Lösung - also das Loslassen des Balken aus der Ruhe heraus - sehen.

Post-Processing

Die FEM-Lösung im Zeitbereich, angepasst an die Anfangsbedingungen, sieht nun so aus:

Lösung im Zeitbereich - hier für I=4.

Animation in Time Domain - hier für I=4.

Links

• ...

Literature

• ...

---

Excel-File mit System-Parametern.

Excel-Parameterdatei: Randbedingungen

Parameter des Modells