Grundgleichungen der Finite-Elemente-Methode (FEM)

Ausgangspunkt ist das Prinzip vom Minimum des elastischen Potentials (vgl. Seite 638):

Es ist die Differenz der in dem elastischen System gespeicherten Formänderungsenergie und der “Endwertarbeit” der äußeren Kräfte, die im Allgemeinen einen Integralausdruck liefert. Gesucht sind die Verschiebungsfunktionen v.

Idee des Verfahrens von Ritz (vgl. Seite 640): Für die unbekannten Verschiebungsfunktionen im elastischen Potential wird ein Ansatz mit n “Vergleichsfunktionen” v_i(z) und unbestimmten Koeffizienten a_i in der Form

gewählt (jede Funktion v_i muss die geometrischen Randbedingungen erfüllen). Die a_i werden so bestimmt, dass das elastische Potential den für den gewählten Ansatz noch möglichen minimalen Wert annimmt. Die dafür notwendigen Bedingungen

bilden ein lineares Gleichungssystem mit n Gleichungen für die n Koeffizienten a_i.
 

2  Grundgleichungen der FEM

Während beim “klassischen” Ritz-Verfahren die Ansatzfunktionen für das gesamte Gebiet des zu untersuchenden Bauteils formuliert werden, wird bei der Finite-Elemente-Methode jeweils ein Ritz-Ansatz für die Verschiebungen eines Teilgebiets (finites Element) formuliert. Dabei ergibt sich das Problem, dass die Verschiebungen an den Grenzen der Elemente natürlich mit denen der Nachbarelemente kompatibel sein müssen.

Dies wird mit einem besonders eleganten “Trick” erreicht: Als Freiwerte des Verschiebungsansatzes (beim klassischen Ritz-Verfahren die a_i) werden die Verschiebungen an den Knoten, an denen benachbarte Elemente Kontakt haben, verwendet. Dann ist an diesen Knoten automatisch die Kompatibilität der Verschiebungen benachbarter Elemente gesichert. Man muss allerdings darauf achten, dass die Ansatzfunktionen so beschaffen sind, dass die Verschiebungskompatibilität an den Knoten automatisch die Kompatibilität an der gesamten Grenzfläche benachbarter Elemente sicherstellt.

Betrachtet wird der allgemeine zwei- oder dreidimensionale elastische Körper, für den das Prinzip vom Minimum des elastischen Potentials so formuliert werden kann:

Darin ist e der Vektor der Verzerrungen (Dehnungen und Gleitungen), s der Vektor der Spannungen (Normal- und Schubspannungen), v ist der Vektor der Verschiebungen, p der Vektor der volumenbezogenen Belastung (z. B.: Eigengewicht, Fliehkraft), q der Vektor der Oberflächenlasten (z. B.: Druck).

Zunächst wird nur ein finites Element betrachtet. Es soll n_k Knoten haben, die Knoten haben k_f Knotenfreiheitsgrade (Verschiebungskomponenten, k_f braucht nicht für alle Knoten gleich zu sein), insgesamt möge das Element n_f Freiheitsgrade haben. Für das Element ist nun ein Ansatz für das Verschiebungsfeld in der Form

erforderlich ist (bezogen auf ein elementeigenes Koordinatensystem). In v_e sind alle n_f Verschiebungen an allen Knoten des Elements zusammengefasst (Elementverschiebungsvektor), und weil das Verschiebungsfeld für alle Punkte des Elements (auch für die Knoten) gilt, muss die Matrix der Ansatzfunktionen G so beschaffen sein, dass beim Einsetzen der Koordinaten eines Knotens das Produkt Gv_e den entsprechenden Knotenverschiebungsvektor ergibt.

Gelingt es nicht, eine solche Matrix aufzuschreiben, so kann zunächst ein Verschiebungsansatz

mit allgemeinen Freiwerten a_i gewählt werden (die Anzahl der Freiwerte und damit die Anzahl der Spalten von N müssen mit der Anzahl der Freiheitsgrade n_f des Elements übereinstimmen). Der Elementverschiebungsvektor kann nun durch Einsetzen der Koordinaten der einzelnen Knoten folgendermaßen aufgeschrieben werden:

Die Matrix A enthält ausschließlich Knotenkoordinaten. Mit ihrer Inversen können die allgemeinen Ansatzparameter entsprechend

durch die Knotenverschiebungen ersetzt werden. Damit kann der allgemeine Ansatz in der Form

geschrieben werden, und daraus liest man die Matrix der Ansatzfunktionen

ab, die die eingangs gestellte Forderung nach der Verknüpfung der Elementknotenverschiebungen mit dem Verschiebungsfeld im Inneren des Elements

erfüllt.

Die elastischen Verzerrungen e in der Formel des elastischen Potentials sind mit den Verschiebungen durch eine Beziehung

verknüpft. Hierin ist D ein Matrixdifferenzialoperator. Einsetzen des Verschiebungsansatzes ergibt:

mit

Die Spannungen s in der Formel des elastischen Potentials sind mit elastischen Verzerrungen über ein Stoffgesetz (z. B. das Hookesche Gesetz) verknüpft:

(an die Stelle der Spannungen können auch Schnittgrößen treten - Kirschoffsche Plattentheorie, Bernoullische Biegetheorie, ... -, dann repräsentiert H deren Zusammenhang mit den Krümmungen).

Mit

kann der Anteil eines einzelnen finiten Elements am elastischen Potential so formuliert werden (man beachte, dass sich beim Transponieren eines Produkts die Reihenfolge der Faktoren umkehrt):

Die Knotenverschiebungsvektoren können aus den Integralen herausgezogen weeerden, für die verbleibenden Integrale werden sinnvolle Abkürzungen eingeführt:

mit der sogenannten Elementsteifigkeitsmatrix

und dem Elementbelastungsvektor infolge der Volumen- und Oberflächenlasten

Das gesamte elastische Potential ergibt sich aus der Summierung der Potentiale aller Elemente:

(hier nur angedeutet mit den beiden Elementen e und f).

Zur Erinnerung: Dies ist (nach einer etwas aufwändigen Zwischenrechnung) das Gesamtpotential, aufgeschrieben mit Ritzschen Ansatzfunktionen. Diese sind über die Integralausdrücke in die Elementsteifigkeitsmatrizen und die Elementbelastungsvektoren eingeflossen, während die Ritzschen Ansatzparamter als Knotenverschiebungen (z. B. in v_e und v_f) in separaten Vektoren konzentriert sind.

Die Ansatzparameter (Knotenverschiebungen) werden nun nach der Vorschrift von Ritz durch das Nullsetzen der partiellen Ableitungen des Gesamtpotentials nach den Ansatzparametern bestimmt. Dies kann als Ableitung nach dem Vektor v, in dem alle Knotenverschiebungen aller Knoten des Gesamtsystems enthalten sind, so formuliert werden:

Dies ist einn großes lineares Gleichungssystem

mit der so genannten Systemsteifigkeitsmatrix K und dem Systembelastungsvektor f (eventuell an den Knoten angreifende äußere diskrete Lasten dürfen auf die entsprechenden Positionen in f einfach aufaddiert werden).

Dabei kann in doppelter Hinsicht ausgenutzt werden werden, dass bei der partiellen Ableitung nach einer Knotenverschiebung nur wenige Elemente betroffen sind (genau die, die an diesem Knoten angeschlossen sind). Einerseits führt dies auf eine nur dünn mit von Null verschiedenen Elementen besetzte Koeffizientenmatrix des großen Gleichungssystems, was bei der Lösung mit Vorteil genutzt werden kann. Andererseits bietet es sich an, zunächst nur zu untersuchen, welchen Anteil eine Elementsteifigkeitsbeziehung bei der Ableitung nach den Knotenverschiebungen der eigenen Knoten liefert, um dann danach die Anteile der anderen Elemente zu addieren.

Die Ableitung des Gesamtpotentials nach den in v_e zusammengefassten Knotenverschiebungen des Elements e entsprechend

liefert die so genannte Elementsteifigkeitsbeziehung

die aus der oben genannten Systemsteifigkeitsbeziehung nur jeweils einen Anteil zu einigen Gleichungen liefert. Alle Elementsteifigkeitsmatrizen zusammen ergeben die Systemsteifigkeitsbeziehung, wenn die Anteile, die in K bzw. f auf die gleichen Positionen fallen, addiert werden. Diese “Einspeicherungsstrategie” entspricht genau dem Algorithmus, der im Abschnitt 15.2 am einfachen Beispiel anschaulich hergeleitet wurde.

Es ist sinnvoll, den Weg bis zur Systemsteifigkeitsbeziehung zunächst formal so zu gehen, wie es hier beschrieben wurde. Es muss deshalb vermerkt werden, dass der Vorteil, für alle Elemente die gleichen Ansatzfunktionen zu verwenden, mit einem Verstoß gegen die Regeln des Ritzschen Verfahrens erkauft wurde, denn zulässige Vergleichsfunktionen müssen die geometrischen Randbedingungen erfüllen. Das ist zumindest für die Knoten nicht eingehalten worden, die durch starre Lager gebunden sind.

Aber dieser eindeutige Regelverstoß kann nachträglich geheilt werden, weil die Ansatzparameter ja die in der Systemsteifigkeitsbeziehung noch sichtbaren Knotenverschiebungen sind. Wenn in

nachträglich alle verhinderten Verschiebungen Null gesetzt werden, dann ist das gleichwertig mit dem nach Ritz korrekten Vorgehen, weil damit genau die Ansatzfunktionen, die mit diesen Ansatzparametern verknüpft waren, ihren Einfluss auf das Ergebnis verlieren.

Das Nullsetzen bestimmter Verschiebungen in der Systemsteifigkeitsbeziehung sollte in jedem Fall so realisiert werden, dass die Symmetrie der Matrix K erhalten bleibt. Dies kann z. B. durch das im Abschnitt 15.2 beim Einfügrungsbeispiel demonstrierte “Zeilen-Spalten-Streichen” erreicht werden. Die dabei herausgenommenen Gleichungen enthalten auf der rechten Seite eine Unbekannte in f (Lagerreaktion), die nach Berechnung aller Verschiebungen mit Hilfe dieser Gleichung berechnet werden kann.

Es soll noch angemerkt werden, dass vor dem Einarbeiten der verhinderten Verschiebungen das Gleichungssystem ohnehin nicht gelöst werden kann, weil die Systemsteifigkeitsmatrix singulär ist.

Für den Biegeträger, der nur durch eine Linienlast q belastet ist, kann das Prinzip vom Minimum des elastischen Potentials so formuliert werden:

Verknüpft mit einem Beispiel, findet sich hier eine allgemeinere Formulierung.

Auf den Seiten 640 bis 642 findet man zwei einfache Biegeträger-Aufgaben, die mit dem Verfahren von Ritz gelöst wurden. Etwas anspruchsvollere Probleme findet man über die folgenden Links:

Für Biegeträger kann der Unterschied zwischen “klassischem” Ritz-Verfahren und der Finite-Elemente-Methode besonders anschaulich gezeigt werden:

Beim “klassischen” Ritz-Verfahren werden Ansatzfunktionen verwendet, die für die gesamte Länge des Trägers gelten, bei der Finite-Elemente-Methode wird zunächst das Element e betrachtet. In einem elementeigenen Koordinatensystem (Koordinate z_e) werden die Ansatzfunktionen für dieses Element aufgeschrieben. Das Element für den geraden Biegeträger soll 2 Knoten mit je zwei Freiheitsgraden haben (Vertikalverschiebung und Biegewinkel). Hier wird also ein Ansatz in der Form

mit vier Ansatzfunktionen benötigt, die beim Einsetzen der Knotenkoordinaten z_e = 0 (Knoten i) bzw. z_e = l_e (Knoten j) die zugehörige Knotenverformung liefern (das bedeutet z. B., dass g₁(0) = 1 und g₂(0) = g₃(0) = g₄(0) = 0 liefern müssen). Man beachte, dass dies gleichzeitig die Ansatzfunktionen für den Biegewinkel sind, denn es gilt natürlich

Hier soll demonstriert werden, wie man mit einem allgemeinen Ansatz zu den Ansatzfunktionen g kommt (Matrizen im links beschriebenen allgemeinen Fall, die für den Biegeträger zu Vektoren degenerieren, werden hier mit Kleinbuchstaben bezeichnet). In

bzw.

müssen entsprechend den 4 Freiheitsgraden des Elements sowohl vier Ansatzfunktionen als auch vier Freiwerte enthalten sein. Nun wird der Elementverschiebungsvektor durch Einsetzen der Knotenkoordinaten in den allgemeinen Ansatz erzeugt:

Für diesen einfachen Fall kann die Matrix A (mit erträglicher Mühe) “von Hand” invertiert werden, und man kann aus

den Vektor der Ansatzfunktionen g berechnen:

mit

Die Analogie zu den nebenstehend beschriebenen allgemeinen Gleichungen soll möglichst deutlich bleiben. Dort gibt es im elastischen Potential den Anteil

Das erste Summand im elastischen Potential für den Biegeträger wird dieser Schreibweise angepasst:

Der Matrix H entspricht hier die Biegesteifigkeit EI, und dem Vektor e entspricht hier die zweite Ableitung der Verschiebung. Der Differentialoperator D wird also besonders einfach, es ist die 2. Ableitung nach der Koordinate z_e:

Die Matrix D_G degeneriert in diesem einfachen Fall zu einem Vektor d_G, der sich nach der Formel

berechnen lässt:

Nun wird der Anteil eines einzelnen Elements am elastischen Potential durch Einsetzen von

formuliert:

Die Knotenverschiebungsvektoren können aus den Integralen herausgezogen werden:

mit der Elementsteifigkeitsmatrix K_e und dem Elementbelastungsvektor f_e:

Das gesamte elastische Potential ergibt sich aus der Summe aller Elementanteile (hier nur mit den beiden Nachbarelementen e und f angedeutet):

Die Ansatzparameter (Knotenverschiebungen) werden nun nach der Vorschrift von Ritz durch das Nullsetzen der partiellen Ableitungen des Gesamtpotentials nach allen Ansatzparametern (Knotenverschiebungen) v₁, φ₁, v₂, φ₂, ... , v_i, φ_i,  v_j, φ_j, ...  bestimmt, was als

formuliert werden kann (v ist der Vektor mit allen Knotenverschiebungen). Die Ableitungen liefern ein lineares Gleichungssystem

mit der Systemsteifigkeitsmatrix K und dem Systembelastungsvektor f (eventuell an den Knoten angreifende äußere diskrete Lasten dürfen auf die entsprechenden Positionen in f einfach aufaddiert werden).

Für den geraden Biegeträger wird besonders deutlich, dass bei den Ableitungen immer nur wenige Elemente betroffen sind. Für die partiellen Ableitungen nach v_j und φ_j z. B. sind es nur die in der Formel des Gesamtpotentials angedeuteten vier Anteile der beiden Elemente e und f, weil nur diese beiden Elemente Kontakt zum Knoten j haben und demzufolge v_j und φ_j nur in den Elementverschiebungsvektoren v_e und v_f vorkommen.

Es ist also sinnvoll, zunächst nur die Ableitung des Gesamtpotentials nach einem Elementverschiebungsvektor v_e auszuführen, um so die Anteile zu bestimmen, den dieses Element zu den Gleichungen für die Knoten i und j liefert und später dann die Anteile der Nachbarelemente hinzuzufügen Aus

erhält man die Elementsteifigkeitsbeziehung

Nach den oben angegebenen Integralformeln können K_e und f_e für diesen einfachen Fall in geschlossener Form angegeben werden (Einsetzen der beiden Vektoren d_G und g und Ausmultiplizieren):

mit   

und   

Für den Sonderfall “Konstante Biegesteifigkeit EI” können die Integrale der k_mn-Formel sogar allgemein gelöst werden, für k₁₂ erhält man z. B.:

Für vorgegebene Funktionen der Linienlast q kann man auch die Integrale der f_m-Formel lösen, z. B. für linear veränderliches q mit der Intensität q_i am Knoten i und q_j am Knoten j. Gemeinsam mit konstanter Biegesteifigkeit ergibt sich dann die Elementsteifigkeitsbeziehung, die im Abschnitt 18.2 als Formel (18.11) auf anderem Wege hergeleitet wurde:

Man beachte, dass im Elementbelastungsvektor auf den Positionen 1 und 3 Kräfte und auf den Positionen 2 und 4 Momente stehen.

3  Beispiel: Aufbau der Systemsteifigkeitsbeziehung, Einbau der geometrischen Randbedingungen

Am Beispiel eines Biegeträgers sollen der Aufbau der Elementsteifigkeitsbeziehungen, der Zusammenbau zur Systemsteifigkeitsbeziehung, der Einbau der geometrischen Randbedingungen und die Berechnung der Verschiebungen demonstriert werden.

Aufgabe: Für den nebenstehend skizzierten Träger mit stückweise konstanter Biegesteifigkeit soll die Biegeverformung berechnet werden.

Gegeben:

l = 1 m ; F = 1 kN ;  q₁ = 2 kN/m ; q₂ = 1 kN/m ;
EI₁ = EI₃ = 1 kNm² ;   EI₂ = 2 kNm² .

Der Träger wird in 3 finite Elemente unterteilt, dabei entstehen 4 Knoten. Die nebenstehende Skizze zeigt die gewählte Nummerierung.

Das Matlab-Interface zum Finite-Elemente-Baukasten Femset gestattet Einblicke in die Zwischenstufen des Algorithmus. Deshalb wird mit dem nebenstehend zu sehenden Matlab-Script zunächst in den Zeilen 5 bis 9 das Berechnungsmodell definiert, danach werden die drei Elementsteifigkeitsbeziehungenen (Zeilen 12 bis 14), die Systemsteifigkeitsbeziehung ohne (Zeile 16) bzw. mit Berücksichtigung der Randbedingungen (Zeile 19) berechnet und in das Command Window ausgegeben. Abschließend werden die Verschiebungen berechnet (Zeile 24).

Das Berechnungsmodell: In xy stehen die Koordinaten der 4 Knoten, bezogen auf ein (beliebiges) hier im Knoten 1 liegendes Koordinatensystem, km (“Koinzidenzmatrix”) enthält die Zuordnung der Knoten zu den Elementen (zu Element 1 gehören Knoten 1 und 2 usw.). In der Elementparametermatrix ep stehen in jeder Zeile 3 Werte für ein Element: Biegesteifigkeit EI, Linienlastintensität am ersten Knoten und Linienlastintensität am zweiten Knoten. Die Matrix der Randbedingungen kr signalisiert für die jeweils zwei Freiheitsgrade (Vertikalverschiebung, Biegewinkel) der Knoten mit einer 1 eine “verhinderte Verschiebung”, mit einer 0 die freie Verschiebungsmöglichkeit. In der Matrix der diskreten Lasten können analog dazu für jeden Knoten eine äußere Einzellast und ein äußeres Moment definiert werden.

Die Zahlenwerte der Aufgabenstellung wurden bewusst so gewählt, dass man in den Elementsteifigkeitsmatrizen die Werte aus der oben angegebenen Formel wiedererkennt (es wurden auch im Berechnungsmodell konsequent die Einheiten kN und m für alle Werte verwendet, so dass sich die Verschiebungen auch mit der Dimension m ergeben).

Nebenstehend sieht man nur den ersten Teil (Ausgabe des Ergebnisses der Berechnung in Zeile 12) des Command Windows (Elementsteifigkeitsbeziehungen des Elements 1). Alle übrigen Ausgaben sind im nachfolgenden Bild so zusammengestellt worden, dass man die Zusammenhänge zwischen Elementsteifigkeitsbeziehungen und Systemsteifigkeitsbeziehung erkennen kann.

Man sieht, dass die Elementsteifigkeitebeziehungen direkt auf die entsprechenden Positionen in die Systemsteifigkeitsbeziehung eingespeichert werden. Au f den Positionen, auf die Anteile aus verschiedenen Elementsteifigkeitsbeziehungen gelangen (hier sind es die Bereiche, wo sich die Anteile des Elements 2 mit einem der beiden anderen Elemente überlappen), werden diese addiert. Anschließend wurde im Systembelastungsvektor die äußere Einzellast am Knoten 2 ergänzt.

In der Systemsteifigkeitsbeziehung (K und f) sind zunächst die geometrischen Randbedingungen (verhinderte Verschiebungen) unberücksichtigt. Am Knoten 1 sind beide Verschiebungen (Vertikalverschiebung, Biegewinkel), am Knoten 4 nur die Vertikalverschiebung verhindert. Das bedeutet, dass im Gleichungssystem

im Vektor v auf den Positionen 1, 2 und 7 jeweils eine 0 steht, und im Vektor f müssen auf den entsprechenden Positionen unbekannte Kräfte (Lagerreaktionen) eingesetzt werden. Bei der Multiplikation Kv sind also (infolge der Nullelemente in v) die Spalten 1, 2 und 7 der Matrix K bedeutungslos. Diese Spalten und die Nullelemente in v könnten also gestrichen werden. Wenn gleichzeitig auch die Zeilen 1, 2 und 7 (in K und f) gestrichen werden, verschwinden diese 3 Gleichungen (und auch die drei Unbekannten in f) aus dem Gleichungssystem.

In der Systemsteifigkeitsbeziehung sind diese Zeilen und Spalten gelb unterlegt. Wenn sie gestrichen werden (Berücksichtigung verhinderter Verschiebungen durch “Zeil-Spalten-Streichen”), verbleibt ein Gleichsungssystem mit nur noch 5 Gleichungen (und 5 unbekannten Verschiebungen), wobei die Symmetrie der Matrix K erhalten bleibt.

In Femset wird eine alternative Strategie zur Berücksichtigung der Randbedingungen realisiert: Die Zeilen und Spalten werden nicht gestrichen, sondern durch die drei Gleichungen v₁ = 0 , φ₁ = 0 und v₄ = 0 ersetzt. Die auf diese Weise modifizierte Systemsteifigkeitsbeziehung ist im nachfolgenden Bild unten zu sehen. Dem geringfügigen Nachteil dieser Strategie, dass sich die Anzahl der Gleichungen nicht reduziert, stehen als Vorteile die Einsparung eines aufwändigen Umspeicherungsprozesses und ein Ergebnisvektor v nach der Lösung des Gleichungssystems gegenüber, der alle Verschiebungen des Systems (auch die Nullverschiebungen) enthält.

Es fällt auf, dass die Gleichungen, die die verhinderten Verschiebungen realisieren, statt der zu erwartenden 1 eine 16 auf der Hauptdiagonale haben. Natürlich ist 16 v₁ = 0 gleichwertig mit v₁ = 0. Weil die Elemente der Steifigkeitsmatrizen außerordentlich groß werden können und bei stark unterschiedlichen Größenordnungen der Elemente der Koeffizientenmatrix die Gefahr numerischer Instabilität besteht, setzt Femset einen Wert in der Größenordnung der übrigen Matrixelemente auf diese Positionen.

Schließlich liefert das Matlab-Script (Zeile 24) noch die sich aus der Lösung des Gleichunggsystem ergebenden Knotenverschiebungen ab (siehe nebenstehenden Ausschnitt aus dem Command Window, in der linken Spalte stehen die Vertikalverschiebungen, in der rechten Spalte die Biegewinkel an den vier Knoten).

Das hier beschriebene Matlab-Script ist als Lösung der Aufgabe 18-14 verfügbar und steht dort gemeinsam mit allen verwendeten DLLs und Matlab-Functions zum Download bereit.

4  Dynamische Probleme (Eigenschwingungen elastischer Systeme)

Bei dynamischen Problemen bietet sich das Hamiltonsche Prinzip

als Ausgangspunkt an: Zwischen zwei Zeitpunkten t₁ und t₂ verläuft die Bewegung eines Massensystems so, dass das Zeitintegral über die (für konservative Systeme geltende) Lagrangesche Funktion

einen stationären Wert annimmt. Π ist das elastische Potential, T die kinetische Energie des Systems, die für ein einzelnes finites Element in der Form

aufgeschrieben werden kann mit

• - Dichte des Materials (ist in der Regel konstant),
• - Geschwindigkeitsfeld der Massenpunkte.

Für das Geschwindigkeitsfeld wird die gleiche Strategie gewählt, mit der im elastostatischen Fall die Verschiebungen durch die Knotenverschiebungen ausgedrückt wurden:

stellt den Zusammenhang der Knotengeschwindigkeiten mit dem Geschwindigkeitsfeld im Inneren des Elements mit den gleichen Ansatzfunktionen her, die auch für den elastostatischen Fall verwendet (und wie dort erzeugt) werden. Damit ergibt sich der Integralausdruck für die kinetische Energie in einem Element zu

(die Knotengeschwindigkeitsvektoren können aus dem Integral herausgezogen werden) mit der so genannten Elementmassenmatrix

Damit kann die Lagrangesche Funktion für ein finites Element aus kinetischer Energie und elastischem Potential zusammengesetzt werden:

Zum Hamiltonschen Prinzip gehören die “Lagrangeschen Gleichungen 2. Art” (vgl. Abschnitt 33.4.3 auf Seite 633)

die mit der Lagrangeschen Funktion für ein Element auf folgendes Differenzialgleichungssystem führt:

Hierin sind K_e und f_e die Elementsteifigkeitsmatrix und der Elementbelastungsvektor wie für das elastostatische Problem, M_e ist die Elementmassenmatrix, für die die gerade entwickelte Formel gilt.

Auch hier sind die Überlegungen zunächst für ein Element angestellt worden. Das Zusammensetzen der Elementmassenmatrizen zur Systemmassenmatrix erfolgt nach dem gleichen “Einspeicherungs -Algorithmus” wie für die Steifigkeitsmatrizen. An den Knoten zusätzlich vorhandene diskrete Massen dürfen auf die entsprechenden Hauptdiagonalelemente der Systemmassenmatrix addiert werden.

Im Gegensatz zum elastostatischen problem entsteht bei dynamischen Problemen ein Differenzialgleichungssystem

Wenn z. B. die freien Schwingungen (f = o) berechnet werden sollen, führt der aus der Schwingungslehre bekannte Ansatz (vgl. Kapitel 32 auf Seite 603)

auf das allgemeine symmetrische Matrizeneigenwertproblem

das die Eigenkreisfrequenzen der Schwingung ω und die zugehörigen Eigenschwingungsformen liefert.

Für den geraden Biegeträger wird angenommen, dass sich (wie in der Skizze angedeutet) die Masseteilchen nur vertikal bewegen.

Zunächst wird auch hier nur das Element e mit den Knoten i und j betrachtet, das zusätzlich zu den Eigenschaften, die es für das elastostatische Problem hatte, noch durch eine Massebelegung ρA charakterisiert wird, die wie die Biegesteifigkeit veränderlich sein darf.

Für das Verschiebungsfeld wird der gleiche Ansatz wie im elastostatischen Fall verwendet, für das Geschwindigkeitsfeld (Geschwindigkeit der Masseteilchen in vertikaler Richtung) wird der folgende Ansatz

 gemacht, wobei die gleichen Funktionen g_i(z_e) verwendet werden wie für das Verschiebungsfeld (also die oben für das elastostatische Problem ermitteltenn Funktionen g₁ bis g₄).

Damit ergibt sich die gleiche Elementsteifigkeitsmatrix wie für das elastostatische Problem, und die Elementmassenmatrix errechnet sich nach der Formel (vgl. nebenstehend beschriebenen allgemeinen Fall)

mit

und

Sie ergibt sich in der gleichen Form wie die Elementsteifigkeitsmatrix:

mit   

Für den Sonderfall “Konstante Massenbelegung ρA” können die Integrale der m_mn-Formel auch hier allgemein gelöst werden, für m₂₄ erhält man z. B.:

Und das ist die komplette Elementmassenmatrix für konstante Massenbelegung des Biegeträgers:

 Der Algorithmus für den Aufbau der Systemmatrizen ist in diesem Fall für Systemsteifigkeitsmatrix und Systemmassenmatrix gleich (Einspeicherungsstrategie mit Addition der Elemente, die auf gleiche Positionen kommen).

Beide Matrizen sind symmetrisch und haben in der Regel eine ausgeprägte Bandstruktur. Dies kann gegebenenfalls bei der Lösung des Matrizeneigenwertproblems für die Berechnung von Eigenschwingungen mit Vorteil genutzt werden. Auch hier müssen natürlich vor der Lösung des Eigenwertproblems die geometrischen Randbedingungen (verhinderte Verschiebungen) eingebaut werden. Dafür wird die Strategie des “Zeilen-Spalten-Streichens” dringend empfohlen.

_{5  Beispiel: Eigenschwingungen eines geraden Biegeträgers}

Am Beispiel eines einfachen Biegeträgers sollen der Aufbau der Elementmatrizen, der Zusammenbau zu den Systemmatrizen, der Einbau der geometrischen Randbedingungen und die Lösung des Eigenwertproblems (Berechnung der Eigenfrequenzen und der Eigenschschwingungsformen) demonstriert werden.

Aufgabe: Für den nebenstehend skizzierten Träger mit konstanter Biegesteifigkeit und konstanter Massebelegung (und einer zusätzlichen Einzelmasse am freien Ende)sollen die beiden kleinsten Eigenfrequenzen der Biegeschwingung und die zugehörigen Eigenschwingungsformen berechnet werden.

Gegeben: EI = 3000 Nm² ; ρA = 3 kg/m ;  m = 2 kg ;  l = 1 m .

Der Träger wird in nur 2 finite Elemente unterteilt, dabei entstehen 3 Knoten. Die nebenstehende Skizze zeigt die gewählte Nummerierung.

Auch für die Eigenschwingungsberechnungen gestattet das Matlab-Interface zum Finite -Elemente-Baukasten Femset Einblicke in die Zwischenstufen des Algorithmus. Deshalb wird mit dem nachstehend zu sehenden Matlab-Script zunächst in den Zeilen 5 bis 9 das Berechnungsmodell definiert, danach werden die jeweils zwei  Elementmatrizen der beiden Elemente (Zeilen 12 und 13), die Systemmatrizen ohne (Zeilen 15 bis 17) bzw. mit Berücksichtigung der Randbedingungen (Zeilen 19 bis 21) berechnet und in das Command Window ausgegeben. Abschließend wird das Matrizeneigenwertproblem gelöst (Zeile 24).

Das Berechnungsmodell: In xy stehen die Koordinaten der 3 Knoten, bezogen auf ein (beliebiges) hier im Knoten 1 liegendes Koordinatensystem, km (“Koinzidenzmatrix”) enthält die Zuordnung der Knoten zu den Elementen (zu Element 1 gehören Knoten 1 und 2, zu Element 2 die Knoten 2 und 3). In der Elementparametermatrix ep stehen in jeder Zeile 2 Werte für ein Element: Biegesteifigkeit EI und Massebelegung ρA. Die Matrix der Randbedingungen kr signalisiert für die jeweils zwei Freiheitsgrade der Knoten (Vertikalverschiebung, Biegewinkel) mit einer 1 eine “verhinderte Verschiebung”, mit einer 0 die freie Verschiebungsmöglichkeit. In der Matrix der diskreten Massen mn können analog dazu für jeden Knoten eine Einzelmasse und ein Massenträgheitsmoment definiert werden.

Weil die Ergebnisse (Eigenkreisfrequenzen) in der Dimension die Zeit enthalten, die bei den Größen des Berechnungsmodells gar nicht vorkommt, gilt die dringende (und hier natürlich eingehaltene) Empfehlung: Für das Berechnungsmodell sollten konsequent die Einheiten kg, N und m für alle Werte verwendet werden. Dann erhält man die Eigenkreisfrequenzen mit der Dimension s ^-1.

Alle Ergebnisse werden in das Command Window ausgegeben. Sie wurden im folgenden Bild mit Erläuterungen so zusammengestellt, dass man die Zusammenhaänge zwischen den einzelnen Matrizen erkennen kann.

Weil beide Elemente (mit gleichen Abmessungen und Materialeigenschaften) identisch sind, sind auch ihre Elementsteifigkeitsmatrizen bzw. Elementmassenmatrizen gleich. Man sieht, dass die Elementmatrizen direkt auf die entsprechenden Positionen in die Systemmatrizen eingespeichert werden. Au f den Positionen, auf die Anteile aus verschiedenen Elementmatrizen gelangen (hier sind es die Bereiche, die zum Knoten 2 gehören), werden diese addiert. Anschließend wurde in der Systemmassenmatrix auf dem ersten der beiden zum Knoten 3 gehörenden Hauptdiagonalelemente die Zusatzmasse m = 2 kg ergänzt.

In den Systemmatrizen (K und M) sind zunächst die geometrischen Randbedingungen (verhinderte Verschiebungen) unberücksichtigt. Am Knoten 1 sind beide Verschiebungen (Vertikalverschiebung, Biegewinkel), am Knoten 2 nur die Vertikalverschiebung verhindert. Weil für Eigenwertprobleme das Berücksichtigen der verhinderten Verschiebungen durch “Zeilen-Spalten-Streichen” zu empfehlen ist, ist diese Strategie auch in Femset realisiert. Das bedeutet, dass in beiden Systemmatrizen K und M die ersten 3 Zeilen und die ersten 3 Spalten zu streichen sind. Es verbleibt das allgemeine symmetrische Matrizeneigenwertproblem mit zwei (3*3)-Matrizen.

Nach der Lösung des Matrizeneigenwertproblem (Zeile 24 im Matlab-Script) werden die Eigenkreisfrequenzen in Zeile 25 auf Frequenzen umgerechnet und in das Command Window ausgegeben. In Zeile 26 werden die beiden Eigenschwingungsformen in ein Graphik-Fenster gezeichnet (nebenstehendes Bild).

Die Ergebnisse im Command Window

können mit den exakten Werten verglichen werden, weil dieses einfache Beispiel einer Lösung der Schwingungsdiffernzialgleichung für den geraden massebehafteten Träger zugänglich ist:

f_1,exakt = 20,78 s ^-1 ;    f_2,exakt = 242,13 s ^-1 .

Die Schlussfolgerungen sind verallgemeinerungsfähig: Während die erste Eigenfrequenz schon mit zwei Elementen hervorragend angenähert wird, zeigt schon die zweite Eigenfrequenz erhebliche Abweichungen vom exakten Wert. Die Finite -Elemente-Methode ist ein Näherungsverfahren, und die Ergebnisse sind immer nur so gut, wie die Ansatzfunktionen in der Lage sind, die exakten Verformungen anzunähern. Bei einer Rechnung mit mehr als zwei Elementen werden auch die Ergebnisse für die höheren Frequenzen sehr schnell besser. Schon bei einer Einteilung der Trägerlänge in 8 Elemente erhält man auch für den 2. Eigenwert mit

 f_2,exakt = 242,21 s ^-1

einen ausgezeichneten  Näherungswert.