Methode der Finiten Elemente - Strukturmechanik · daher hat FEM ohne adequate Rechneruntersttzung...
Transcript of Methode der Finiten Elemente - Strukturmechanik · daher hat FEM ohne adequate Rechneruntersttzung...
Methode der Finiten Elemente
Skript zur Vorlesung
fur den Studiengang Produktionstechnik
Dr.-Ing. M. Mehrafza
WS 2006/07
i
Inhaltsverzeichnis
Literaturverzeichnis v
1 Grundkonzepte 3
1.1 Losung der Randwertprobleme . . . . . . . . . . . . . . . . 4
1.2 Integralformulierungen fur die numerische Losung . . . . . 7
1.2.1 Variationsmethode . . . . . . . . . . . . . . . . . . 8
1.2.2 Punktweise Kollokation . . . . . . . . . . . . . . . . 8
1.2.3 Bereichsweise Kollokation . . . . . . . . . . . . . . 9
1.2.4 Galerkin-Verfahren . . . . . . . . . . . . . . . . . . 10
1.2.5 Verfahren der kleinsten Quadrate . . . . . . . . . . 10
1.3 Potentialenrgie-Formulierung . . . . . . . . . . . . . . . . . 11
1.4 Die Methode der Finiten Elemente . . . . . . . . . . . . . 12
2 Eindimensionales lineares Element 15
2.1 Aufteilung des Gebiets in Elemente . . . . . . . . . . . . . 15
2.2 Das lineare Element . . . . . . . . . . . . . . . . . . . . . . 16
2.3 Eine kontinuierliche stuckweise glatte Funktion . . . . . . . 18
2.4 Ein Beispiel zur FE-Formulierung . . . . . . . . . . . . . . 19
2.4.1 Gewichtsfunktionen . . . . . . . . . . . . . . . . . . 20
2.4.2 Das gewichtete Restintegral . . . . . . . . . . . . . 21
2.4.3 Auswertung der Integrale . . . . . . . . . . . . . . . 22
2.4.4 Analyse eines Biegebalkens . . . . . . . . . . . . . . 25
iii
iv INHALTSVERZEICHNIS
2.4.5 Matrizendarstellung . . . . . . . . . . . . . . . . . . 27
2.5 Elementmatrizen: Galerkin-Formulierung . . . . . . . . . . 28
2.5.1 Elementmatrizen . . . . . . . . . . . . . . . . . . . 28
2.5.2 Direkte Steifigkeitsmethode . . . . . . . . . . . . . 30
2.5.3 Analyse eines Biegebalkens . . . . . . . . . . . . . . 31
2.5.4 Eigenschaften der globalen Steifigkeitsmatrix . . . . 34
3 Zweidimensionale Elemente 35
3.1 Zweidimensionale Netze . . . . . . . . . . . . . . . . . . . 35
3.2 Das lineare Dreieckselement . . . . . . . . . . . . . . . . . 36
3.3 Das bilineare Rechteckselement . . . . . . . . . . . . . . . 39
3.4 Eine kontinuierliche stuckweise glatte Funktion . . . . . . . 41
4 Koordinatensysteme 45
4.1 Lokale Koordinatensysteme . . . . . . . . . . . . . . . . . 45
4.2 Naturliche Koordinatensysteme . . . . . . . . . . . . . . . 48
4.3 Das Rechteckselement . . . . . . . . . . . . . . . . . . . . 50
4.4 Das Dreieckselement: Flachenkoordinaten . . . . . . . . . . 51
4.5 Kontinuitat . . . . . . . . . . . . . . . . . . . . . . . . . . 56
5 Zweidimensionale Feldprobleme 59
5.1 Problembeschreibende Differentialgleichungen . . . . . . . 60
5.2 Integralgleichungen fur die Elementmatrizen . . . . . . . . 62
5.3 Elementmatrizen: Dreieckselement . . . . . . . . . . . . . . 65
5.4 Elementmatrizen: Rechteckselement . . . . . . . . . . . . . 67
5.5 Torsion eines nichtkreiformigen Querschnitts . . . . . . . . 70
5.5.1 Allgemeine Theorie . . . . . . . . . . . . . . . . . . 70
5.5.2 Verdrehung des quadratischen Torsionsstabs . . . . 72
5.5.3 Komponenten der Schubspannung . . . . . . . . . . 76
5.5.4 Berechnung des Torsionsmoments . . . . . . . . . . 78
INHALTSVERZEICHNIS v
6 2D-Feldprobleme mit gemischten Randbedingungen 81
6.1 Ableitungsenthaltene Randbedingungen . . . . . . . . . . . 82
6.1.1 Auswertung von Elementintegralen . . . . . . . . . 85
6.2 Punktquellen und Punktsenken . . . . . . . . . . . . . . . 88
6.3 Drehfreie Stromung idealer Flussigkeit . . . . . . . . . . . 91
6.3.1 Stromung einer idealen Flussigkeit . . . . . . . . . 91
6.3.2 Grundwasserstromung . . . . . . . . . . . . . . . . 95
6.4 Warmeubertragung durch Leitung und Strahlung . . . . . 96
6.4.1 Eindimensionale Schaufel . . . . . . . . . . . . . . . 97
6.4.2 Verbundwand . . . . . . . . . . . . . . . . . . . . . 99
6.4.3 Zweidimensionale Schaufel . . . . . . . . . . . . . . 100
6.4.4 Lange zweidimensionale Korper . . . . . . . . . . . 102
6.5 Akustische Vibrationen . . . . . . . . . . . . . . . . . . . . 106
6.5.1 Eindimensionale Vibrationen . . . . . . . . . . . . . 106
6.5.2 Zweidimensionale Vibrationen . . . . . . . . . . . . 109
Vorwort
Methode der Finiten Elemente ist eine sehr verbreitete numerische Pro-zedur fur die Losung der Differentialgleichungen in der Ingenieurwissen-
schaften und Physik und stellt die rechnerische Basis von vielen CAD-Systemen (Computational Aided Design Systems) dar. Das Lernen der
Grundlagen dieser Methode ist langst fur die jenigen, die sich mit derLosung von Problemen in Struktur- und Kontinuumsmechanik, Warme-
leitung und Hydromechanik oder andere Feldprobleme beschaftigen, un-abdingbar geworden. In dieser Vorlesung geht es um die Grundlagen derFEM und ihre Anwendung in der Analyse von ebenen Strukturen und
zweidimensionalen Kontinuumsproblemen der Elastizitat.
1
Kapitel 1
Grundkonzepte
Die Methode der Finiten Elemente (ab jetzt nur durch FEM bezeichnet)kann man in zwei Primare Methoden aufteilen. Die erste Methode macht
von diskreten Elementen Gebrauch, um die Verschiebungen und Element-krafte von Struktursystemen zu berechnen. Die zweite Methode setzt
Kontinuumselemente ein, um Naherungslosungen fur Festkorpermechanik,Warmeleitung und andere Feldprobleme zu bekommen. Die erste Formu-
lierung wird auch Matritzen-Analyse von Struktrukturen genannt und lie-fert Ergebnisse, die mit den Ergebnissen aus der klassischen Strukturme-chanik identisch sind. Das zweite Verfahren ist die echte FE-Methode, und
liefert approximative Werte fur gewnschte Großen an bestimmten Punk-ten innerhalb des Systems. Ein allgemeines FE-Programm ist in der Lage
beide Aufgabeformen zu losen und der Begriff FEM wird oft fur beidenFormulierungen (mit diskreten und mit Kontinuumselemeneten) verwen-
det. Die FE-Methode kombiniert einigen mathematischen Konzepte umein System von linearen oder nichtlinearen Gleichungen aufzustellen. Die
Zahl der Gleichungen ist meist sehr groß- kann mehrere ZehntausendeGleichungen erreichen - und benotigt entsprechend hohe Rechnerleistung,daher hat FEM ohne adequate Rechneruntersttzung kaum praktischen
Wert. Die Vorteile der FEM liegen auf der Hand: Sie kann fur Objek-te eingesetzt werden, die aus mehreren Materialien bestehen, gemischte
Randbedingungen besitzen, keine regulare Form haben, materielle Nicht-linearitat aufweisen, zeitunabhangig oder zeitabhangig sind. Daher stellt
FEM die rechnerische Basis vieler CAD-Programme dar und macht es furpraktische Ingenieure notwendig zu wissen, wie diese Methode funktio-
niert.
3
4 KAPITEL 1. GRUNDKONZEPTE
1.1 Losung der Randwertprobleme
Wenn ein Problem durch eine Differentiagleichung (Dgl) beschrieben wird,ist die analytische Losung der Dgl der beste Weg zur Problemlosung, aber
es gibt zahlreiche Situationen, bei denen man eine analytische Losungschwer finden kann. Die betrachtet Problemregion kann so irregular sein,
dass eine mathematische Beschreibung der Grenzen unmoglich ist. DasProblemgebiet kann aus Kombination verschiedener Materialien beste-
hen, deren Bereichen mathematisch schwer beschreibbar sind. Problememit nicht isotropen Materialien haben in der Regel Dgl mit nicht linea-
ren Termen, die analytische schwer zu losen sind. Wenn eine analytischeLosung schwer zu finden ist,kann eine numerische Methode uns eine ap-proximative Losung anbieten. Alle numerische Losungsmethoden erzeugen
fur einen bestimmten Satz von Parametern Werte an bestimmten Punktendes Aufgabengebiets. Bei jeder Anderung der Parameter wird der kom-
plette Losungsprozedur wiederholt. Die berechneten Werte geben wichtigeInformationen ber den physikalischen Prozess, auch wenn sie nur zu dis-
kreten Punkten gehoren. Es gibt mehrere methoden fur die numerischeLosung einer Dgl. Diese Methoden kann man in drei Gruppen aufteilen:(1) Finite Differenzen Methode, (2) Variationsmethode, und (3) Methode
des gewichteten Restes. Diese Methoden werden hier kurz angesprochen.
Finite Differenzen Methode
Diese Methode approximiert Ableitungen der Dgl mit Hilfe der Diffe-
renzengleichungen. Sie ist fur die Losung der Warmeleitungsproblemeund Fluidmechanik besonders geeignet und funktioniert am besten, wenn
die Grenzen des Problemgebiets zu den Korrdinatenachsen paralell sind.FurProbleme mit gekrmmten oder irregularen Randern ist diese Metho-de sehr lastig und darber hinaus ist es sehr schwer diese Methode in ein
allgemeines Computer-Programm umzusetzen.
Variationsmethode
Die Herangehensweise der Variationsmethode beinhaltet das Integral einer
Funktion, die eine Zahl erzeugt. Jede neue Funktion erzeugt eine neueZahl. Eine Funktion, die die kleinste Zahl erzeugt, besitzt die zusatzliche
Eigenschaft fur die Erfllung einer bestimmten Dgl. Um dieses Konzept zu
1.1. LOSUNG DER RANDWERTPROBLEME 5
erlautern, betrachten wir das folgende Integral
Π =
∫ H
0
[
D
2
(
dy
dx
)2
−Qy
]
dx . (1.1)
Der numerische Wert von Π kann fur eine gegebene Gleichung y = f(x)gerechnet werden. Die Variationsrechnung zeigt, dass eine bestimmte Glei-
chung y = f(x), die den kleinsten numerischen Wert fur Π liefert, dieLosung der Dgl
Dd2y
dx2+Q = 0 (1.2)
mit den folgenden Randbedingungen ist: y(0) = y0 und y(H) = yH .
Diesen Prozess kann man auf ungekehrte Weise nutzen. Ist eine Dgl ge-geben, dann kann eine Naherungslosung erreicht werden, indem man un-terschiedliche Versuchsfunktionen in das Integral einsetzt. In diesem Fall
bezeichnet man den Integralausdruck als Approimationsfunktional odereinfach Funktional. Die Versuchsfunktion, die den kleinsten Wert (Mini-
mum) fur Π liefert, ist die Naherungslosung der Dgl 1.2. Die Variations-methode ist die Basis vieler FE-Formulierungen, aber hat einen großen
Nachteil: Sie ist fur Differentialgleichungen (Dgln), die einen Term miterster Ableitung haben, nicht einsetztbar.
Methoden des gewichteten Restes; Ritz’sches Verfahren
Diese Methoden arbeiten auch mit einem Integral. Hier wird eine Nahe-rungslosung in die Dgl eingesetzt. Da diese Losung der Dgl nicht ganz
erfllt, bleibt ein Fehler oder ein Rest brig. Nehmen wir an, dass y = h(x)eine Naherungslosung fur 1.2 darstellt. Einsetzen ergibt
Dd2h(x)
dx2+Q = R(x) 6= 0 . (1.3)
Die Mehtoden des gewichteten Restes verlangen, dass die Gleichung
∫ H
0
Wi(x)R(x)dx = 0 (1.4)
erfllt wird. Der Rest R(x) ist mit einer Gewichtsfunktion Wi(x) multipli-ziert und das Integral des Produktes muss gleich Null sein. Die Anzahl der
Gewichtsfunktionen entspricht der Anzahl der unbekannten Koeffitienten
6 KAPITEL 1. GRUNDKONZEPTE
in der Naherungslosung. FurGewichtsfunktionen stehen viele Varianten
zur Wahl. Einige der bekanntesten sind:
Punktweise Kollokation : Impulsfunktionen Wi(x) = δ(x−Xi) werden
als Gewichtsfunktionen gewahlt. Das entspricht der Aufforderung, dassder Rest an bestimmten Punkten verschwinden muss. Die Anzahl der
gewahlten Punkte gleicht der Anzahl der unbekannten Koeffizienten inder Naherungslosung. Damit bekommt man die gewnschte Anzahl derGleichungen, die man fur die Bestimmung der unbekannten Koeffizienten
braucht.
Kollokation im Teilbereich : Der gesamte Losungsbereich wird in n
Teilbereiche untergliedert. Die Gewichtsfunktion wird ber jeden Teilbe-reich gleich eins gesetzt, Wi(x) = 1. Das entspricht der Aufforderung,
dass das Integral des Restes ber jeden Teilbereich verschwindet. So erhohtman wiederum n Gleichungen fur die Bestimmung von n unbekannten
Koeffizienten der Naherungslosung.
Galerkin-Verfahren : Bei diesem Verfahren werden fur die Ge-
wichtsfunktionen Wi(x) die gleichen Funktionen eingesetzt, die in derNaherungslosung verwendet worden sind. Diese Vorgehensweise ist die
Basis der FE-Methode fur Problemen, die die erste Ableitung beinhalten.Furselbstadjugierte Dgln liefert diese Methode Ergebnisse, die mit denenaus Variationsmethode identisch sind. Das Galerkin-Verahren wird fur
die Entwicklung der FE-Gleichungen bei Feldproblemen in diesem Kursverwendet.
Verfahren der kleinsten Quadrate : Diese Verfahren setzt den Restals Gewichtsfunktion ein und erhoht einen neuen Fehler, der wie folgt
definiert wird:
Er =
∫ H
0
[R(x)]2 dx (1.5)
Dieser Fehler muss nun in Bezug auf unbekannten Koeffizienten der Nahe-rungslosung minimiert werden (Ableitung nach Koeffizienten gleich Null
setzen). Dieses Verfahren wurde zur Formulierung der FE-Methode ver-
1.2. INTEGRALFORMULIERUNGEN FUR DIE NUMERISCHE LOSUNG 7
wendet, aber hat sich gegen Galerkin-Verfahren und Variationsverfahren
nicht durchgesetzt.
1.2 Integralformulierungen fur die numerische
Losung
Nun bleibt zu zeigen, wie die genannten Integralmethoden zur Gewinnungeiner Naherungslosung fur ein physikalisches Problem eingesetzt werden.
Als Beispiel nehmen wir einen einfach gelagerten Balken mit einem Ein-zelmoment an jedem Ende als Belastung. Der Balken und der Verlauf desBiegemoments sind in Bild ?? dargestellt. Wir suchen die Verformung des
Balkens y(x) (als unbekannte Große) unter der Last. Das System gehorchtder folgenden Dgl
EId2y
dx2−M(x) = 0 , (1.6)
die die Verformung mit der Last verknpft und den Rabdbedingungen
y(0) = y(H) = 0 unterliegt. Der Koeffizient EI beschreibt die Biegesteifig-keit des Balkens und ist bekannt (E: Elastizitatsmodul und I: Tragheits-
moment). M(x) stellt den Verlauf des Biegemoments entlang des Balkensdar (hier ist berall M(x) = M0). Auf der Suche nach einer Losung fur
y(x) kann man die folgende Gleichung als eine Naherung fur ihren Verlaufwahlen:
y(x) = A sinπx
H(1.7)
wobei A ein unbekannter Koeffizient ist. Diese Gleichung erfllt die Rand-
bedingungen y(0) = y(H) = 0 und weist eine ahnliche Form wie dieVerformung des Balkens auf. Sie ist daher eine akzepable Kandidatin fur
eine Naherung. Die exakte Losung der Dgl 1.6 lautet:
y(x) =M0x
2EI(x−H) (1.8)
8 KAPITEL 1. GRUNDKONZEPTE
1.2.1 Variationsmethode
Der Funktionalausdruck fur die Dgl 1.6 ist:
Π =
∫ H
0
[
EI
2
(
dy
dx
)2
+M0y
]
dx (1.9)
Die Gl. 1.7 stellt sich als beste Naherung heraus, wenn ein bestimmter Be-trag von A den Wert von Π minimiert. Um dieses A ermitteln zu konnen,
mussen wir zuerst Π als eine Funktion von A umschreiben und dann siebezogen auf A minimieren. Mit der Beziehung
dy
dx=Aπ
Hcos
πx
H
bekommen wir
Π =
∫ H
0
[
EI
2
(
Aπ
Hcos
πx
H
)2
+M0A sinπx
H
]
dx
oder
Π =
(
EIπ2
4H
)
A2 +
(
2M0H
π
)
A (1.10)
Die Minimierung von Π bedeutet
∂Π
∂A= 2
(
EIπ2
4H
)
A+2M0H
π= 0 (1.11)
und ergibt
A = −4M0H
2
π3EI(1.12)
Die Naherungslosung sieht dann so aus:
y(x) = −4M0H
2
π3EIsin
πx
H(1.13)
1.2.2 Punktweise Kollokation
Die punktweise Kollokation verlangt, dass der Rest an bestimmter An-
zahl (gleich der unbekannten Koeffizienten der Naherungsfunktion) von
1.2. INTEGRALFORMULIERUNGEN FUR DIE NUMERISCHE LOSUNG 9
Punkten im Problembereich gleich Null wird. Der Rest wird durch Ein-
setzen von Naherungsfunktion (Gl. 1.7) in die Dgl des Problems (Gl. 1.6)berechnet:
R(x) = EId
dx2
(
A sinπx
H
)
−M0 = −EIAπ2
H2sin
πx
H−M0 (1.14)
Da es nur einen unbekannten Koeffizient, A, gibt, wird R(x) irgendwozwischen 0 und H gleich Null gesetzt. Einfachheithalber nehmen wir x =H2
und bekommen damit
R(H
2) = −EI
Aπ2
H2sin
π
2−M0 = 0
und
A = −M0H
2
EIπ2. (1.15)
Damit ergibt sich die Naherungslosung zu
y(x) = −M0H
2
EIπ2sin
πx
H(1.16)
Wenn wir einen anderen Punkt statt x = H2
wahlen, bekommen wir eine
andere Naherungslosung.
1.2.3 Bereichsweise Kollokation
Dieses Kollokationsverfahren verlangt, dass das Integral des Restes ber je-
den Teilbereich von n-Teilbereichen, in denen das Problemgebiet aufgeteiltwurde, gleich Null ist, wobei n die Anzahl der unbekannten Koeffizienten
der Naherungslosung darstellt. Die Große der Teilbereiche steht dem An-wender frei zu wahlen. Hier haben wir nur einen Unbekannten und daherbrauchen wir auch nur einen Teilbereich und der muss [0, H] sein. Der
Rest ist identisch mit Gl. 1.14; daher∫ H
0
R(x)dx =
∫ H
0
[
−EIAπ2
H2sin
πx
H−M0
]
dx = 0
Integrieren liefert
−
(
2EIπ
H
)
A−M0H = 0 und A = −M0H
2
2πEI(1.17)
10 KAPITEL 1. GRUNDKONZEPTE
Die Naherungslosung lautet somit:
y(x) = −M0H
2
2πEIsin
πx
H(1.18)
1.2.4 Galerkin-Verfahren
Bei Galerkin-Verfahren wird der Integralausdruck∫
Wi(x)R(x)dxmit den
gleichen Funktionen fur Wi(x) ausgewertet, die fur die Naherungslosungverwendet wurden. In diesem Beispiel gibt es nur eine Gewichtsfunktion,
Wi(x) = sinπx/H. Der Ausdruck fur den Rest ist wieder die Gl 1.14 undder Funktionalausdruck lautet:
∫ H
0
sinπx
H
[
−EIAπ2
H2sin
πx
H−M0
]
dx = 0
Das Integrieren ergibt:
−EIπ2A
2H+
2M0H
π= 0
Losung dieser Gleichung fhrt zu
A = −4M0H
2
π3EI(1.19)
und die Naherungslosung ist somit
y(x) = −4M0H
2
π3EIsin
πx
H(1.20)
Diese Losung ist mit der Losung aus Variationsmethode identisch.
1.2.5 Verfahren der kleinsten Quadrate
Bei diesem Verfahren wird ein Neuer Fehlerausdruck, Er =∫
[R(x)]2 dx,aufgestellt. Einsetzen des Restes in diesen Ausdruck ergibt
Er =
∫ H
0
[
−EIπ2
H2A sin
πx
H−M0
]2
dx .
1.3. POTENTIALENRGIE-FORMULIERUNG 11
Die Integration liefert
Er =A2H
2
(
EIπ2
H2
)2
+4M0EIπ
HA+M2
0H
Die Minimierung des Fehlers in Bezug auf A bedeutet
∂Er
∂A= AH
(
EIπ2
H2
)2
+4M0EIπ
H= 0 (1.21)
Nach der Losung fur A ergibt sich die Naherungsloung zu
y(x) = −4M0H
2
π3EIsin
πx
H. (1.22)
Diese Losung ist identisch mit den Losungen aus Variationsmethode und
Galerkin-Verfahren. Das ist nicht moglich zu sagen, welches Verfahrenist das genaueste. Der Fehler hangt von der Wahl der Naherungsfunk-
tion und der zu losenden Gleichung. Verlaufe der Fehlerquoten von un-terschiedlichen Verfahren wurden in Bild ?? dargestellt. Anscheinend istdie Naherung mit der Gl. 1.22 genauer als mit den Gln. 1.16 und 1.18,
obwohl es moglich ist, eine Kollokationspunkt zu finden, der eine Verfor-mung liefert, die der exakten Losung am besten entspricht. Die gewahl-
ten Kollokationspunkte oder -teilbereiche beeinflssen die Genauigkeit derNaherungslosung. Der wichtige Punkt aus diesem Beispiel ist, dass die nu-
merische Losung einer Dgl in einem Integralausdruck formuliert werdenkann. Die Integralformulierung ist die Grundeigenschaft der FE-Methode.
1.3 Potentialenrgie-Formulierung
Probleme aus dem Bereich der elastischen Festkorpermechanik, wieBalken-, Platten-, Scheiben- und Schalenprobleme, konnen auf verschie-
denen Weisen gelost werden. Wenn eine klassische Losung (Aufstellungder Dgl. vom Problem und ihre analytische Losung) aufgrund schwierige
mathematische Beschreibung der Aufgabe nicht moglich ist, wird oft eineAlternativlosung verwendet. Diese Alternative basiert auf einem Konzept,
das besagt, dass ein Korper im Gleichgewichtzustand solche Verformungenhat, die einem Minimum von Potentialenergie entsprechen. Dieses Kon-
zept wird auch das Prinzip des Minimums von Potentialenergie genannt.
12 KAPITEL 1. GRUNDKONZEPTE
Die Potentialenergie besteht zum Teil aus der Verzerrungsenergie, die in-
folge der Verformungen des Korpers entsteht und im Korper gespeichertwird. Verzerrungsenergie ist ein Volumenintegral, das das produkt von
Spannungs- und Verzerrungskomponenten beinhaltet. Zum Beispiel lasstsich fur einen Stab die Verzerrungsenergie wie folgt schreiben:
Λ =
∫
V
σxxεxx
2dV . (1.23)
Uber die Verzerrungsenergie und das Prinzip des Minimums von Poten-
tialenergie werden wir in weiteren Kapiteln ausfhrlicher reden. Wichtigwar hier die Vorstellung, dass die Verformungsanalyse der Struktur- und
Festkorpermechanik die Verzerrungsenergie mit einem Minimierungspro-zess kombiniert. Aus rechnerischer Sicht sieht die Analyse eines Fach-
werks oder einer Platte mit Hilfe dieses Prinzips sehr ahnlich aus wie dieVariations- oder Galerkin-Methode. Diese ahnlichkeit lasst sich jedoch erstspater offenbaren.
1.4 Die Methode der Finiten Elemente
Diese Methode ist eine numerische Prozedur fur die Losung physikali-
scher Probleme, die einer Dlg oder einem Energietheorem gehorchen. Siehat zwei Kennzeichen, die sie von anderen numerischen Methoden untzer-scheidet:
1. Sie setzt eine Integralformulierung fur die Aufstellung eines Systems
von algebraischen Gleichungen ein.
2. Sie nutzt kontinuierliche, stckweise glatte Funktionen um die unbe-
kannte(n) Große(n) zu approximieren.
Das zweite Kennzeichen unterscheidet die FEM von anderen numerischen
Methoden, die ebenfalls eine Integralformuliereung zu Grunde legen. Erin-nert wird an die Naherungslosung im vorigen Abschnitt, y = A sinπx/H.
Diese Funktion hat unendlich kontinuierliche Ableitungen. Fur den Ein-satz in FEM braucht eine kontinuierliche Funktion nur soviel Kontinuitat
in ihrer Ableitungen vorzuweisen, die fur die Auswertung der Integrale
1.4. DIE METHODE DER FINITEN ELEMENTE 13
notwendig sind. Fur eine Integralformulierung wie die Variationsmethode
(s. Gl. 1.9) ist keine Kontinuitat der ersten Ableitung notig. Das Integralkann auch dann ausgewertet werden, wenn die erste Ableitung stckweise
kontinuierlich ist. Eine Gleichung zusammengesetzt aus mehreren linea-ren Teilen kann da genauso gut als Naherungsfunktion verwendet werden.Um die FE-Methode einsetzen zu konnen, muss man das Problemgebiet in
kleine Teilbereiche aufteilen, die sich an gemeinsamen Grenzen und Kno-ten trefen. Dieser Schritt wird Diskretisierung genannt. Dann wird das
Problemgebiet mit einem Netz bespannt, dessen Netzaugen mit den Teil-bereichen Deckungsgleich sind. Dieser Schritt ist als Vernetzung bekannt
und kann unter Umstanden zugleich die Diskretisierung beinhalten. EinFE-Modell der Balkenaufgabe kann wie in Bild ?? aussehen. Das Bild ??a
zeigt eine Diskretisierung des Balkens in 6 Teilbereiche und eine Vernet-zungt mit 6 Elementen, die jeweils zwei Knoten besitzen und durch lineareTerme der Konotenwerte, y = f(x), definiert sind. Die Verformung des
Balkens wird durch Geradensegmente approximiert. Das Bild ??b zeigteine Diskretisierung des Balkens in 3 Teilbereiche und eine Vernetzung
mit 3 Elementen, die jeweils drei Knoten besitzen und durch quadratischeTerme der Knotenwerte, y = g(x), definiert sind. Die Verformung des Bal-
kens wird diesmal durch quadratische Teilkurven approximiert. Keine derbeiden Gleichung, y = f(x) oder y = g(x), hat zwischen zwei benach-barten Elementen eine kontinuierliche erste Ableitung. Funktionen ohne
kontinuierliche erste Ableitung konnen auch in Galerkin-Verahren verwen-det werden. In diesem Fall muss die Integration des Termes mit zweiter
Ableitung, d2y/dx2, stuckweise geschehen.Die FE-Methode kann in funf Schritten aufgeteilt werden, die heir aufge-
listet und in den nachsten Kapiteln erlautert werden:
1. Diskretisierung des Gebiets. Das bedeutet die Lage der Knoten fest-
zulegen (inkl. Koordinaten) und sie zu numerieren.
2. Die Naherungsfunktion bestimmen. Das bedeutet die Ordnung der
Approximation, linear oder quadratisch oder hohere Ordnung, fest-zulegen und fur jedes Element eine Gleichung in Termen der unbe-
kannten Knotenwerte zu schreiben.
3. Das Gleichungssystem fur das Gesamtgebiet aufstellen. Bei Nutzung
14 KAPITEL 1. GRUNDKONZEPTE
der Galerkin-Methode wird die Gewichtsfunktion fur jede unbkann-
te Knotengroße definiert, dann wird das Restintegral ausgwertet.Dies erzeut fur jede unbekannte Knotengroße eine Gleichung. Bei
Potentialenergie-Formulierung wird die Potentialenergie des Systemsin Termen von Knotenverschiebungen geschrieben und dann mini-miert. Dies ergibt fur jede unbekannte Knotenverschiebung eine Glei-
chung.
4. Losung des Gleichungsystems
5. Berechnung der gesuchten Großen. Diese Großen sind normalerweisemit den Ableitungen der Parameter verwandt.
Kapitel 2
Eindimensionales lineares Element
Hier diskutieren wir die Aufteilung eines eindimensionalen Gebiets inlinearen Elementen und die Entwicklung einer Elementgleichung. Die
Elementgleichung wird dann verallgemeinert, sodass eine kontinuierlichestuckweise glatte Gleichung fur das Gebiet geschrieben werden kann. Das
lineare Element wird benuzt um eine Naherungslosung fur die Dgl.
Dd2φ
dx2+Q = 0 (2.1)
zu erreichen. Dieses Element wird spater benutzt um Verschiebungen imGebiet zu berechnen.
2.1 Aufteilung des Gebiets in Elemente
Das eindimensionale Gebiet ist ein Liniensegment und die Aufteilung in
Teilbereichen ist einfach. Das Liniensegment wird durch Knotenpukte inkleinere Segmente geteilt. Danach folgt eine Nummerierung der Knoten
und Elemente, die an Knoten mitenander verbunden sind (s. Bild ??).Hier ware angebracht auf einige Regeln zu achten:
• In Bereichen des Gebiets, wo sich unbekannte Großen wahrscheinlich
starker andern, sind die Knoten dichter zu einander zu legen, damitdie rapide Anderungen besser erfasst werden konnen. Hier mischt sich
die Erfahrung in den Losungsprozess ein.
15
16 KAPITEL 2. EINDIMENSIONALES LINEARES ELEMENT
• Wo es scharfe Anderung in den Werten der Koeffizienten D und Q in
Dgl. (2.1) gibt, mussen Knoten platziert werden (in eindimensiona-lem Fall ein Knoten). Dies vereinfacht die Auswertung der Integral-
ausdrucke mit diesen Termen.
• Dort, wo ein numerischer Wert von φ gesucht wird, muss auch ein
Knoten vorhanden sein.
2.2 Das lineare Element
Das eindimensionale lineare Element ist ein Liniensegment mit der Lnge
L und zwei Knoten, jeder an einem Ende (s. Bild ??). Die Knoten werdenmit i und j bezeichnet und die Knotengroßen durch Φi und Φj dargestellt.Der Parameter φ1 hat zwischen den Knoten einen linearen Verlauf, und
damit hat die Ansatzfunktion fur φ die folgende Form
φ = a1 + a2x . (2.2)
Die Koeffizienten a1 und a2 konnen durch geltende Bedingungen an Kno-
ten bestimmt werden:
φ = Φi an x = Xi
φ = Φj an x = Xj (2.3)
Und damit hat man
Φi = a1 + a2Xi
Φj = a1 + a2Xj , (2.4)
woraus sich die unbekannten Koeffizienten a1 und a2 berechnet lassen:
a1 =ΦiXj − ΦjXi
Xj −Xi
a2 =Φj − Φi
Xj −Xi. (2.5)
1Das Symbol φ wird uberall in diesem Text als allgemeine skalare Große verwendet. Symbole inGroßbuchstaben wie X , Y , Φ und U bezeichnen die Knotengroßen
2.2. DAS LINEARE ELEMENT 17
Einsetzen von a1 und a2 in die Ansatzfunktion (2.2) ergibt:
φ =
(
Xj − x
L
)
Φi +
(
x−Xi
L
)
Φj , (2.6)
wobei Xj−Xi durch die Elementlange L ersetzt wurde. Die Gl. (2.6) stellteine Standardform der finiten Elemente dar. In dieser Gl. sind die Kno-
tengroße mit linearen Funktionen von x multipliziert, die Formfunktion
oder Interpolationsfunktion genannt werden. Formfunktionen werden mit
N bezeichnet und ein Indiz zeigt mit welchem Knoten eine Formfunkti-on assoziiert ist. Somit stehen die Formfunktionen fur die Verknopfung
der Knotengroßen mit der Große φ an beliebiger Stelle x innerhalb desElementes. Die Formfunktionen in (2.6) werden mit Ni und Nj wie folgtbezeichnet
Ni =Xj − x
Lund Nj =
x−Xi
L(2.7)
Damit lasst sich die Gl. (2.6) wie folgt umschreiben
φ = NiΦi +NjΦj (2.8)
oder in Matrizendarstellung
φ = [N ]{Φ} , (2.9)
wobei [N ] = [Ni Nj] ein Zeilenvektor mit Formfunktionen ist und
{Φ} =
{
Φi
Φj
}
ein Spaltenvektor mit Knotengroßen des Elements. Jede Formfunktionerhalt an ihrem assoziierten Knoten den Wert 1 und an dem andren Kno-
ten den Wert 0 (s. Bild ??). Die Summe der beiden Formfunktionen ergibtden Wert 1. Die dritte Eingenschaft der Formfunktionen ist, dass sie im-
mer Polynomen desselben Typs sind wie die ursprungliche Ansatzfunktion.Die Ansatzfunktion (2.2) ist linear, so sind auch die Formfunktionen. Ware
die Ansatzfunktion quadratisch gewesen, so waren unsere Formfunktionenauch quadratisch. Eine Weitere Eigenschaft der Formfunktionen ist, dassdie Summe ihrer Ableitungen in Bezug auf x den Wert Null ergibt.
Beispiel 2.1
18 KAPITEL 2. EINDIMENSIONALES LINEARES ELEMENT
2.3 Eine kontinuierliche stuckweise glatte Funktion
Eine kontinuierliche stuckweise glatte Funktion fur ein eindimensionalesGebiet kann konstruiert werden, indem mann mehrere lineare Funktio-
nen mit den genannten Eigenschaften im vorigen Abschnitt miteinanderverknopft. Jede solcher Funktionen kann man wie folgt darstellen
φ(e) = N(e)i Φi +N
(e)j Φj , (2.10)
wobei
N(e)i =
Xj − x
Xj −Xiund N
(e)j =
x−Xi
Xj −Xi(2.11)
Der Superskript (e) kennzeichnet eine Elementgroße. Diese sind alle
notig, um jedem Element die richtigen Werte von i, j und e zuweisen zukonnen. Die mit einem Element e korrespondierenden Werte von i und j
erhalt man aus dem Netz der Elemente
e i j
1 1 2
2 2 33 3 4
4 4 5
Die Gleichung for jedes Element in Bild ?? ist
φ(1) = N(1)1 Φ1 +N
(1)2 Φ2
φ(2) = N(2)2 Φ2 +N
(2)3 Φ3
φ(3) = N(3)3 Φ3 +N
(3)4 Φ4 (2.12)
φ(4) = N(4)4 Φ4 +N
(4)5 Φ5
Beachte, dass die Gln. N(1)2 und N
(2)2 unterschiedlich sind, obwohl beide
den Knoten 2 beinhalten. Diese zwei Gleichungen sind
N(1)2 =
x−X1
X2 −X1und N
(2)2 =
X3 − x
X3 −X2
Es ist zu beachten, dass jede Gl. in (2.12) steht nur fur ein einziges Element
und ist außerhalb dieses Element nicht gultig. Die erste Gl. gilt damit nur
2.4. EIN BEISPIEL ZUR FE-FORMULIERUNG 19
fur X1 ≤ x ≤ X2. Die Geltungsbereiche der Gln werden i.d.R. nicht
explizit genannt.
Ein Kommentar uber die Notation
Die Bezeichnung der Elementgroßen wird in diesem Text sehr oft notigsein. Die folgende Notationsregel wird benutzt, sodass ein Superskript (e)
nicht auf jeden Koeffizient zu platzieren ist.
1. Wenn Klammern einen Superskript (e) haben, wie z.B. (Gφ+Q)(e),das bedeutet, dass jeder Term im Klammer auf Elementbasis zu in-
terpretieren ist.
2. Eine Große auf der linken Seite des Gleichzeichens mit einem Su-
perskript (e) bedeutet, dass die Großen auf der rechten Seite desGleichzeichens auf ein bestimmtes Element beziehen. Z. B.
φ(e) = NiΦi +NjΦj
bedeutet, dass Ni und Nj in der Tat N(e)i und N
(e)j sind und Φi und
Φj die Knotengroßen des Elements.
2.4 Ein Beispiel zur FE-Formulierung
In diesem Abschnitt wird die FE-Methode anhand der Entwicklung einer
Naherungslosung fur eine eindimensionale Dgl. gezeigt. Dabei setzen wirdie Erkenntnisse uber Formfunktionen aus vorigem Abschnitt ein. Die zulosende Dgl.
Dd2φ
dx2+Q = 0 (2.13)
soll folgende Randbedingungen erfullen
φ(0) = φ0 und φ(H) = φH (2.14)
Zwei physikalische Probleme lassen sich durch die Dgl. (2.13) beschrei-ben: 1. Die Verformung eines einfach gelagereten Balkens mit bekanntem
Momentenverlauf und 2. Warmeleitung durch eine Verbundwand mit be-kannter Temperatur auf beiden Seiten. Die FE-Gleichungen werden mit-
tels Galerkin-Methode aufgestellt. Die Auswertung des Restintegralsfuhrt
20 KAPITEL 2. EINDIMENSIONALES LINEARES ELEMENT
zu einer Knotengleichung, die wir zur Aufstellung des Gleichungssytems
auf rekursive Art nutzen. Die Knotengleichung wird dann zur Losung einesBalkenverformungsproblems verwendet.
2.4.1 Gewichtsfunktionen
Ein System von linearen Gleichungen wird durch Auswertung des Restin-
tegrals 2
−
∫ H
0
W (x)
(
Dd2φ
dx2+Q
)
dx = 0 (2.15)
fur jeden Knoten aufgestellt. Dabei wird fur jeden Knoten eine neue Ge-wichtsfunktion verwendet. Im Integralausdruck ist φ die unbekannte Nahe-
rungslosung.Galerkin-Verfahren bei gewichteten Resten verlangt, dass die Gewichts-
funktionen mit Formfunktionen Ni und Nj konstruiert werden. Die Ge-wichtsfunktionen werden in der Galerkin’schen FE-Formulierung wie folgt
definiert: Die Gewichtsfunktion fur den s-ten Knoten, Ws, beinhaltet die
Formfunktionen, die mit dem Knoten s assoziiert sind. Die Gewichtsfunk-tion fur den dritten Knoten des linearen Netzes (s. Bild ??) beinhaltet die
Formfunktionen fur diesen Knoten:
W3(x) =
{
N(2)3 X2 ≤ x ≤ X3
N(3)3 X3 ≤ x ≤ X4
(2.16)
Im Allgemeinen gilt
Ws(x) =
{
N(e)s Xr ≤ x ≤ Xs
N(e+1)s Xs ≤ x ≤ Xt
(2.17)
Die Gweichtsfunktionen fur den ersten und den letzten Knoten sind in
Bild ??a und c dargestellt. Die entsprechenden Gleichungen sind:
W1(x) = N(1)1 und Wp(x) = N
(P−1)P (2.18)
Die Gewichtsfunktionfur jeden Knoten beinhaltet entweder Ni oder Nj,oder eine Kombination der beiden.
2Das Integral wurde mit −1 multipliziert, um die Ergebnisse in einfacherer Form darstellen zu konnen
2.4. EIN BEISPIEL ZUR FE-FORMULIERUNG 21
2.4.2 Das gewichtete Restintegral
Nach Definition der Wichtungsfunktionen kommt die Auswertung des Re-
stintegrals (2.15). Mit Berucksichtugung der Knotenreihenfolge r, s und tin Bild ??b bekommen wir fur das Restintegral am Knoten s
Rs = R(e)s +R(e+1)
s = −
∫ Xs
Xr
[
Ns
(
Dd2φ
dx2+Q
)](e)
dx
−
∫ Xt
Xs
[
Ns
(
Dd2φ
dx2+Q
)](e+1)
dx , (2.19)
weil fur x < Xr und x > Xt gilt: Ws = 0. Das Restintegral ist in zwei
Teilen geteilt, weil Ws(x) durch zwei separate Funktionen innerhalb des
Intervals Xr ≤ x ≤ Xt definiert worden ist. Die Termen R(e)s und R
(e+1)s
stellen die Beitrage der Elemente (e) und (e+ 1) zum Rest am Knoten sdar.Es gibt ein Problem bei jedem Integral in Gl. (2.19). Die Naherungslosung
ist in ihrer ersten Ableitung, dφ/dx, nicht kontinuierlich und daher ist dasIntegral von d2φ/dx2 nicht definiert. Dieses Problem kann man umgehen,
indem man den Term d2φ/dx2 in einen neuen Term umschreibt. Wir be-trachten zuerst das erste Integral in Gl. (2.19) und weisen darauf hin,
dass
d
dx
(
Nsdφ
dx
)
= Nsd2φ
dx2+dNs
dx
dφ
dx. (2.20)
Daher gilt
Nsd2φ
dx2=
d
dx
(
Nsdφ
dx
)
−dNs
dx
dφ
dx. (2.21)
einsetzen ins erste Integral ergibt
−
∫ Xs
Xr
(
NsDd2φ
dx2
)(e)
dx = −
(
DNsdφ
dx
)(e)∣
∣
∣
∣
∣
Xs
Xr
+
∫ Xs
Xr
(
DdNs
dx
dφ
dx
)(e)
dx (2.22)
22 KAPITEL 2. EINDIMENSIONALES LINEARES ELEMENT
Ahnliche Operationen fur den ersten Term des zweiten Integrals in Gl.
(2.19) ergibt
−
∫ Xt
Xs
(
NsDd2φ
dx2
)(e+1)
dx = −
(
DNsdφ
dx
)(e+1)∣
∣
∣
∣
∣
Xt
Xs
+
∫ Xt
Xs
(
DdNs
dx
dφ
dx
)(e+1)
dx (2.23)
Der erste Term in Gl. (2.22) und in Gl. (2.23) lasst sich vereinfachen, weil
die Formfunktionen an entsprechenden Knoten entweder 0 oder 1 sind.Die komplette Gleichung des Restes ergibt sich zu
Rs = R(e)s + R(e+1)
s = −
∫ H
0
Ws
(
Dd2φ
dx2+Q
)
dx
= −
(
Ddφ
dx
)(e)∣
∣
∣
∣
∣
x=Xs
+
∫ Xs
Xr
(
DdNs
dx
dφ
dx−NsQ
)(e)
dx (2.24)
+
(
Ddφ
dx
)(e+1)∣
∣
∣
∣
∣
x=Xs
+
∫ Xt
Xs
(
DdNs
dx
dφ
dx−NsQ
)(e+1)
dx = 0
Die beiden Terme, die an der Stelle x = Xs ausgewertet werden, legen eineinterelemntare Forderung fest. Der Rest kann nicht Null werden, solange
die Differenz zwischen diesen beiden Termen nicht verschwindet.
2.4.3 Auswertung der Integrale
Die Auswertung der Integrale in Gl. (2.24) fuhrt zu der Restgleichung fur
einen Knoten im Innenbereich. Fangen wir mit dem Element (e) (s. Bild??b) an. Unter Berucksichtigung der Gl. (2.6) haben wir
φ(e) = NrΦr +NsΦs
φ(e) =
(
Xs − x
L
)
Φr +
(
x−Xr
L
)
Φs . (2.25)
Damit bekommen wir fur den Knoten s
N (e)s =
x−Xr
L;
dN(e)s
dx=
1
L(2.26)
2.4. EIN BEISPIEL ZUR FE-FORMULIERUNG 23
unddφ(e)
dx=
1
L(−Φr + Φs) (2.27)
Einsetzen der entsprechenden Terme und Auswertung der Integrale ergibt
∫ Xs
Xr
DdNs
dx
dφ
dxdx =
D
L(−Φr + Φs) (2.28)
und∫ Xs
Xr
QNsdx =QL
2. (2.29)
Fur die Berechnung von R(e)s sind die Integrale (2.28 und 2.29) mit dem
interelemntaren Beitrag fur das Element (e) zu kombieren:
R(e)s = −
(
Ddφ
dx
)(e)∣
∣
∣
∣
∣
x=Xs
+D
L(−Φr + Φs) −
QL
2. (2.30)
Wir konnen nun mit dem Element (e+1) fortfahren (s. Bild ??b) und daszweite Integral in Gl. (2.24) auswerten
φ(e+1) = NsΦs +NtΦt
φ(e+1) =
(
Xt − x
L
)
Φs +
(
x−Xs
L
)
Φt (2.31)
Daher bekommen wir
N (e+1)s =
Xt − x
L;
dN(e+1)s
dx= −
1
L(2.32)
unddφ(e+1)
dx=
1
L(−Φs + Φt) (2.33)
Damit ergibt die Auswertung der Integrale
∫ Xt
Xs
DdNs
dx
dφ
dxdx =
D
L(Φs − Φt) (2.34)
∫ Xt
Xs
QNsdx =QL
2. (2.35)
24 KAPITEL 2. EINDIMENSIONALES LINEARES ELEMENT
Das Element (e + 1) leistet somit den folgenden Beitrag zum Rest am
Knoten s
R(e+1)s = D
dφ
dx
∣
∣
∣
∣
x=Xs
+D
L(Φs − Φt) −
QL
2(2.36)
Setzen wir die beiden Beitrage R(e)s und R
(e+1)s zusammen und dann gleich
Null, erhalten wir die Restgleichung fur den Knoten s:
Rs =
(
Ddφ
dx
)(e+1)∣
∣
∣
∣
∣
x=Xs
−
(
Ddφ
dx
)(e)∣
∣
∣
∣
∣
x=Xs
−
(
D
L
)(e)
Φr +
[
(
D
L
)(e)
+
(
D
L
)(e+1)]
Φs −
(
D
L
)(e+1)
Φt
−
(
QL
2
)(e)
−
(
QL
2
)(e+1)
= 0 (2.37)
Die ubliche Losungsprozedur ist die Aufstellung des Gleichungssystems
ohne die interelemntaren Terme. Nach der Losung des Gleichungssytemslasst sich der Ausdruck
(
Ddφ
dx
)(e+1)∣
∣
∣
∣
∣
x=Xs
−
(
Ddφ
dx
)(e)∣
∣
∣
∣
∣
x=Xs
(2.38)
einfach rechnen. Theoretisch konnte der Wert des Ausdrucks (2.38) alsein Maß fur die Qualitat des Netzes werden oder als ein Indikator, der
bestimmt wo das Netz zu verfeinern ist.Wenn D(e) = D(e+1) ist, dann reduziert sich die interelemntare Forderungauf das Verschwinden von
(
dφ
dx
)(e+1)∣
∣
∣
∣
∣
x=Xs
−
(
dφ
dx
)(e)∣
∣
∣
∣
∣
x=Xs
. (2.39)
Ohne Kontinuitat in der ersten Ableitung wird diese Differrenz (2.39) nicht
verschwinden, und die Kontinuitat kann niemals mit linearen Elementenerreicht werden, es sei denn, die Losung ist eine gerade Linie. Der Betrag
von (2.39) wird kleiner je das Netz der finiten Elemente feiner wird, aberer wird nie fur alle Knoten Null sein. Die Differenz (2.39) kann als einen
Fehler angesehen werden, der im Gleichungssytem nicht integriert ist, aber
2.4. EIN BEISPIEL ZUR FE-FORMULIERUNG 25
uns stets daran erinnert, dass die Losung eine Naherung ist.
Wenn wir diese interelementare Forderung aus der Gl. (2.37) heraus neh-men, eine aufeinander folgende Numerierung der Knoten und Elemente
vornehmen und alle Ausdrucke in Termen von s schreiben, dann reduziertsich die Restgleichung fur den Knoten s auf
Rs = −
(
D
L
)(s−1)
Φs−1 +
[
(
D
L
)(s−1)
+
(
D
L
)(s)]
Φs −
(
D
L
)(s)
Φs+1
−
(
QL
2
)(s−1)
−
(
QLK
2
)(s)
= 0 (2.40)
2.4.4 Analyse eines Biegebalkens
Hier wollen wir die allgemeine Restgleichung (2.40), die wir durch Aus-wertung des Restintegrals gewonnen haben, fur die Approximation von
Verfgormungen eines einfach gelagerten Balkens einsetzen. Der Balken istim mittleren Bereich verstarkt worden und wird durch gleich große jedochgegensinnige Momente an beiden Enden belastet. Detailangaben sind dem
Bild (??) zu entnehmen. Die herrschende Dgl. fur die Verformung ist
EId2y
dx2−M(x) = 0 , (2.41)
wobei EI die Biegesteifigkeit, y die Verformung und M(x) den Verlauf des
Biegemomentes darstellt und die Randbedingung in y(0) = y(800) = 0erfasst sind.ein Vergleich mit der Gl. (2.13) zeigt, dass D = EI und Q = −M(x) =
−106 ist. Die Verstarkung in der Mitte erzeugt einen Sprung in EI bzw.in D am anfang und Ende des verstarkten Bereiches. Daher setzen wir bei
Diskretisierung Knoten an diesen Punkten. Das Interesse an Verformungin der Mitte des Balkens ist Grund genug fur einen weiteren Knoten in
der Mitte des Balkens. Somit haben wir die einfachste Diskretisierungund Vernetzung mit 5 Knoten und 4 Elementen. Die Elementdaten sind
in dieser Tabelle zusammengefasst
26 KAPITEL 2. EINDIMENSIONALES LINEARES ELEMENT
e D Q L
1 2,4e10 -1,0e6 200
2 4,0e10 -1,0e6 2003 4,0e10 -1,0e6 200
4 2,4e10 -1,0e6 200
Da Q und L konstant bleiben, lasst sich Gl. (2.40) wie folgt vereinfachen
Rs =−Ds−1Ys−1 + (Ds−1 +Ds)Ys −DsY(s+ 1)
L−QL = 0 , (2.42)
wobei Y fur die Knotenverschiebungen steht (anstelle des Φ). Schreibenwir die Restgleichungen fur die mittleren Knoten 2, 3 und 4, bekommen
wir
R2 = −1, 2Y1 + 3, 2Y2 − 2, 0Y3 + 2 = 0
R3 = −2, 0Y2 + 4, 0Y3 − 2, 0Y4 + 2 = 0
R4 = −2, 0Y3 + 3, 2Y4 − 1, 2Y5 + 2 = 0 . (2.43)
Durch Einarbeiten der Randbedingungen, Y1 = Y5 = 0, erhalten wir
R2 = 3, 2Y2 − 2, 0Y3 = −2
R3 = −2, 0Y2 + 4, 0Y3 − 2, 0Y4 = −2
R4 = −2, 0Y3 + 3, 2Y4 = −2 . (2.44)
Losung dieses Gleichungssystems ergibt: Y2 = −2, 5 cm, Y3 = −3, 0 cmund Y4 = −2, 5 cm. Nun kann man anhand der bekannten Knotenver-
schiebungen andere Großen berechnen, z.B. die Verformung an der Stellex = 300 cm, oder die Neigung des Balkens am Knoten 1.
Fur die Verformung an x = 300 cm mussen wir das Element 2 betrachten,weil sich der Punkt innerhalb dieses Elementes befindet. Furs Element 2gilt
y(2) = N(2)2 Y2 +N
(2)3 Y3 =
(
X3 − x
X3 −X2
)
Y2 +
(
x−X2
X3 −X2
)
Y3
Einsetzen von x = 300 cm und allen anderen Werten ergibt
y =
(
400 − 300
400 − 200
)
(−2, 5) +
(
300 − 200
400 − 200
)
(−3, 0) = −2, 75 cm
2.4. EIN BEISPIEL ZUR FE-FORMULIERUNG 27
Fur die Berechnung der Neigung am Knoten 1 mussen wir die erste Ab-
leitung der Verformungsfunktion (2.6) heranziehen und zwar fur das ersteElement:
y(1) = N(1)1 Y1 +N
(1)2 Y2 =
(
X2 − x
X2 −X1
)
Y1 +
(
x−X1
X2 −X1
)
Y2
dy(1)
dx=
1
L(−Y1 + Y2) =
−2, 5
200= −0, 0125 cm/cm
Wie man sieht ist die Neigung fur das Element 1 konstant (unabhangig
von x innerhalb des Elementes). Dieses gilt auch fur andere Elemente,obwohl sich die Betrage unterscheiden. Konstante Neigung ist ein großer
Nachteil von linearen Elementen.Das Balkenproblem wurde hier wegen seiner Dgl. gewahlt und soll nicht
den Eindruck erwecken, dass alle Balkenaufgaben auf dieser Weise zu losensind. Es gibt eine leistungsfahigere Methode mit einem spezifischen Bal-
kenelement, die wir in dem strukturmechanischen Kapitel kennen lernen.
2.4.5 Matrizendarstellung
Es ist einfacher und ordentlicher bei der Behandlung von Gleichungssyste-men (Gls.) die Matrizendarstellung zu verwenden. Das Gleichungssytem
(2.44) kann z.B. in Matrizenschreibweise so dargestellt werden
R2
R3
R4
=
3, 2 −2 0
−2 4 −20 −2 3, 2
Y2
Y3
Y4
−
−2
−2−2
=
0
00
(2.45)
Symbolisch kann das Gls. so geschrieben werden
{R} = [K] {Y } − {F} = {0} , (2.46)
wobei {R} fur den Restevektor steht,
[K] =
3, 2 −2 0−2 4 −2
0 −2 3, 2
(2.47)
die globale Steifigkeitsmatrix genannt wird und
{F} =
−2−2
−2
(2.48)
28 KAPITEL 2. EINDIMENSIONALES LINEARES ELEMENT
den globalen Kraftevektor darstellt.
2.5 Elementmatrizen: Galerkin-Formulierung
In diesem Abschnitt wollen wir sehen welche Elementbeitrage in dasendgultigen Gls. einfließen und wo sie plaziert werden. Dabei nutzen
wir die Matrizenschreibweise und definieren eine Elementsteifigkeitsmatrix
und einen Elementkraftevektor fur die Dgl. aus dem vorigen Abschnitt.Fur die Entwicklung der Elementmatrizen sind stets diese drei Punkte zu
beachten:
• Restgleichungen werden immer in numerischer Reihenfolge geordnet:
R1, R2, · · · , RP−1, RP , wobei P s fur Knoten-Nr. stehen.
• Innerhalb einer Gl. werden die Knotengroßen,Φ1,Φ2, · · · ,ΦP , in nu-merischer Reihenfolge geordnet.
• Fur jeden Knoten wird eine Gl. aufgestellt.
Die Randbedingungen werden erst dann eingearbeitet, wenn alle Gln. ent-wickelt worden sind.
2.5.1 Elementmatrizen
Wir fangen mit der Definition eines Spaltenvektors {R} an. Jede Kompo-
nente von {R} reprasentiert eine Restgleichung. Dieser Vektor ist
{R} =
R1
R2...RP−1
RP
, (2.49)
Darin reprasentiert Rβ die Restgleichung fur den Knoten β. Ferner wird
die Restgleichung jedes Knoten in Elementbeitrage aufgeteilt; z.B. R(e)β
stellt den Beitrag vom Element (e) zur Restgleichung des Knoten β dar.
Aus Knotengesichtspunkt besteht die Knotenrestgleichung aus Beitrage
2.5. ELEMENTMATRIZEN: GALERKIN-FORMULIERUNG 29
der Nachbarelemente (s. Abschnitt 2.4.2 Gl. 2.19). Aus Elementengesichts-
punkt leiste jedes Element Beitrage zu den Restgleichungen seines eigenenAnfangs- und Endknoten. Fur ein beliebiges Element mit den Knoten i
und j finden wir den Beitrag zum Knoten i als
R(e)i = −
∫ Xj
Xi
Ni(x)
(
Dd2φ
dx2+Q
)
dx (2.50)
und den Beitrag zum Knoten j als
R(e)j = −
∫ Xj
Xi
Nj(x)
(
Dd2φ
dx2+Q
)
dx . (2.51)
Diese Integrale haben wir bereits in Abschnitt 2.4.3 ausgewertet: Fur s = i
und t = j ist R(e)i aquivalent mit R
(e+1)s . Fur r = i und s = j ist R
(e)j
aquivalent mit R(e)s . Unter Berucksichtigung von Gln. (2.30) und (2.36)
erhalten wir
R(e)i = D
dφ
dx
∣
∣
∣
∣
x=Xi
+D
L(Φi − Φj) −
QL
2(2.52)
R(e)j = −D
dφ
dx
∣
∣
∣
∣
x=Xj
+D
L(−Φi + Φj) −
QL
2(2.53)
Die Gln. (2.52) und (2.53) konnen in Matrizenschreibweise so
{
R(e)i
R(e)j
}
=
{
I(e)i
I(e)j
}
+D
L
[
1 −1
−1 1
]{
Φi
Φj
}
−QL
2
{
1
1
}
(2.54)
oder so geschrieben werden
{
R(e)}
={
I(e)}
+[
k(e)]{
Φ(e)}
−{
f (e)}
, (2.55)
wobei{
R(e)}
den Beitrag vom Element (e) zu dem endgultigen Glei-
chungssytem darstellt. Dieser Beitrag beinhaltet eine Elementsteifigkeits-matrix
[
k(e)]
und einen Elementkraftevektor{
f (e)}
. Andere Vektoren dar-
in sind{
Φ(e)}
=
{
Φi
Φj
}
,
30 KAPITEL 2. EINDIMENSIONALES LINEARES ELEMENT
der der Spaltenvektor der Knotengroßen ist und
{
I(e)}
=
{
I(e)i
I(e)j
}
=
D dφdx
∣
∣
∣
x=Xi
−D dφdx
∣
∣
∣
x=Xj
, (2.56)
der der Beitrag des Elements zu der interelementaren Forderung wieder-
gibt. Den Term (2.56) schließen wir aus unserer Diskussion aus, es sei dennes gibt am Knoten 1 oder P Randbedingung mit Ableitung von φ.
Elementbezogene Großen sind hier die wichtigen Ergebnisse und sie sind:Die Elementsteifigkeitsmatrix
[
k(e)]
=D
L
[
1 −1
−1 1
]
(2.57)
und der Elementkraftevektor
{
f (e)}
=QL
2
{
11
}
. (2.58)
Der Vektor R reprasentiert ein Gleichungssytem, das symbolisch wie folgt
geschrieben wird
R = [K]{Φ} − {F} = {0} (2.59)
Die Gl. (2.54) zeigt, dass die Koeffizienten in der ersten Reihe von[
k(e)]
und{
f (e)}
sich in der i-ten Reihe von [K] und {F} befinden, weil die i-te
Reihe mit{
R(e)i
}
assoziiert ist. Ahnlicherweise die Koeffizienten in der
zweiten Reihe von[
k(e)]
und{
f (e)}
befinden sich in der j-ten Reihe von
[K] und {F}, weil diese Reihe mit{
R(e)j
}
assiziiert ist. Die Koeffizienten
von[
k(e)]
befinden sich in der i-ten und j-ten Spalte von [K], weil die
Koeffizienten in der ersten Spalte Φi multiplizieren und die in der zweitenSpalte Φj multiplizieren.
2.5.2 Direkte Steifigkeitsmethode
Dieser Namen wird einer Prozedur gegeben, die Elementmatrizen direkt in
das endgultige Gleichungssystem einarbeitet. Fur ein bestimmtes Element
2.5. ELEMENTMATRIZEN: GALERKIN-FORMULIERUNG 31
werden die numerische Werte von i und j uber die Spalten von[
k(e)]
und
auf der Seite von[
k(e)]
und{
f (e)}
entlang der Zeilen geschrieben.
[
k(e)]
=
i j[
k11
k21
k12
k22
]
ij
{
f (e)}
=
{
f1
f2
}
ij
(2.60)
Um diesen Schritt zu erlautern, nehmen wir einen hypothetischen Satz
von Matrizen[
k(e)]
=
[
4 65 7
]
,{
f (e)}
=
{
89
}
fur ein lineares Element zwischen den Knoten 2 und 3 (i = 2, j = 3) an.Nach direkter Steifigkeitsmethode schreiben wir
[
k(e)]
=
2 3[
45
67
]
23
{
f (e)}
=
{
89
}
23.
Damit wird die Position jeder Koeffizient in [K] und F wie folgt bestimmt
4 addiert auf K22 7 addiert auf K33
6 addiert auf K23 8 addiert auf F2
5 addiert auf K32 9 addiert auf F3
Hier betonen wir auf den Begriff addiert auf, weil es zu K22, K23, K32,K33, F2 und F3 Beitrage von anderen Elementen geben werden kann.
2.5.3 Analyse eines Biegebalkens
Hier behandeln wir die gleiche Aufgabe wie im Abschnitt 2.4.4 mit der
direkten Steifigkeitsmethode.Die Elementsteifigkeitsmatrix ist durch Gl. (2.57)und der Elementkrafte-
vektor durch die Gl. (2.58) gegeben. die Elementdaten sind in dieser Ta-belle zusammengefasst
e i j DL
QL2
1 1 2 1,2e8 -1,0e8
2 2 3 2,0e8 -1,0e83 3 4 2,0e8 -1,0e8
4 4 5 1,2e8 -1,0e8
32 KAPITEL 2. EINDIMENSIONALES LINEARES ELEMENT
Setzen wir DL und QL
2 in Gln (2.57) und (2.58) ein, erhalten wir
[
k(1)]
= 108
1 2[
1, 2−1, 2
−1, 21, 2
]
12
,{
f (1)}
= −108
{
11
}
12
[
k(2)]
= 108
2 3[
2
−2
−2
2
]
2
3,
{
f (1)}
= −108
{
1
1
}
2
3
[
k(3)]
= 108
3 4[
2−2
−22
]
34
,{
f (1)}
= −108
{
11
}
34
[
k(3)]
= 108
4 5[
1, 2
−1, 2
−1, 2
1, 2
]
4
5,
{
f (1)}
= −108
{
1
1
}
4
5
Die globalen Matrizen [K] und {F} werden zuerst mit Nullen belegt,anschließend werden die Komponenten der Elementmatrizen auf ihre zu-
gehorigen Stellen addiert und zwar Element fur Element. Das Addierenvom Element 1 ergibt
[K] = 108
1, 2 −1, 2 0 0 0−1, 2 1, 2 0 0 0
0 0 0 0 00 0 0 0 0
0 0 0 0 0
, {F} = −108
11
00
0
Das Addiren vom Element 2 ergibt
[K] = 108
1, 2 −1, 2 0 0 0−1, 2 3, 2 −2 0 0
0 −2 2 0 00 0 0 0 0
0 0 0 0 0
, {F} = −108
12
10
0
2.5. ELEMENTMATRIZEN: GALERKIN-FORMULIERUNG 33
Das Addiren vom Element 3 ergibt
[K] = 108
1, 2 −1, 2 0 0 0−1, 2 3, 2 −2 0 0
0 −2 4 −2 0
0 0 −2 2 00 0 0 0 0
, {F} = −108
122
10
Das Addieren vom Elementt 4 schließt den Prozess ab und ergibt dasGleichungssystem
1, 2 −1, 2 0 0 0−1, 2 3, 2 −2 0 0
0 −2 4 −2 0
0 0 −2 3, 2 −1, 20 0 0 −1, 2 1, 2
Y1
Y2
Y3
Y4
Y5
−
−1−2−2
−2−1
=
000
00
(2.61)
Das Endresultat ist eine 5×5 Steifigkeitsmatrix [K] und ein 5×1 Spalten-
vektor {F}. Die Anzahl der Gln. in diesem Gleichungssytem unterscheidetsich vom (2.44), weil die Randbedingungen noch nicht eingearbeitet wor-
den sind. Die Berucksichtigung der Randbedingungen bedeutet die Eli-minierung der ersten und der funften Gl. in (2.61), weil keine Gleichung
fur bekannte Knotengroßen, Y1 und Y5, aufzustellen ist. Dies fuhrt zu dreiGleichungen mit funf Knotengroßen:
R2 = −1, 2Y1 + 3, 2Y2 − 2, 0Y3 + 2 = 0
R3 = −2, 0Y2 + 4, 0Y3 − 2, 0Y4 + 2 = 0
R4 = −2, 0Y3 + 3, 2Y4 − 1, 2Y5 + 2 = 0 . (2.62)
Das Einsetzen der bekannten Werte von Y1 und Y5 fuhrt zu einem Systemvon drei Gln. mit drei unbekannten Knotengroßen
R2 = 3, 2Y2 − 2, 0Y3 + 2 = 0
R3 = −2, 0Y2 + 4, 0Y3 − 2, 0Y4 + 2 =
R4 = −2, 0Y3 + 3, 2Y4 + 2 = 0 , (2.63)
die mit Gln. (2.44) identisch sind.
34 KAPITEL 2. EINDIMENSIONALES LINEARES ELEMENT
2.5.4 Eigenschaften der globalen Steifigkeitsmatrix
Fur Strukturprobleme und fur selbst adjugierte Dgln. ist die Steifigkeits-
matrix [K] immer symmetrisch und positiv definiert. Jeder Koeffizientauf der Hauptdiagonale, Kii, ist im Vergleich zu den anderen auf der sel-
ben Reihe relativ groß. FE-Gleichungen werden oft mit der Methode derGaus’schen Elimination oder einer ihren Erweiterungen gelost, weil das
Gleichungssytem nicht immer diagonal dominant ist; d.h. Kii kann klei-ner als die Summe der restlichen Komponente der selben Reihe, i, sein.Die relativ großen Diagonalkomponenten ermoglichen, dass die Gaus’sche
Elimination ohne Reihenwechsel (Pivoting) erfolgen kann. Aus rechen-technischer Sicht ist dies wichtig, weil damit nur die Koeffizienten zu spei-
chern sind, die vom Null unterschieden sind. Die Symmetrieeingenschaftist ebenso wichtig, weil sie die Speicherung der Koeffizienten unter der
Hauptdiagonale uberflussig macht.
Die Symmetrieeigenschaft und die positive Definiertheit ergeben sich ausder mathematischen Formulierung. Die globale Steifigkeitsmatrix hat aberauch eine andere Eigenschaft, namlich die Bandformigkeit; d.h. die vom
Null unterschiedenen Koeffizienten dieser Matrix befinden sich innerhalbeines Bandes, dessen Breite vom Elementnetz und der Numerierung der
Knoten abhangt. Wahrend es innerhalb des Bandes Koeffizienten gleichNull geben kann, sind alle Koeffizienten außerhalb des Bandes gleich Null.
Die Bandformigkeit ergibt sich aus der Konstruktion des Gleichungssy-tems, die hier nicht naher erlautert wird. Je kleiner die Bandbreite, um-
so einfacher und kosteneffektiver wird die Losung des Gleichungssytemssein. Um die Bandbreite klein zu halten, ist es notwendig die Kontennu-merierung von Elementen so zu gestalten, dass die Differenz zwischen den
Knotennummern jedes Elementes moglichst gering bleibt. Wenn BW (e)
die Differenz zwischen den Knotennummern eines Elementes im FE-Netz
darstellt, lasst sich die Bandbreite aus folgender Gleichung berechnen
NBW = max[
BW (e)]
+ 1 . (2.64)
Kapitel 3
Zweidimensionale Elemente
Ein primarer Vorteil der FE-Methode ist die Einfachheit ihrer Verallge-
meinerung fur die Losung zweidimensionaler Probleme, die aus vielen un-terschiedlichen Materialien bestehen und nicht regulare Rander besitzen.
Es gibt eine Vielzahl von allgemeinen FE-Programmen fur die Losungzweidimensionaler Probleme. Alle diese Programme verwenden Dreiecks-
und Rechtseckelemente oder die generalisierten Formen dieser Elemen-te. Wir beginnen mit dem linearen Dreieckselement und dem bilinearenRechteckselement.
3.1 Zweidimensionale Netze
Das lineare Dreieckselement hat gerade Seiten und einen Knoten an jeder
Ecke. Die Interpolationsfunktion fur eine skalare Große im Element lautet
φ = α1 + α2x+ α3y . (3.1)
Diese Gl. ist ein vollstandiges lineares Polynom, weil sie einen konstantenTerm und alle mogliche lineare Terme, namlich x und y besitzt. Folglich
kann das Dreieckselement jede Orientierung haben und die Kontinuitats-bedingungen fur die benachbarten Elemente erfulen.
Das bilineare Rechteckselement hat ebenso gerade Seiten und einen Kno-
ten an jeder Ecke. Die Interpolationsfunktion fur eine skalare Große imElement ist
φ = C1 + C2x+ C3y + C4xy . (3.2)
35
36 KAPITEL 3. ZWEIDIMENSIONALE ELEMENTE
(a)
i
j
ky
xi
y
xj
km
(b)
Abb. 3.1: (a) Das lineare Dreieckselement, (b) Das bilineare Rechteckselement.
Diese Gl. hat nur eine aus drei moglichen Termen der zweiten Ordnung,
xy. Das Rechteck kann nicht beliebig orientiert sein, weil die quadratischenTermen, x2 und y2, fehlen. Daher mussen die Seiten des Rechteckelementes
parallel zu den Koordinatenachsen sein. Ein Netz von Rechteckelementenist leicht zu konstruieren. Alle Elemente in einer Reihe sollen die gleiche
Hohe und alle Elemente in einer Spalte die gleiche Breite haben. Durchpassende Kombination von Dreiecks- und Rechteckselmeneten kann manauch irregulare Gebiete vernetzen. Debei teilt man zuerst das Gebiet in
groben Teilgebieten, die man sichtbar leichter mit Dreiecks- oder Recht-eckselementen diskretisieren kann. Die Knoten an Grenzen der Teilgebie-
ten mussen identisch sein (in Anzahl, Nummer und Position), damit dieKontinuitat in φ gewahrleistet ist. Die gekrummten Gebietsrander sind
durch eine Reihe von Geraden zu ersetzen.
3.2 Das lineare Dreieckselement
Das lineare Dreieckselement in Bild (3.2) hat gerade Seiten und drei Kno-ten, die an Ecken liegen. Zur konsistenten Bezeichnung der Knoten und
Seiten fangen wir mit einem beliebigen Knoten an und bezeichnen ihnmit i, dann schreiten wir gegen Uhrzeigersinn fort. Die Knotenkoordina-
ten bezeichnen wir mit (Xi, Yi), (Xj, Yj) und (Xk, Yk). Die Knotengroßensind: Φi, Φj und Φk. Die Ansatzfunktion fur dieses Element wurde durch
Gl. (3.1) gegeben. Damit gilt fur die Knotengroßen:
3.2. DAS LINEARE DREIECKSELEMENT 37
(X ,Y )j j
(X ,Y )k k
i i(X ,Y )
φ=α +α +α21 x 3y
x
y
φ
Φj
Φk
Φi
j
k
i
Abb. 3.2: Parameter des linearen Dreieckselementes.
φ = Φi an x = Xi, y = Yi
φ = Φj an x = Xj, y = Yj
φ = Φk an x = Xk, y = Yk
Durch Einsetzen dieser Werte in Gl. (3.1) erhalten wir das Gleichungssy-tem
Φi = α1 + α2Xi + α3Yi
Φj = α1 + α2Xj + α3Yj (3.3)
Φk = α1 + α2Xk + α3Yk
Die Auflosung nach unbekannten Koeffizienten ergibt
α1 =1
2A[(XjYk −XkYj)Φi + (XkYi −XiYk)Φj + (XiYj −XjYi)Φk]
α2 =1
2A[(Yj − Yk)Φi + (Yk − Yi)Φj + (Yi − Yj)Φk]
α3 =1
2A[(Xk −Xj)Φi + (Xi −Xk)Φj + (Xj −Xi)Φk] ,
38 KAPITEL 3. ZWEIDIMENSIONALE ELEMENTE
wobei A der Flacheninhalt des Dreiecks ist und das zweifache davon der
Determiniante der Koeffizientenmatrix gleich ist:∣
∣
∣
∣
∣
∣
1 Xi Yi
1 Xj Yj
1 Xk Yk
∣
∣
∣
∣
∣
∣
= 2A (3.4)
Das Einsetzen von α1, α2 und α3 in Gl. (3.1) ergibt eine Gleichung fur φ,die man in Termen von Knotengroßen, Φi, Φj und Φk, und Formfunktio-
nen, Ni, Nj und Nk, ausdrucken kann:
φ = NiΦi +NjΦj +NkΦk (3.5)
Darin sind
Ni =1
2A[ai + bix+ ciy] (3.6)
Nj =1
2A[aj + bjx+ cjy] (3.7)
Nk =1
2A[ak + bkx+ cky] (3.8)
und
ai = XjYk −XkYj , bi = Yj − Yk und ci = Xk −Xj
aj = XkYi −XiYk , bj = Yk − Yi und cj = Xi −Xk
ak = XiYj −XjYi , bk = Yi − Yj und ck = Xj −Xi
Die skalare Große φ wurde durch Formfunktionen, die in x und y linearsind, mit den Knotengroßen verknopft. Dies bedeutet, dass die Gradienten∂φ/∂x und ∂φ/∂y innerhalb des Elementes konstant sind. Zum Beispiel
∂φ
∂x=∂Ni
∂xΦi +
∂Nj
∂xΦj +
∂Nk
∂xΦk , (3.9)
und mit∂Nβ
∂x=
bβ2A
β = i, j, k
erhalt man∂φ
∂x=
1
2A[biΦi + bjΦj + bkΦk] . (3.10)
Die Ableitung ∂φ/∂x hat einen konstanten Wert, weil bi, bj und bk Kon-
stanten und Φi, Φj und Φk unabhangig von Raumkoordinaten sind. Ein
3.3. DAS BILINEARE RECHTECKSELEMENT 39
konstanter Gradient innerhalb des Elementes bedeutet, dass viele klei-
ne Elemente notwendig sind, um starke Anderungen in φ relativ genauapproximieren zu konnen.
Beispiel 3.1
3.3 Das bilineare Rechteckselement
Das bilineare Rechteckselement hat eine Lange von 2b und eine Hohe von
2a. Sein Knoten sind durch i, j, k und m bezeichnet, wobei i immer denKnoten unten links bezeichnet und die fortlaufende Bezeichnung gegen den
Uhrzeigersinn lauft. Das Element und die wichtigen Koordinatensystemesind in Bild (3.3) dargestellt. Die Ansatzfunktion in Gl. (3.2) ist hier in
Termen der lokalen Koordinaten, s und t, umgeschrieben:
φ = C1 + C2s+ C3t+ C4st (3.11)
Es gibt mindestens zwei weitere Formen der Umschreibung, die den Termst entweder durch s2 oder durch t2 ersetzen. Die Gl. (3.11) wird hier be-
nutzt, weil φ entlang jeder Linie mit konstantem t in s linear ist und ent-lang jeder Linie mit konstantem s in t linear ist. Wegen dieser Eigenschaf-ten bezeichnet man das Element als bilinear. Die Ansatzfunktion (3.11)
Φi
Φm
Φj
Φk
φ =C +C +C +C 1 2 3s t 4 st
x
y
q
φ
2b
2a
ij
k
mt
s
r
Abb. 3.3: Parameter des bilinearen Rechteckselementes.
40 KAPITEL 3. ZWEIDIMENSIONALE ELEMENTE
wurde bezogen auf ein lokales Koordinatensytem geschrieben, dessen Ur-
sprung auf dem Knoten i liegt, weil die Auswertung der Formfunktionenin diesem Referenzsystem einfacher ist. Ein anderes wohl bekanntes Be-
zugssytem ist das qr-System, dessen Ursprung auf dem Mittelpunkt desElementes liegt und seine Achsen parallel zu den Elementseiten laufen.Die Koeffizienten C1 bis C4 in Gl. (3.11) werden durch Einsetzen der Kno-
tenwerten von φ und Knotenkoordinaten (in st-System) gerechnet. Manerhalt vorerst die folgenden vier Gln.:
Φi = C1
Φj = C1 + (2b)C2
Φk = C1 + (2b)C2 + (2a)C3 + (4ab)C4
Φm = C1 + (2a)C3
Die Losung nach unbekannten Koeffizienten ergibt
C1 = Φi
C2 =1
2b(Φj − Φi)
C3 =1
2a(Φm − Φi) (3.12)
C4 =1
4ab(Φi − Φj + Φk − Φm)
Durch Einsetzen der Koeffizienten in Gl. (3.11) und Umschreiben erhalten
wirφ = NiΦi +NjΦj +NkΦk +NmΦm . (3.13)
Darin sind
Ni =(
1 −s
2b
)
(
1 −t
2a
)
Nj =s
2b
(
1 −t
2a
)
Nk =st
4ab(3.14)
Ni =t
2a
(
1 −s
2b
)
Die Formfunktionen fur das bilineare Rechteckselement haben ahnlicheEigenschaften wie die vom Dreieckselement. Jede Formfunktion lauft ent-
lang zweier Elementkanten zwischen ihrem zugehorigen Knoten und zwei
3.4. EINE KONTINUIERLICHE STUCKWEISE GLATTE FUNKTION 41
benachbarten Knoten linear. Zum Beispiel lauft Ni entlang der Seiten ij
und im linear. Jede Formfunktion ist an den Seiten, die nicht mit ihremKnoten verbunden sind, gleich Null, d. h. Ni ist enlang der Seiten jk und
km gleich Null. Der lineare Verlauf von φ an Seiten des Rechteckselementesund an Seiten des Dreieckselementes bedeutet, dass diese Elemente zuein-ander kompatible sind und nebeneinander verwendet werden konnen.
Die Transformationsgleichungen zwischen qr-Koordinatensystem und st-
Koordinatensystem sind
s = b+ q und t = a+ r . (3.15)
Einsetzen dieser Transformationen in Gln. (3.14) ergibt
Ni =1
4
(
1 −q
b
)(
1 −r
a
)
Nj =1
4
(
1 +q
b
)(
1 −r
a
)
Nk =1
4
(
1 +q
b
)(
1 +r
a
)
(3.16)
Ni =1
4
(
1 −q
b
)(
1 +r
a
)
Diese Formfunktionen sind sehr nutzlich, weil sie zu einem naturlichen
Koordinatensystem fuhren, das dem Rechteckselement erlaub die Formeiner allgemeinen Viereckelement anzunehmen.Eine Kontur- oder Profillinie ist in einem Rechteckselement im Allgemei-
nen gekrummt. Die zwei Schnittstellen einer Konturlinie mit den Ele-mentkanten konnen aus linearen Interpolationen gerechnet werden. die
einfachste Methode zur Berechnung eines dritten Punktes ist, s oder t inden Formfunktionen gleich Null setzen und dann die Gl. (3.13) fur die
andere Koordinate losen.
Beispiel 3.2
3.4 Eine kontinuierliche stuckweise glatte Funktion
Die in Gl. (3.5) und (3.13) definierten Interpolationsfunktionen kann man
fur jedes Dreiecks- bzw. Rechteckselement verwendet, wenn man die ent-sprechenden Knotenwerte von i, j und k bzw. von i, j, k und m ein-
setzt. Jeder Knoten eines Dreieckselementes kann der Knoten i sein.
42 KAPITEL 3. ZWEIDIMENSIONALE ELEMENTE
1
2
3
4
5
6
7
8
(1)
(2)
(3)
(4)
*
*
*
*
Abb. 3.4: Ein FE-Netz ausDreiecks-und Rechteckselemen-ten. * zeigt den Knoten i vonjedem Element
Bei einem Rechteckselement liegt der Knoten i stets am Ursprung desst-Koordinatensystems. Die folgende Tabelle stellt die Element-Knoten-Daten fur das FE-Netz in Bild (3.4) dar, in dem der Knoten i mit einem
∗ bezeichnet wurde:
e i j k m
1 1 4 5 22 2 5 6 3
3 3 6 74 8 3 7
Die Interpolationsfunktion fur das Element 1 ist
φ(1) = N(1)1 Φ1 +N
(1)4 Φ4 +N
(1)5 Φ5 +N
(1)2 Φ2 (3.17)
Hier sind die Knotennummer eines Elements nicht mehr in Reihenfolge,und das ist in zweidimensionalem Fall selbstverstandlich. Die Formfunk-
tionen in Gl. (3.5) sind Funktionen der globalen Koordinaten, in demSinne, dass es gilt:
2b = Xj −Xi = X4 −X1
und2a = Ym − Yi = Y2 − Y1
Die Interpolationsfunktion fur das Element 4 ist
φ(4) = N(4)8 Φ8 +N
(4)3 Φ3 +N
(4)7 Φ7 (3.18)
3.4. EINE KONTINUIERLICHE STUCKWEISE GLATTE FUNKTION 43
Die Formfunktionen in (3.18) sind Funktionen der globalen Koordinaten
und die Bezeichnungen i, j, und k machen sofort deutlich, welche Koor-dinaten einzusetzen sind. Betrachten wir zum Beispiel N
(4)8 . Aus Gl. (3.6)
erhalten wir
N(4)8 =
1
2A(a
(4)8 + b
(4)8 x + c
(4)8 y) ,
wobei mit j = 3 und k = 7 gilt
a(4)8 = X3Y7 −X7Y3
b(4)8 = Y3 − Y7
c(4)8 = X7 −X3
Die Flache A bezieht sich auf das Element 4. Formfunktionen wieN(4)8 sind
kontinuierlich und in x und y linear. Dies Macht die Interpolationsfunktionin (3.17) oder (3.18) zu einer kontinuierlichen stuckweise glatten Funktion.
Kapitel 4
Koordinatensysteme
Alle FE-Losungen sind an Auswertung von Integralen angewiesen, die
nicht immer leicht zu machen ist. Manche Integrale konnen gar nicht ana-lytisch ausgewertet werden und mussen mit numerischen Techniken be-handelt werden. Die Schwierigkeiten einer Integralauswertung konnen oft
durch Variablenanderung gemindert werden. Dies schließt die Umschrei-bung des Integrals in einem neuen Koordinatensystem ein. Hier lernen
wir ein paar Koordinatensysteme kennen, die die Auswertung von FE-Integralen erleichtern konnen. Betrachtet werden lokale Koordinatensy-
steme und naturliche Koordinatensysteme fur das lineare eindimensionaleElement, das lineare Dreieckselement und das bilineare Rechteckselement.
4.1 Lokale Koordinatensysteme
Fur ein lineares Element wurden im Abschnitt 2 folgende Formfunktionenabgeleiete
Ni =Xj − x
Lund Nj =
x−Xi
L, (4.1)
wobei der Ursprung des KS irgendwo auf der linken Seite des Knoten i
gesetzt wurde. Diese Formfunktionen sind fur alle lineare 1D-Elementeunabhangig von ihrer Lage im Raum gultig. Ihrer Nachteil zeigt sich erst,wenn man Integrale auswerten will, die Produkten solcher Formfunktionen
beinhalten, wie z.B.∫ Xj
Xi
Ni(x)Nj(x)dx oder
∫ Xj
Xi
N2i (x)dx (4.2)
45
46 KAPITEL 4. KOORDINATENSYSTEME
i j
s
x L
(a)
i jx
L
L / 2
q
(b)
Abb. 4.1: Lokale Koordinatensysteme fur das lineare eindimensionale Element
Solche Integrale treten bei Feldproblemen und Festkorpermachanik auf.
Legt man den Ursprung des KS auf dem Element fest, kann man neueFormfunktionen ableiten, die die Integrale in (4.2) vereinfachen. Diese
Art von KS wird als lokales Koordinatensystem, LKS, bezeichnet. Zweider bakanntesten LKS fur das lieare 1D-Element haben ihren Ursprung
auf dem Knoten i bzw. auf dem Mittelpunkt des Elementes (s. Bild 4.1).Die Formfunktionen fur das LKS mit Ursprung auf dem Knoten i erhaltman durch den Austausch von x mit Xi +s in Gln. (4.1). Die Substitution
ergibt:
Ni(s) =Xj − x
L=Xj − (Xi + s)
L= 1 −
s
L(4.3)
und
Nj(s) =x−Xi
L=
(Xi + s) −Xi
L=s
L(4.4)
Die Koordinate s befindet sich im Bereich [0, L. Man sieht, dass jedeFormfunktion weiterhin an eigenem Knoten den Wert 1 und am anderen
den Wert Null hat.
Die Formfunktionen fur das LKS mit Ursprung auf dem Mittelpunkt desElementes erhalt man dadurch, dass man x in Gln. (4.1) durchXi+(L/2)+q ersetzt. Die Formfunktionen bezogen auf dieses LKS sind:
Ni(q) =
(
1
2−q
L
)
und Nj(q) =
(
1
2+q
L
)
(4.5)
Die Koordinate q befindet sich im Bereich [−L/2, L/2].
Die Formfunktionen in (4.3), (4.4) und (4.5) sind erst dann nutzlich, wenn
eine Variablenanderung im Integral vollzogen wird. Die Formel fur Varia-
4.1. LOKALE KOORDINATENSYSTEME 47
blenanderung ist nach Integralrechnung wie folgt
∫ b
a
f(x)dx =
∫ p2
p1
f(g(p))
[
d(g(p))
dp
]
dp , (4.6)
wobei p die neue Koordinatenvariable und g(p) die Beziehung zwischen
x und p darstellt: x = g(p). Die Interpretation der Formel (4.6) bezogenauf die LKS in Bild 4.1 lauft wie im Folgenden. Fur die Koordinate s gilt:x = Xi + s, und hier bekommt man
∫ Xj
Xi
f(x)dx =
∫ s2
s1
f(Xi + s)d(Xi + s)
dsds =
∫ L
0
h(s)ds , (4.7)
wobei h(s) die f(x) in Termen von s wiedergibt. Zur Umrechnung desIntegralbereiches wurde x in x = Xi + s durch Xi bzw. Xj ersetzt und dieGleichung fur s gelost.
Fur die Koordinate q gilt: x = Xi + L/2 + q, und hier erhalt man
∫ Xj
Xi
f(x)dx =
∫ q2
q1
r(q)d(Xi + L/2 + q)
dqdq =
∫ L/2
−L/2
r(q)dq , (4.8)
wobei, r(q) die f(x) in Termen von q beschreibt.Der Vorteil von Gln. (4.7) und (4.8) zeigt sich bei der Auswertung von
Integralen wie∫ Xj
Xi
N2i dx .
In diesem Fall bekommen wir mit der Koordinate s∫ Xj
Xi
N2i dx =
∫ L
0
N2i (s)ds =
∫ L
0
(
1 −s
L
)2
ds =L
3.
Und mit der Koordinate q erhalten wir
∫ Xj
Xi
N2i dx =
∫ L/2
−L/2
N2i (q)dq =
∫ L/2
−L/2
(
1
2−q
L
)2
dq =L
3.
Versuchen Sie bitte selbst die Auswertung des oberen Integrals ohne Va-
riablenanderung.
48 KAPITEL 4. KOORDINATENSYSTEME
4.2 Naturliche Koordinatensysteme
Ein naturliches KS ist ein lokales KS, das erlaubt einen Punkt innerhalbdes Elements durch eine dimensionslose Zahl festzulegen, derer absoluterWert nie großer als 1 wird. Die lokale Koordinatensysteme s und q konnen
zur naturlichen KS konvertiert werden.
Wir fangen mit der q-Koordinate in Bild 4.1 an und bilden das Langen-verhaltnis q/(L/2) = 2q/L = ξ. Die Koordinate ξ variiert von -1 bis +1 (s.
Bild 4.2a). Die Formfunktionen in Gln. (4.5) konnen mittels Substitutionvon q durch q = ξL/2 in Termen von ξ geschrieben werden:
Ni(ξ) =1
2(1 − ξ) und Nj(ξ) =
1
2(1 + ξ) (4.9)
Die Variablenanderung bei Integration ergibt:
∫ L/2
−L/2
r(q)dq =
∫ ξ2
ξ1
g(ξ)d(ξL/2)
dξ=L
2
∫ 1
−1
g(ξ)dξ , (4.10)
wobei g(ξ) die r(q) in Termen von ξ wiedergibt. Der Vorteil von Koordi-natenvariable ξ liegt im daraus resultierenden Integrationsbereich [−1, 1].
Die meisten Rechenprogramme nutzen numerische Integrationstechnikenzur Auswertung von Elementmatrizen. Ein numerisches Integrationsver-
fahrern, das in FE-Programmen verwendet wird, ist das Gauß-Legendre-Verfahren, wobei Stutzpunkte und Gewichtsfaktoren im Interval [−1, 1]
eingesetzt werden.
Ein zweites naturliches KS beinhaltet ein Paar von Langenverhaltnissen(s. Bild 4.2b). Wenn s den Abstand vom Knoten i angibt, dann werden
i jL / 2
L
ξ = 1ξ = −1 ξ
(a)
l 2 l1
i js
L
(b)
Abb. 4.2: Naturliche Koordinatensysteme fur das lineare eindimensionale Element
4.2. NATURLICHE KOORDINATENSYSTEME 49
ℓ1 und ℓ2 als die Langenverhaltnisse wie folgt definiert
ℓ1 =L− s
Lund ℓ2 =
s
L(4.11)
Die Koordinaten ℓ1 und ℓ2 sind von einander abhangig, denn es gilt
ℓ1 + ℓ2 = 1 (4.12)
Die wichtigste Eigenschaft von Gln. (4.11) und (4.12) ist, dass ℓ1 undℓ2 mit den Formfunktionen in (4.3) und (4.4) identisch sind. Der Vorteil
dieser Koordinaten zeigt sich bei der Auswertung von Integralen der Art
∫ L
0
Nai (s)N b
j (s)ds , (4.13)
die das Produkt von Formfunktionen beinhalten. Die Regel der Varia-blenanderung und die Beziehungen Ni(s) = ℓ1, Nj(s) = ℓ2, s = Lℓ2 und
ds/dℓ2 = L ergeben
∫ L
0
Nai (s)N b
i (s)ds =
∫ 1
0
ℓa1ℓb2Ldℓ2 (4.14)
Das Integral auf der rechten Seite von (4.14) lasst sich mit Hilfe von Gl.(4.12) wie folgt umschreiben
L
∫ 1
0
(1 − ℓ2)aℓb2dℓ2 (4.15)
Diese Integral ist von der gleichen Form wie
∫ 1
0
tz−1(1 − t)w−1dt =Γ(z)Γ(w)
Γ(z + w, (4.16)
wobei Γ(n+ 1) = n! ist. Daher
L
∫ 1
0
ℓa1ℓb2dℓ2 = L
Γ(a+ 1)Γ(b+ 1)
Γ(a+ b+ 1 + 1
= La!b!
(a+ b+ 1)!(4.17)
50 KAPITEL 4. KOORDINATENSYSTEME
Koordinaten- Koordinaten- Integrations-system variable Formfunktionen bereich
Global x Ni =Xj−x
L, Nj = x−Xi
L[Xi, Xj]
Lokal s Ni =1 − s
L, Nj = s
L[0, L]
Lokal q Ni =(
12− q
L
)
, Nj =(
12
+ q
L
) [
−L
2, L
2
]
Naturlich ξ Ni =12(1 − ξ) , Nj = 1
2(1 + ξ) [−1, 1]
Naturlich ℓ2 Ni = ℓ1 , Nj = ℓ2 [0, 1]
Tabelle 4.1: Koordinatensysteme und Integrationsbereiche fur das lineare 1D-Element
Diese Gleichung besagt, dass ein relativ kompliziertes Integral mit Hilfeeiner Gleichung, die die Lange des Elements und ie Ordnungen der Pro-dukten beinhaltet, ausgewertet werden kann. Dies wird durch Auswertung
einiger Integrale gezeigt. Begonnen mit∫ Xj
Xi
N2i (x)dx =
∫ L
0
N2i (s)ds
erhalten wir von (4.12)∫ L
0
N2i (s)ds = L
∫ 1
0
ℓ21ℓa2dℓ2 = L
2!0!
(2 + 0 + 1)!=L
3.
Ein anderes Beispiel ist∫ L
0
N3i (s)N2
j (s)ds = L
∫ 1
0
ℓ31ℓ22dℓ2 = L
3!2!
(3 + 2 + 1)!=
L
60
Die Koordinatensysteme, Formfunktionen und Integrationsbereiche furdas eindimensionale lineare Element sind in Tabelle 4.2 aufgelistet.
4.3 Das Rechteckselement
Fur das Rechteckselement haben wir bereits im Abschnitt 3.2 zwei LKS,
st− und qr−KS, und daraus abgeleiteten Formfunktionen vorgestellt.
4.4. DAS DREIECKSELEMENT: FLACHENKOORDINATEN 51
η =1
=−1η
ξ=
−1
=1
ξ
i j
km
η
ξ Abb. 4.3: Ein naturliches Koor-dinatensystem fur das Rechtecks-element
Das NKS fur das Rechteckselement hat seinen Ursprung am Mittelpunkt
des Elementes (s. Bild 4.3) und seine Koordinaten sind die Langenverhalt-nisse:
ξ =q
bund η =
r
a, (4.18)
wobei q und r die lokalen Koordinaten sind (s. Bild 3.3). Somit konnen
die Formfunktionen in (3.16) leicht zu naturlichen koordinaten konvertiertwerden. Das Ergebnis lautet:
Ni =1
4(1 − ξ)(1 − η), Nj =
1
4(1 + ξ)(1 − η)
Nk =1
4(1 + ξ)(1 + η), Nm =
1
4(1 − ξ)(1 + η) (4.19)
Es ist zu beachten, dass ξ und η zwischen -1 und 1 variieren, d.h.
−1 ≤ ξ ≤ 1 und − 1 ≤ η ≤ 1
4.4 Das Dreieckselement: Flachenkoordinaten
Ein naturliches Koordinatensystem fur das Dreieckselement ist durch De-
finition von drei Langenverhalnissen L1, L2 und L3 zu erhalten (s. Bild4.4a). Jede Koordinate ist ein Verhaltnis, das zwischen dem Abstand s
von einer Seite und der Hohe h fur die gleiche Seite besteht (s. Bild 4.4b).Damit variiert jede Koordinate zwischen 0 und 1. Die Linien der konstan-
ten L1 sind in Bild 4.4c dargestellt. Diese sind alle parallel zur Bezugsseite
52 KAPITEL 4. KOORDINATENSYSTEME
L 3
L1
L 2
i
j
k
(a)
L1
s
b
B
h
k
i
j
(b)
L1
L =11
L =1/21
1L =1/4L =3/41
L =01
k
i
j
(c)
Abb. 4.4: Die drei Flachenkoordinaten fur ein Dreieckselement
jk, von der L1 gemessen wurde.Die Koordinaten L1, L2 und L3 werden Flachenkoordinaten genannt, weil
ihre Werte das Flachenverhaltnis zwischen einem Unterdreieck und demganzen Dreieckelement angeben. Betrachten wir den Punkt B im Drei-eckselement in Bild 4.5. Die Flache des ganzen Dreiecks ist A und gegeben
durch
A =bh
2,
wahrend die Flache des schraffierten Dreiecks (B, j, k) ist
A1 =bs
2(4.20)
Das Flachenverhaltnis A1/A2 ergibt
A1
A=s
h= L1 (4.21)
Die Koordinate L1 stellt das Verhaltnis der schraffierten Flache in Bild4.5 zur Gesamtflache dar. Ahnliche Gleichungen konnen fur L2 und L3
geschrieben werden:
L2 =A2
Aund L3 =
A3
A(4.22)
Da es gilt A1 +A2 + A3 = A, erhalt man
L1 + L2 + L3 = 1 (4.23)
Diese Gleichung beschreibt die Abhangigkeit der drei Koordinaten. Der
Ort eines Punktes lasst sich mit nur zwei Koordinaten festlegen.
4.4. DAS DREIECKSELEMENT: FLACHENKOORDINATEN 53
�����������������������������������
�����������������������������������
A1A2
A3
bh
k
i
j
s
B
x
y
Abb. 4.5: Aufteilung eines Drei-eckselementes in Teilflachen ent-sprechend der Flachenkoordina-ten
Die Gl. (4.21) kann in eine andere Form wiedergeben werden: Multiplika-
tion von Zahler und Nenner mit 2 ergibt
L1 =2A1
2A(4.24)
Die Determinante-Erweiterung fur 2A1 liefert:
2A1 =
∣
∣
∣
∣
∣
∣
1 x y
1 Xj Yj
1 Xk Yk
∣
∣
∣
∣
∣
∣
oder
2A1 = (XjYk −XkY j) + (Yj − Yk)x+ (Xk − Yj)y , (4.25)
wobei x und y die Koordinaten von B in Bild 4.5 sind. Einsetzen von
(4.25) in GL. (4.24) ergibt
L1 =1
2A[(XjYk −XkY j) + (Yj − Yk)x+ (Xk − Yj)y] . (4.26)
Diese Gl. ist mit der Gl. (3.6) identisch; daher gilt
L1 = Ni (4.27)
Ahnlische Verfahren fur L2 und L3 fuhren zu
L2 = Nj und L3 = Nk (4.28)
54 KAPITEL 4. KOORDINATENSYSTEME
Die Flachenkoordinaten fur das lineare Dreieckselement sind mit den
Formfunktionen identisch und konnen als zwei tauschbare Wertesatze ver-wendet werden.
Der Vorteil von Flachenkoordinaten liegt in Existenz einer Integralglei-chung, die die Auswertung von Flachenintegralen vereinfacht (Eisenbergund Malvern, 1973). Diese Integralgleichung ist mit der Gl. (4.17) ver-
wandt und lautet∫
A
La1L
b2L3cdA =
a!b!c!
(a+ b+ c+ 2)!2A (4.29)
Die Nutzung von (4.29) kann durch Auswertung vom Integral eines Pro-
duktes von Formfunktion∫
A
Ni(x, y)Nj(x, y)dA (4.30)
uber die Flache eines Dreiecks erortert werden. Die Flachenintegral ist∫
A
NiNjdA =
∫
A
L11L
12L
03dA
=1!1!0!
(1 + 1 + 0 + 2)!2A =
2A
4!=A
12
Die Flachenkoordinaten L1 und L2 konnen anstelle Ni bzw. Nj eingesetztwerden. Da Nk nicht im Produkt vorhanden war, wurde L3 mit der Potenz
0 eingesetzt. Die 0 Fakultat ist als 1 definiert.
Die Berucksichtigung von abgeleiteten Randbedingungen oder von
Flachenlasten in FE-Analyse fordert, dass ein Integral entlang einer Ele-mentkante ausgewertet wird. Diese Auswertung ist dann leicht, wenn man
weiß, wie sich die Flachenkoordinate auf dieser Kante verhalt. Betrachtewir den Punkt B auf der Seite ij (s. Bild 4.6). Die Koordinate L3 istgleich Null und L1 gleicht dem Verhaltnis der schraffierten Flache zur Ge-
samtflache. Definiert wird hier die Koordinatenvariable s, die zur Seite ijparallel ist und vom Knoten i gemessen wird. Wenn die Koordinate von
Punkt B gleich s ist und die Lange der Seite gleich b, dann erhalt man
L1 =2A1
2A=
2h(b−s)2
2bh2
=b− s
b= 1 −
s
b(4.31)
4.4. DAS DREIECKSELEMENT: FLACHENKOORDINATEN 55
����������������������������������������
����������������������������������������
A1
L 1
k
j
sb
i
B
h
x
y
Abb. 4.6: Die Flachenkoordina-ten fur einen Punkt auf einer Sei-te des Dreieckselementes
Die Flachenkoordinate L2 bekommt man als
L2 =s
b(4.32)
Die Flachenkoordinaten L1 und L2 reduzieren sich auf eindimensionale
Formfunktionen Ni(s) und Nj(s), die in Gln. (4.3) und (4.4) definiertworden sind. Wenn wir von eindimensionalen naturlichen Koordinaten ℓ1und ℓ2 gebrauch machen, die in Gl. (4.11) definiert wurden, dann werdendie Flachenkoordinaten zu
L1 = ℓ1 und L2 = ℓ2 side i→ j (4.33)
Die beziehungen fur die zwei anderen Seiten sind
L2 = ℓ1 und L3 = ℓ2 side j → k (4.34)
L3 = ℓ1 und L1 = ℓ2 side k → i (4.35)
Das Wichtige an Gln. (4.33), (4.34) und (4.35) ist, dass jedes Integral ubereine Seite eines Dreieckselementes durch ein Linienintegral in Termen von
s und ℓ2,∫
Γ
f(L1, L2, L3)dΓ =
∫ L
0
g(s)ds = L
∫ 1
0
h(ℓ2)dℓ2 , (4.36)
ersetzt und mittels Fakultatformel (4.17) ausgewertet werden kann. Der
Rand eines zweidimensionalen Elementes wurde hier durch Γ bezeichnet.
Beispiel 4.1
56 KAPITEL 4. KOORDINATENSYSTEME
4.5 Kontinuitat
Die Naherung von φ(x, y) in einem Bereich beinhaltet eine Reihe von kon-tinuierlichen stuckweise glatten Funktionen, die jeweils uber ein Element
definiert sind. Der Bedarf an Integration stellt eine Forderung an die Kon-tinuitatsordnung zwischen den Elementen. Das Integral
∫ H
0
dnφ
dxn
ist nur dann definiert, wenn φ eine Kontinuitat der Ordnung (n − 1)hat. Dies stellt sicher, dass es nur Diskontinuitat in Form von endlichen
Sprunge gibt. Diese Forderung bedeutet, dass zwischen den Elementen dieerste Ableitung der Naherungsfunktion kontinuierlich sein muss, wenn dasIntegral Termen mit zweiter Ableitung (n = 2) beinhaltet. Alle Integra-
le in diesem Skript, ausgenommen fur das Balkenelement, haben Termenmit der ersten Ableitung. Daher muss φ zwischen den Elementen konti-
nuierlich sein, aber ihre Ableitungen nicht. Kontinuitat in der Ableitungist jedoch fur das Balkenelement erforderlich.
Kontinuitat von φ ist fur das eindimensionale Element gewahrleistet,weil zwei benachbarte Elemente einen gemeinsamen Knoten haben. Kon-tinuitat von φ entlang einer gemeinsamen Seite von zwei Rechtecksele-
menten ist relativ leicht nachzuweisen. Entlang einer gemeinsamen Grenz
Φ1
Φ2
Φ4
Φ3
12
3
4(1)
(2)
φ
x
y
**
Abb. 4.7: Ein Netz von zweiDreieckselementen
4.5. KONTINUITAT 57
zwischen zwei beliebig orientierten Dreieckselementen ist es etwas schwe-
rer und wird hier behandelt.Betrachten wir die zwei benachbarten Elemente in Bild 4.7 mit Koordi-
natenursprung auf den Knoten 1. Die Knotengroßen sind Φ1, φ2, φ3 undφ4. Die Gln. fur φ sind
φ(1) = N(1)1 Φ1 +N
(1)3 Φ3 +N
(1)4 Φ4
φ(2) = N(2)1 Φ1 +N
(2)2 Φ2 +N
(2)2 Φ3 (4.37)
Aus Eigenschaften der Formfunktionen wissen wir, dass entlang der ge-meinsamen Grenze N
(2)2 = N
(1)4 = 0 ist. Unter Berucksichtigung der
Gleichwertigkeit zwischen den Formfunktionen und den Flachenkoordi-naten (4.27) und (4.28) schreiben wir die Gln. (4.37) so um:
φ(1) = L(1)1 Φ1 + L
(1)2 Φ3
φ(2) = L(2)1 Φ1 + L
(2)3 Φ3 . (4.38)
Es sei daraf hingewiesen, dass die Indizien von Flachenkoordinaten mit
den Knotennummern nicht zusammenhangen. Da L(1)3 = L
(2)2 = 0 ist,
ergeben sich die Gln. (4.38) unter Berucksichtigung der Gl. (4.23) in dieser
����������������������������������������������������������������������
����������������������������������������������������������������������
���������������������������������������������������������������������������������
���������������������������������������������������������������������������������
h(2)
h(1)
(1)L 1
(2)L 1
A1
(1)
A1
(2)
y
x
4
1
2
3
c
b
B
(2)
(1)Abb. 4.8: Die Flachenkoordi-
naten L(1)1 und L
(2)1 entlang einer
gemeinsamen Grenze zweier Drei-eckselemente
58 KAPITEL 4. KOORDINATENSYSTEME
Form
φ(1) = L(1)1 Φ1 + (1 − L
(1))1 Φ3
φ(2) = L(2)1 Φ1 + (1 − L
(2))1 Φ3 . (4.39)
Ein Punkt auf der gemeinsamen Grenze der beiden Elemente wurde in Bild
4.8 mit den entsprechenden Teilflachen (schraffiert), die mit L(1)1 und L
(2)1
korrespondieren, dargestellt. Wenn c den Abstand zwischen den Punkt B
und den Knoten 3 ist und b die Lange der Seite 1-3, dann erhlten wir furdie Flachenkoordinaten
L(1)1 =
2A(1)1
2A(1)=
2ch(1)
22bh(1)
2
=c
b
und
L(2)1 =
2A(2)1
2A(2)=
2ch(2)
22bh(2)
2
=c
b= L
(1)1
Damit ist die Kontinuitat von φ auf der gemeinsamen Grenze bewiesen.
Kapitel 5
Zweidimensionale Feldprobleme
In doesem Kapitel befassen wir uns mit der FE-Losung von stationarenFeldproblemen. Spezielle Anwendungsbereiche beinhalten Warmeleitung,
rotationsfreie Stromung und akustistische Vibrationen, die wir nach derEinfuhrung der gemischten Randbedingungen im nachsten Kapitel behan-
deln.
Die Implementierung von FE-Methode kann hier in drei Schritte einge-teilt werden: 1. Festlegung der Eigenschaften von Elementinterpolatio-
nen, 2. Auswertung der Elementmatrizen und 3. Losung eines konkretenProblems. Uber die Eigenschaften von zwei zweidimensionalen Elementenhaben wir bereits im Kapitel 3 gesprochen. Hier konzentrieren wir uns auf
die Elementmatrizen fur die zwei dimensionale Feldgleichung
Dx∂2φ
∂x2+Dy
∂2φ
∂y2−Gφ+Q = 0 . (5.1)
Wir beginnen mit einer kurzen Diskussion uber Dgln, die in (5.1) bein-haltet sind und mehrere physikalische Probleme beschreiben. Im Rest des
Kapitels werden wir den Schwerpunkt auf die Herleitung der Integralglei-chungen fur die Elementmatrizen und ihre Auswertung fur das lineare
Dreieckselement und das bilineare Rechteckselement legen. Anschließendbehandeln wir ein Anwendungsbeispiel mit der einfachsten Randbedin-
gung.
59
60 KAPITEL 5. ZWEIDIMENSIONALE FELDPROBLEME
5.1 Problembeschreibende Differentialgleichungen
Die allgemeine Feldgleichung (5.1) hat viele wichtige Anwendungen in der
physikalischen Wissenschaften, von denen hier einige besprochen werden.Die Randbedingungen werden spater bei den Anwendungen diskutiert.
Der erste Anwendungsbereich ist die Torsion eines nichtkreiformigen Quer-schnittes. In diesem Fall ist die herrschende Dgl.
1
g
∂2φ
∂x2+
1
g
∂2φ
∂y2+ 2θ = 0 , (5.2)
wobei g den Schubmodul des Materials und θ den Verdrehungswinkel dar-stellt. Die Gl. (5.2) erhalt man aus Gl. (5.1) mit Dx = Dy = 1/g,G = 0
und Q = 2θ. Die Variable φ ist eine Spannungsfunktion, und die Schub-spannungen im Schaft hangen mit den Ableitungen von φ nach x und y
zusammen.Verschiedene stromungsmechanische Probleme sind in Dgl. (5.1) einge-
bettet. Fur die ideale wirbelfreie Flussigkeit werden die Potential- und dieStromlinienformulierung durch
∂2φ
∂x2+∂2φ
∂y2= 0 (5.3)
und bzw.∂2ψ
∂x2+∂2ψ
∂y2= 0 (5.4)
beschrieben. Die Stromlinien, ψ, stehen zu den konstanten Potentiallini-
en, φ, senkrecht, und die Komponenten der Geschwindigkeit sind von denAbleitungen von φ oder ψ nach x und y abhangig. In Gln. (5.3) und (5.4)sind Dx = Dy = 1 und G = Q = 0.
Die Stromung von Grundwasser wird ebenso durch Gln. beschrieben, diein Gl. (5.1) eingebettet sind. Das Durchsickern von Wasser unter einem
Damm oder einer Spundwand innerhalb eines geschlossenen Bodenspei-chers wird durch
Dx∂2φ
∂x2+Dy
∂2φ
∂y2= 0 (5.5)
beschrieben, wobei Dx und Dy die Durchlassigkeit des Bodens und φ der
Grundwasserspiegel im Bodenspeicher darstellen. Der Wasserspiegel um
5.1. PROBLEMBESCHREIBENDE DIFFERENTIALGLEICHUNGEN 61
einen Brunnen wahrend des Pumpbetriebes wird durch
Dx∂2φ
∂x2+Dy
∂2φ
∂y2+Q = 0 (5.6)
beschrieben, wobei Q der Term fur eine Punktsenke ist und andere Koef-
fizienten wie in (5.5) sind.
Es gibt zwei Warmeleitungsgleichungen, die in Dgl. (5.1) eingebettet sind:Die Warmeleitung durch Konvektion von einer zweidimensionalen Flosse
zu ihrer umgebenden Flussigkeit wird durch
Dx∂2T
∂x2+Dy
∂2T
∂y2−
2h
tT +
2h
tTf = 0 (5.7)
beschrieben, wobei Dx und Dy die Koeffizienten der Warmeleitfahigkeit
in x- und y-Richtung sind; h ist der Konvektionskoeffizient; t ist die Dickeder Flosse; Tf ist die Temperatur der umgebenden Flussigkeit; und T die
Temperatur der Flosse. Aus dem Vergleich zwischen (5.7) und (5.1) stelltman fest, dass G = 2h/t und Q = 2hTf/t ist.
Die der Gl. (5.7) entsprechende Flosse wurde als dunn angenommen undauf den Warmeverlust uber die Kanten wurde verzichtet. Wenn der Korper
in z-Richtung lang ist, und die Temperatur nur von x- und y-Richtungabhangt, dann wird die Warmeleitung im Allgemeinen von Dgl.
Dx∂2T
∂x2+Dy
∂2T
∂y2= 0 (5.8)
beherrscht, wobei die Koeffizienten wie bei der Gl. (5.7) definiert sind.Konvektive Warmeleitung steht durch Randbedingungen in Verbindungmit der Gl. (5.8).
Die Dgl. (5.1) wird Helmholz’sche Dgl. genannt, wenn G < 0 und Q = 0ist. Ein negatives G fuhrt zur Losung eines Eigenwertproblems. Ein paarphysikalische Probleme, die durch Helmholz’sche Dgl. beschrieben wer-
den, sind die Seiche-Bewegung des Wassers und die akustische Vibratio-nen. Seiche-Bewegung ist eine stehende Wellenbewegung in geschlossenen
flachen Gewassern und wird durch
h∂2w
∂x2+ h
∂2w
∂y2+
4π2
gT 2w = 0 (5.9)
62 KAPITEL 5. ZWEIDIMENSIONALE FELDPROBLEME
beschrieben, wobei h die Wassertiefe im Stillstand, w die Wellenhohe uber
dem Ruheniveau, g die Erdbeschleunigung und T die Schwingungsperiodeist.
Die Vibration einer Flussigkeit innerhalb eines geschlossenen Volumenswird durch
∂2P
∂x2+∂2P
∂y2+w2
c2P = 0 (5.10)
beschrieben, wobei P der zusatzliche Druck uber den Umgebungsdruck,
w die Wellenfrequenz und c die Wellengeschwindigkeit in dem Mediumist. Die in Gl. (5.10) beschriebenen Wellen sind keine Funktion von z-Richtung.
Die Gln. (5.2) bis (5.10) beschreiben neun physikalisch unterschiedlichenProbleme, die in der allgemeinen Dgl. (5.1) enthalten sind. Dies macht
eine Diskussion uber die FE-Losung von Dgl. (5.1) interessant.
5.2 Integralgleichungen fur die Elementmatrizen
Unser unmittelbares Ziel ist die Herleitung von Integralgleichungen, die die
Elementmatrizen fur die Probleme definieren, die in Dgl. (5.1) beinhaltetsind. Der Elementbeitrag zu dem Gleichungssytem ist durch
{
R(e)}
= −
∫
A
[N ]T(
Dx∂2φ
∂x2+Dy
∂2φ
∂y2−Gφ +Q
)
dA (5.11)
gegeben, wobei [N ] der Zeilenvektor von Formfunktionen ist. Da die In-terpolationsfunktion φ(x, y) keine kontinuierliche Ableitung zwischen den
Elementen hat, mussen die Termen mit zweiter Ableitung in (5.11) durchTermen mit erster Ableitung ersetzt werden. Dafur soll man die Produkt-
regel fur Ableitung verwenden. Betrachte die Große
∂
∂x
(
[N ]T∂φ
∂x
)
(5.12)
Ableitung ergibt
∂
∂x
(
[N ]T∂φ
∂x
)
= [N ]T∂2φ
∂x2+∂[N ]T
∂x
∂φ
∂x(5.13)
5.2. INTEGRALGLEICHUNGEN FUR DIE ELEMENTMATRIZEN 63
Umsortieren und Einsetzen fur [N ]T∂2φ/∂x2 in (5.11) ergibt
−
∫
A
[N ]TDx∂2φ
∂x2dA = −
∫
A
Dx∂
∂x
(
[N ]T∂φ
∂x
)
dA+
∫
A
Dx∂[N ]T
∂x
∂φ
∂xdA
(5.14)Nutzung des Green’schen Theorems lasst das erste Integral auf der rechtenSeite von (5.14) durch ein Integral uber den Rand ersetzen:
∫
A
∂
∂x
(
[N ]T∂φ
∂x
)
dA =
∫
Γ
[N ]T∂φ
∂xcos θdΓ (5.15)
Darin ist Γ der Elementrand und θ der Winkel zu der Normale nach außen.
Einsetzen von (5.15) in (5.14) ergibt die endgultige Beziehung fur denTerm mit der zweiten Ableitung als
−
∫
A
Dx[N ]T∂2φ
∂x2dA = −
∫
Γ
Dx[N ]T∂φ
∂xcos θdΓ +
∫
A
Dx∂[N ]T
∂x
∂φ
∂xdA
(5.16)Ahnliche Operationen, angefangen mit
∂
∂y
(
[N ]T∂φ
∂y
)
(5.17)
ergeben
−
∫
A
Dy[N ]T∂2φ
∂y2dA = −
∫
Γ
Dy[N ]T∂φ
∂ysin θdΓ +
∫
A
Dy∂[N ]T
∂y
∂φ
∂ydA
(5.18)Einsetzen von Gln. (5.16) und (5.18) in die Gl. (5.11) ergibt
{
R(e)}
= −
∫
Γ
[N ]T(
Dx∂φ
∂xcos θ +Dy
∂φ
∂ysin θ
)
dΓ
+
∫
A
(
Dx∂[N ]T
∂x
∂φ
∂x+Dy
∂[N ]T
∂y
∂φ
∂y
)
dA
+
∫
A
G[N ]TφdA−
∫
A
Q[N ]TdA (5.19)
Durch Einsetzen von
φ(e) = [N ]{
Φ(e)}
(5.20)
64 KAPITEL 5. ZWEIDIMENSIONALE FELDPROBLEME
fur φ in Gl. (5.19) und Umordnung erhalt man
{
R(e)}
= −
∫
Γ
[N ]T(
Dx∂φ
∂xcos θ +Dy
∂φ
∂ysin θ
)
dΓ
+
(∫
A
(
Dx∂[N ]T
∂x
∂[N ]
∂x+Dy
∂[N ]T
∂y
∂[N ]
∂y
)
dA
)
{
Φ(e)}
+
(∫
A
G[N ]T [N ]dA
)
{
Φ(e)}
−
∫
A
Q[N ]TdA , (5.21)
dass die allgemeine Form{
R(e)}
={
I(e)}
+[
k(e)]{
Φ(e)}
−{
f (e)}
(5.22)
hat, wobei
{
I(e)}
= −
∫
Γ
[N ]T(
Dx∂φ
∂xcos θ +Dy
∂φ
∂ysin θ
)
dΓ , (5.23)
[
k(e)]
=
∫
A
(
Dx∂[N ]T
∂x
∂[N ]
∂x+Dy
∂[N ]T
∂y
∂[N ]
∂y
)
dA+
∫
A
G[N ]T [N ]dA
(5.24)und
{
f (e)}
=
∫
A
Q[N ]TdA (5.25)
ist. Die Variable φ in Gl. (5.23) wurde nicht ersetzt, weil die Große
Dx∂φ
∂xcos θ +Dy
∂φ
∂ysin θ (5.26)
in ableitungsenthaltenen Randbedingungen auftritt und spater ausfuhrlich
behandelt wird. Das erste Integral in Gl. (5.24) kann durch die Definition
[D] =
[
Dx 00 Dy
]
(5.27)
und den Gradientenvektor
{gv} =
{ ∂φ∂x
∂φ∂y
}
=
∂[N ]∂x
∂[N ]∂y
{
Φ(e)}
= [B]{
Φ(e)}
(5.28)
5.3. ELEMENTMATRIZEN: DREIECKSELEMENT 65
in kompackterer Form geschrieben werden. Die erste Reihe von [B] ist die
Ableitung von [N ] nach x und die zweite Reihe die Ableitung von [N ]nach y. Die Matrix [B]T ist gegeben durch
[B]T =[
∂[N ]T
∂x∂[N ]T
∂y
]
(5.29)
Mit den Gln. (5.27), (5.28) und (5.29) erhat man
∫
A
[B]T [D][B]dA =
∫
A
(
Dx∂[N ]T
∂x
∂[N ]
∂x+Dy
∂[N ]T
∂y
∂[N ]
∂y
)
dA . (5.30)
Die Steifigkeitsmatrix fur das Feldproblem mit der Dgl. (5.1) ergibt sich
als[
k(e)]
=
∫
A
[B]T [D][B]dA+
∫
A
G[N ]T [N ]dA (5.31)
Die einzelnen Matrizen werden hier durch [k(e)D ] und [k
(e)G ] bezeichnet. So-
mit erhalten wir[
k(e)]
=[
k(e)D
]
+[
k(e)G
]
. (5.32)
5.3 Elementmatrizen: Dreieckselement
Nach der Erkennung der allgemeinen Form von Steifigkeitsmatrix fur zwei-dimensionale Feldprobleme beschaftigen wir uns mit der Auswertung der
Elementmatrizen fur das Dreieckselement. Die Große φ uber eine Drei-ecksregion ist durch
φ(e) =[
Ni Nj Nk
]
{
Φ(e)}
(5.33)
definiert, wobei die Formfuktionen
Ni =1
2A(ai + bix+ ciy)
Nj =1
2A(aj + bjx+ cjy)
Nk =1
2A(ak + bkx+ cky)
66 KAPITEL 5. ZWEIDIMENSIONALE FELDPROBLEME
und die Koeffizienten a, b und c aus dem Abschnitt 3.1 bekannt sind. Der
Gradinetenvektor fur dieses Element ist
{gv} =
∂Ni
∂x∂Nj
∂x∂Nk
∂x
∂Ni
∂y∂Nj
∂y∂Nk
∂y
{
Φ(e)}
(5.34)
oder
{gv} =1
2A
[
bi bj bkci cj ck
]
{
Φ(e)}
= [B]{
Φ(e)}
(5.35)
Die zwei Matrizen [B] und [D] bestehen nur aus Konstanten, weil bβ, cβ,β = i, j, k Konstanten und Dx und Dy auch Materialkonstanten sind.
Daher ist das erste Integral in Gl. (5.31) leicht auszuwerten. Es ergibt sichzu
[
k(e)D
]
=
∫
A
[B]T [D][B]dA = [B]T [D][B]
∫
A
dA
oder[
k(e)D
]
= [B]T [D][B]A (5.36)
Die Erweiterung des Matrizenproduktes ergibt
[
k(e)D
]
=Dx
4A
b2i bibj bibkbibj b2j bjbkbibk bjbk b2k
+Dy
4A
c2i cicj cickcicj c2j cjckcick cjck c2k
(5.37)
Das zweite Integral in Gl. (5.31) beinhaltet ebenso Formfunktionen. Wennwir annehmen, dass G innerhalb des Elementes konstant ist, dann wird
dieses Integral zu
[
k(e)G
]
=
∫
A
G[N ]T [N ]dA = G
∫
A
Ni
Nj
Nk
[
Ni Nj Nk
]
dA
= G
∫
A
N2i NiNj NiNk
NiNj N2j NjNk
NiNk NjNk N3k
dA
= G
∫
A
L21 L1L2 L1L3
L1L2 L22 L2L3
L1L3 L2L3 L33
dA , (5.38)
5.4. ELEMENTMATRIZEN: RECHTECKSELEMENT 67
weil fur lineares Dreieckelement gilt: Ni = L1, Nj = L2 und Nk = L3. Die
Nutzung der Integrationsformel (4.29) zur Auswertung dieses Integralsliefert
[
k(e)G
]
=GA
12
2 1 11 2 1
1 1 2
(5.39)
Die Elementsteifigkeitsmatrix des Dreieckselements ergibt sich fur diezweidimensionale Feldgleichung (5.1) aus der Summe von (5.37) und
(5.39).Der Kraftevektor des Elementes beinhaltet ebenso Formfunktionen unddie Auswertung von (5.25) ist ganz ahnlich wie [k
(e)G ]. Mit der Annahme,
dass Q innerhalb des Elementes konstant ist, ergibt das Einsetzen
∫
A
Q[N ]TdA = Q
∫
A
Ni
Nj
Nk
dA = Q
∫
A
L1
L2
L3
dA . (5.40)
Nutzung der Integrationsformel (4.29) ergibt
{
f (e)}
=QA
3
1
11
(5.41)
Beispiel 5.1
5.4 Elementmatrizen: Rechteckselement
Die Auswertung der Elementmatrizen fur das Rechteckselement ist nichtso einfach wie fur das Dreieckselement. Jeder Koeffizient der Element-matrizen wird durch Integration eines Polynoms uber ein Gebiet gerech-
net. Die Integrale konnen mittels Formfunktionen in (3.14) oder (3.16)ausgewertet werden. Aus Grund der Ahnlichkeit zwischen den lokalen st-
Koordinaten und den globalen xy-Koordinaten nutzen wir die Formfunk-tionen in (3.14), die im st-Koordinaten entwickelt worden sind.
Dass alle Integrale im xy-Koordinaten definiert wurden, bereitet ein klei-nes Problem. Insbesondere die Gradientenmatrix mit Koeffizienten, die
Ableitungen von Formfunktionen nach x und y beinhalten.
68 KAPITEL 5. ZWEIDIMENSIONALE FELDPROBLEME
Die Gleichung fur Variablenanderung bei Doppelintegralen wird hier nicht
diskutiert, aber ihrer Anwendung fur ein Rechteckselement, das in st-Koordinaten definiert wurde, kann wie folgt zusammengefasst werden. Da
das st-Koordinatensystem parallel zum xy-Koordinatensystem ist und ei-ne Einheitslange in s oder t die gleiche Einheitslange in x und y ist, erhaltman
∫
A
f(x, y)dxdy =
∫
A
f(s, t)dsdt . (5.42)
Genauso wichtig sind die Beziehungen zwischen den Ableitungen. Die Ket-
tenregel ergibt∂Nβ
∂x=∂Nβ
∂sund
∂Nβ
∂y=∂Nβ
∂t(5.43)
Die Formfunktionen in st-Koordinaten haben wir aus (3.14):
Ni = 1 −s
2b−
t
2a+
st
4ab, Nj =
s
2b−
st
4ab
Nk =st
4ab, Nm =
t
2a−
st
4ab. (5.44)
Die Auswertung von [k(e)] und {f (e)} wird durch Betrachtung eines spe-
ziellen Integrals in jedem Fall erlautert. Die einfachsten Integrale gehorendem Elementkraftevektor {f (e)}:
{
f (e)}
=
∫
A
Q[N ]TdA =
∫ 2b
0
∫ 2a
0
Q
Ni
Nj
Nk
Nm
dtds (5.45)
Fur den dritten Koeffizient bekommt man∫ 2b
0
∫ 2a
0
Nkdtds =
∫ 2b
0
∫ 2a
0
st
4abdtds
=
∫ 2b
0
st2
8ab
∣
∣
∣
∣
2a
0
ds =
∫ 2b
0
as
2bds =
A
4(5.46)
Die drei weiteren Koeffizienten liefern das gleiche Ergebnis und daher ist
{
f (e)}
=QA
4
1
11
1
. (5.47)
5.4. ELEMENTMATRIZEN: RECHTECKSELEMENT 69
Das Integral im Zusammenhang mit [k(e)G ] ist
[
k(e)G
]
=
∫
A
G[N ]T [N ]dA
=
∫
A
G
N2i NiNj NiNk NiNm
NiNj N2j NjNk NjNm
NiNk NjNk N2k NkNm
NiNm NjNm NkNm N2m
dA (5.48)
Fur den ausgewahlten Term N2k erhalten wir
∫ 2b
0
∫ 2a
0
(
st
4ab
)2
dtds =
∫ 2b
0
∫ 2a
0
s2t2
16a2b2dtds =
4ab
9=A
9(5.49)
Den kompletten Satz von Koeffitzienten erhalten wir als
[
k(e)G
]
=GA
36
4 2 1 22 4 2 1
1 2 4 22 1 2 4
(5.50)
Die Auswertung von [k(e)D ] erfordert Ableitungen von Formfunktionen. Die
Gradientenmatrix lautet
[B] =
∂Ni
∂x∂Nj
∂x∂Nk
∂x∂Nm
∂x
∂Ni
∂y∂Nj
∂y∂Nk
∂y∂Nm
∂y
(5.51)
Setzen wir die Beziehungen in (5.43) hier ein, lasst sich [B] in Termen vons und t ausdrucken. Die Ableitung von Formfunktionen fuhrt dann zu
[B] =1
4ab
[
−(2a− t) (2a− t) t −t−(2b− s) −s s (2b− s)
]
(5.52)
Den Koeffizient in der ersten Reihe und der ersten Spalte von [k(e)D ] erhalt
man nach Durchfuhrung des Produktes [B]T [D][B] in folgender Form
Dx
16a2b2(2a− t)2 +
Dy
16a2b2(2b− s)2 , (5.53)
und das dazu assoziierte Integral lautet∫ 2b
0
∫ 2a
0
Dx
16a2b2(2a− t)2dtds+
∫ 2b
0
∫ 2a
0
Dy
16a2b2(2b− s)2dtds . (5.54)
70 KAPITEL 5. ZWEIDIMENSIONALE FELDPROBLEME
Durchfuhrung der Integration resultiert in
Dxa
3b+Dyb
3a. (5.55)
Das vollstandige Ergebnis fur [k(e)D ] ist dann
[
k(e)D
]
=Dxa
6b
2 −2 −1 1
−2 2 1 −1−1 1 2 −2
1 −1 −2 2
+Dyb
6a
2 1 −1 −2
1 2 −2 −1−1 −2 2 1−2 −1 1 2
(5.56)
Die Elemntsteifigkeitsmatrix [k(e)] fur das Rechteckselement ergibt sich
aus der Summe von (5.50) und (5.56).
5.5 Torsion eines nichtkreiformigen Querschnitts
Nach der Auswertung von Elementmatrizen fur das lineare Dreiecksele-
ment und das bilineare Rechteckselement befassen wir uns nun mit ihrerAnwendung zur numerischen Losung eines praktischen Problems. Wir wol-
len die Schubspannungen in einem Stab mit quadratischem Querschnittermitteln, der durch Torsion belastet wird. Diese Aufgabe wurde gewahlt,
weil sie die einfachste Randbedingung hat: φ ist am gesamten Rand gleichNull.
5.5.1 Allgemeine Theorie
Zur Berechnung von Schubspannungen im Torsionsstab mit nicht-
kreisformigem Querschnitt gibt es zwei Theorien: Die erste ist von St.
Venant und die zweite von Prandtl. Hier wird die Theorie von Prandtl
zu Grunde gelegt. Nach dieser Theorie konnen die Schubspannungen ineinem nichtkreisformigen Stab, der mit einem Torsionsmoment T um die
z-Achse belastet wird, mit
τzx =∂φ
∂yund τzy = −
∂φ
∂x(5.57)
5.5. TORSION EINES NICHTKREIFORMIGEN QUERSCHNITTS 71
τ zxy
x
z
ϕ
τ zxτ zy
(a)
(b)
y
z
x
T
Abb. 5.1: (a) Komponenten derSchubspannung in einem nicht-kreisformigen Querschnitt unterTorsionslast (b) Die Flache desSpannungsfeldes φ(x, y) und dieSchubkomponente
berechnet werden, wobei φ(x, y) eine Spannungsfunktion ist. Die Problem-
beschreibende Dgl. ist
1
g
∂2φ
∂x2+
1
g
∂2φ
∂y2+ 2θ = 0 (5.58)
und uberall am Rand (Querschnittsumfang) gilt
φ = 0 . (5.59)
Die physikalischen Parameter in Dgl. (5.58) sind der Schubmodul,
g [N/cm2], und der Verdrehwinkel, θ [rad/cm]. Die Prandtl-Formulierunghat kein Torsionsmoment, T [N.cm], in der herrschenden Dgl., stattdessen
wird das Torsionsmoment durch
T = 2
∫
A
φdA (5.60)
berechnet, wenn die Spannungsfunktion φ(x, y) bekannt ist.Die Spannungsfunktion stellt eine Flache dar, die uber den Stabquer-
schnitt verlauft und ihn abdeckt (s. Bild 5.1). Das Torsionsmoment istzu dem Volumen unter dieser Flache proportional, wahrend die Schub-
spannungen mit den Flachengradienten in x- und y-Richtung zusam-menhangen.
Wenn der Stab aus einem Material besteht, schreibt man die Gl. (5.58)
72 KAPITEL 5. ZWEIDIMENSIONALE FELDPROBLEME
normalerweise in der Form
∂2φ
∂x2+∂2φ
∂y2+ 2gθ = 0 , (5.61)
die im Vergleich mit der Gl. (5.1) folgende Konstanten besitzt:Dx = Dy =
1, G = 0 und Q = 2gθ.
5.5.2 Verdrehung des quadratischen Torsionsstabs
Anhand dieses Beispiels wollen wir die Auswertung von Elementmatri-
zen und ihrer Assembelierung zu einem System von linearen Gleichungenerlautern. Der Torsionsstab (s. Bild 5.2) hat vier Symmetrieachsen und
deshalb reicht es, ein achtel von seinem Querschnitt zu analysieren. Die-ser Teilquerschnitt wird durch drei Elemente diskretisiert (Bild 5.2b). DreiElemente reichen fur eine genaue Losung nicht aus, fur die Darstellung der
Berechnungsweise schon.Die Element-Knoten-Numerierung ist wie folgt
e i j k m
1 1 2 4
2 2 3 5 43 4 5 6
Dreieckselemente 1 und 3 haben dieselbe Orientierung und dieselben Di-
mensionen, daher sind ihre Elementmatrizen identisch. Die Elementma-trizen fur ein Dreieckselement sind laut Gln. (5.37) und (5.41)
[
k(e)D
]
=Dx
4A
b2i bibj bibkbibj b2j bjbkbibk bjbk b2k
+Dy
4A
c2i cicj cickcicj c2j cjckcick cjck c2k
(5.62)
und
{
f (e)}
=2gθA
3
11
1
. (5.63)
5.5. TORSION EINES NICHTKREIFORMIGEN QUERSCHNITTS 73
θ = 0,01 Grad/cmg = 8E+6 N/cm2
������������������������������
������������������������������
0,25
cm
0,25
cm
(1)(2)
(3)
1 2 3
4 5
6
1 cm
(a) (b)
1 cm 0,25 cm 0,25 cm
Abb. 5.2: (a) Der Teilquerschnitt aus Symmetriegrunden (b) Netz der finiten Elemente
Fur das Rechteckselement mit 2a = 2b bekommt man aus Gln. (5.56) und(5.47)
[
k(e)]
=1
6
4 −1 −2 −1−1 4 −1 −2
−2 −1 4 −1−1 −2 −1 4
(5.64)
und
{
f (e)}
=2gθA
3
1
11
1
. (5.65)
Nun werten wir die Gl. (5.62) fur die Elemente 1 und 3 aus. Die Element-flache betragt 1
32und daher 4A(1) = 1
8. Die Koeffizienten b und c sind
b(1)1 = Y2 − Y4 = −0, 25, c
(1)1 = X4 −X2 = 0
b(1)2 = Y4 − Y1 = 0, 25 c
(1)2 = X1 −X4 = −0, 25
b(1)4 = Y1 − Y2 = 0 c
(1)4 = X2 −X1 = 0, 25
Durch Einsetzen von Koeffizienten in (5.62) unter Berucksichtigung von
1/4A(1) = 8 und der Tatsache, dass die Produktformen bibj und cicj ent-
74 KAPITEL 5. ZWEIDIMENSIONALE FELDPROBLEME
weder Null oder ± 116 ergeben, bekommen wir
[
k(1)]
=8
16
1 −1 0−1 1 0
0 0 0
+8
16
0 0 00 1 −1
0 −1 1
. (5.66)
Addition der beiden Matrizen ergibt
[
k(1)]
=1
2
1 −1 0−1 2 −1
0 −1 1
=[
k(3)]
. (5.67)
Fur den Elementkraftevektor {f (e)} mussen wir zuerst den Term 2gθ be-
rechnen. Einsetzen von Werten aus dem Bild (5.2) ergibt
2gθ = 2(8.106)(0, 01)( π
180
)
= 2790.
Substitution von 2790 und A(1) in (5.63) fuhrt zu
{
f (1)}T
=[
29, 1 29, 1 29, 1]
={
f (3)}T
(5.68)
Fur das Element 2 ist der Elementkraftevektor mit A(2) = 116 und 2gθ =
2790 gleich
{
f (2)}T
=[
43, 6 43, 6 43, 6 43, 6]
. (5.69)
Die Elementmatrizen sind hier zusammengefasst. Die Knotennummer be-
zeichnen die Zeilen und Spalten von [K] und {F} auf die der jeweilige
5.5. TORSION EINES NICHTKREIFORMIGEN QUERSCHNITTS 75
Koeffizient zu addieren ist.
[
k(1)]
=1
2
1 2 4
1−1
0
−12
−1
0−1
1
,{
f (1)}
=
29, 129, 1
29, 1
12
4
(5.70)
[
k(2)]
=1
6
2 3 5 4
4−1
−2−1
−14
−1−2
−2−1
4−1
−1−2
−14
,{
f (1)}
=
43, 643, 6
43, 643, 6
23
54
(5.71)
[
k(1)]
=1
2
4 5 6
1−1
0
−12
−1
0−1
1
,{
f (1)}
=
29, 129, 1
29, 1
45
6
(5.72)
Das Aufaddieren der Elementbeitrage nach direkter Steifigkeitsmethode
und Multiplizierung mit 6 fuhrt zum Gleichungssystem des Problems
3 −3 0 0 0 0−3 10 −1 −4 −2 0
0 −1 4 −2 −1 00 −4 −2 10 −4 0
0 −2 −1 −4 10 −30 0 0 0 −3 3
Φ1
Φ2
Φ3
Φ4
Φ5
Φ6
−
175436
262611
436175
=
00
00
00
(5.73)
Die Knotengroßen Φ3, Φ5 und Φ6 sind gleich Null, weil sie auf dem außeren
Rand liegen. Daher konnen die Gleichungen drei, funf und sechs gestri-chen werden. Der Einfluss von Spalten drei, funf und sechs muss auf den
Kraftvektor [F ] berucksichtigt werden. Da aber Φ3 = Φ5 = Φ6 = 0 ist,liefern diese Spalten keinen Beitrag und das modifizierte Gleichungssytem
ist
3 −3 0
−3 10 −40 −4 10
Φ1
Φ2
Φ4
−
175
436611
=
0
00
.
Losung des Gleichungssystems liefert
Φ1 = 217, Φ2 = 159, und Φ4 = 125
76 KAPITEL 5. ZWEIDIMENSIONALE FELDPROBLEME
zyτ
zxτ
(b)
1 2 3
4 5
6
y
232
136
568
68
500
0
636
x
Φ = 1592
Φ = 05
Φ = 06
Φ = 03
Φ = 1254
Φ = 2171
xy
φ
(a)
Abb. 5.3: (a) Die φ(x, y)-Flache fur ein achtel des Querschnittes vom Torsionsstab mitquadratischem Querschnitt (b) Die Scherspannungswerte in den Elementmittelpunkten
Die Flache von φ fur diesen Satz von Knoten ist in Bild 5.3 dargestellt.
Bestimmung der Knotengroßen ist der Hauptteil der Losung, aber es gibtauch eine Reihe von Elementgroßen, die anhand Knotengroßen zu berech-
nen sind. In diesem Beispiel sind die Schubspannungen und das Torsions-moment von Interesse.
5.5.3 Komponenten der Schubspannung
Die Komponenten der Schubspannung hangen mit den Gradienten derSpannungsfunktion wie folgt zusammen:
τzx =∂φ
∂yund τzy = −
∂φ
∂x(5.74)
Der Gradientenvektor ist fur das lineare Dreieckselement durch Gl. (5.35)
als
{gv} =1
2A
[
bi bj bkci cj ck
]
Φi
Φj
Φk
(5.75)
5.5. TORSION EINES NICHTKREIFORMIGEN QUERSCHNITTS 77
gegeben und ist uber das Element konstant. Der Gradientenvektor ist fur
das Rechteckselement durch Gl. (5.52) als
{gv} =1
4ab
[
−(2a− t) (2a− t) t −t−(2b− s) −s s (2b− s)
]
Φi
Φj
Φk
Φm
(5.76)
gegeben und ist uber das Element nicht konstant. Die Elementflache und
die Koeffizienten b und c sind fur beide Dreieckselemente gleich und bereitsvor der Gl. (5.66) berechnet worden. Mit Berucksichtigung der Knoten-
großen Φ1, Φ2 und Φ4 erhalten wir
{
gv(1)}
=16
4
[
−1 1 0
0 −1 1
]
217159125
=
{
−232
−136
}
und
τ (1)zx =
∂φ
∂y= −136 N/cm2
τ (1)zy = −
∂φ
∂x= 232 N/cm2
fur das Element 3 bekommt man auf gleicher Art und Weise
τ (3)zx = 0 und τ (3)
zy = 500 N/cm2
Die Gradiantenwerte sind innerhalb des Rechteckselements nicht konstant.Wir berechnen die Komponeneten der Schubspannung am Knoten 3, wo
τxy am großten ist. Die lokalen Koordinaten vom Knoten 3 sind: s = 2bund t = 0, dabei ist 2a = 2b = 0, 25. Der Gradientenvektor ist damit
{
gv(2)}
=16
4
[
−1 1 0 00 −1 1 0
]
Φ2
Φ3
Φ5
Φ4
=
[
−4 4 0 00 −4 4 0
]
159
00
125
=
{
−6360
}
78 KAPITEL 5. ZWEIDIMENSIONALE FELDPROBLEME
und
τ (2)zx = 0 und τ (2)
zy = 636 N/cm2
Die berechneten Schubspannungen sind anhand der Werte von τzx und τzy
am Mittelpunkt der Elemente in Bild ?? dargestellt.Es gibt mindestens zwei Wege fur die Verbesserung der Spannungswer-
te bei diesem Beispiel. Der erste Weg geht uber eine hohere Anzahl vonElementen (feinere Diskretisierung). Wenn die Große des Elementes ab-
nimmt, wird die Annahme uber die konstante Spannung uber das Elementrealistischer. Der zweite Weg geht uber den Einsatz von Elementen mit
Ansatzfunktionen von hoherer Ordnung (quadratisch oder qubisch). Imletzten Fall erhalt man durch Ableitung Gradienten, die Funktionen vom
Ort innerhalb des Elementes sind.
5.5.4 Berechnung des Torsionsmoments
Eine weitere interessante Große ist das Torsionsmoment T , das in Gl.
(5.60) definiert wurde und fur unser Beispiel mit
T =n
∑
e=1
2
∫
A(e)
φ(e)dA . (5.77)
aquivalent ist. Das Integral∫
A(e)
φ(e)dA
ist fur das lineare Dreieckselement gleich∫
A(e)
φ(e)dA =A
3(Φi + Φj + Φk) (5.78)
und fur das Rechteckselement gleich∫
A(e)
φ(e)dA =A
4(Φi + Φj + Φk + Φm). (5.79)
Das Torsionsmoment T ist die Summe aller Elementbeitrage:
T = T (1) + T (2) + T (3)
5.5. TORSION EINES NICHTKREIFORMIGEN QUERSCHNITTS 79
Die Elementbeitrage sind hier
T (1) = 2(1
32)(
1
3)(Φ1 + Φ2 + Φ4)
=1
48(217 + 159 + 125) = 10, 4
T (2) = 2(1
16)(
1
4)(Φ2 + Φ3 + Φ5 + Φ4) = 8, 88
T (3) = (1
48)(Φ4 + Φ5 + Φ6) = 2, 6
und das Torsionsmoment fur den Querschnittsachtel ist
T = 10, 4 + 8, 88 + 2, 6 = 21, 9 N.cm (5.80)
Auf den Gesamtquerschnitt wirkendes Torsionsmoment ist damit gleich8(21, 9) ≃ 175 N.cm.
Ein Torsionsmoment von 175 N.cm ist notwendig um einen Stab mit derLange 100 cm und dem quadratischen Querschnitt der Kantenlange 1 cm
um 1◦ zu verdrehen. Die Genauigkeit dieses Ergebnisses ist wegen derGrobheit des FE-Netzes fraglich. In der Tat liegt unsere Losung um 11%unter dem theoretischen Wert1 von 196 N.cm.
Solange wir den Verdrehungswinkel, θ, nur nach Schatzung anneh-men, kann der rechnerische Wert vom Torsionsmoment nur zufallig dem
tatsachlich aufgebrachten Moment entsprechen. Die korrekten Werte vonτzx, τzy und θ sind durch Skalierung mit dem Faktor Tact/Tcal zu berech-
nen. Als Beispiel nehmen wir an, dass das belastende Torsionsmoment inunserer Aufgabe 250 N.cm ist, so wird der wahre Verdrehungswinkel
θtrue =
(
Tact
Tcal
)
θassumed =
(
250
175
)
(0, 01) = 0, 0143 rad/cm (5.81)
sein. Die Schubspannungen werden auf gleicher Art skaliert. Der großte
Wert von τzy fur das tatsachliche Torsionsmoment ergibt sich dann als
τzy =
(
Tact
Tcal
)
τzy =
(
250
175
)
(636) = 909 N/cm2. (5.82)
1Die Beziehung zwischen dem Torsionsmoment und der Verdrehung ist fur einen quadratischen Quer-schnitt der Kantenlange 2a durch T = 0, 1406gθ(2a)4 gegeben (Timoshenko and Goodier, 1970).
Kapitel 6
2D-Feldprobleme mit gemischten
Randbedingungen
Im Gegensatz zur Torsion eines nichtkreisformigen Querschnittes, wobei
der Spannungsverlauf, φ(x, y), uber den Querschnitt am gesamten Randdefiniert ist, haben die meisten physikalischen Feldprobleme eine Mischung
von Randbedingungen. An einem Teil des Randes ist die Feldvariableφ(x, y) definiert und an einem anderen Teil sind Randbedingungen vor-
gegeben, die von den Ableitungen der Feldvariable, ∂φ/∂x und ∂φ/∂y,abhangen. In diesem Abschnitt wollen wir uns mit den ableitungsenthal-tenen Randbedingungen befassen. Hier beschranken wir uns weiterhin auf
zweidimensionale Feldprobleme, fur die wir die allgemeine Form der her-schenden Dgl., namlich die Gl. (5.1), kennen gelernt haben. Fur diese
Problemklasse haben wir bereits in den Abschnitten 5.3 und 5.4 die Ele-mentmatrizen fur das lineare Dreieckselement und das bilineare Recht-
eckselement hergeleitet.
Unter den Elementmatrizen (Steifigkeitsmatrix, Kraftevektor und Inter-elementvektor) haben wir bis jetzt auf den Einfluss von Interelementvek-tor,
{
I(e)}
, in der Restegleichung (5.22) verzichtet, weil die Unabhangig-
keit der Randbedingungen von den Ableitungen von φ vorausgesetzt war.Nun ist diese Voraussetzung nicht mehr erfullt. Da der Vektor
{
I(e)}
die
Moglichkeit zur Implementierung der gemischten Randbedingungen bie-tet, beschaftigen wir uns zuerst mit diesem Vektor und daraus resultieren-
den Einflussen auf die Steifigkeitsmatrix und den Kraftevektor im Glei-chungssystem. Anschließend behandeln wir verschiedene Anwendungsbe-
reiche.
81
82KAPITEL 6. 2D-FELDPROBLEME MIT GEMISCHTEN RANDBEDINGUNGEN
6.1 Ableitungsenthaltene Randbedingungen
Die zwei Randbedingungsarten fur zweidimensionale Feldprobleme sind
im Bild 6.1 schematisch dargestellt. Uber einen Teil des Randes, Γ1, ist φdefiniert, und uber den Restteil des Randes, Γ2, ist eine Randbedingung
der Art
Dx∂φ
∂xcos θ +Dy
∂φ
∂ysin θ = −Mφb + S (6.1)
vorgegeben. Fur Dx = Dy reduziert sich diese Gleichung auf
Dx∂φ
∂n= −Mφb + S, (6.2)
wobei ∂φ/∂n die Ableitung von φ bezogen auf die Normale zum Randdarstellt. In beiden Fallen steht φb fur den Wert von φ auf dem Rand Γ2
und ist vorerst unbekannt. Die Gln. (6.1) und (6.2) vereinfachen sich mitM = S = 0 zu
Dx∂φ
∂xcos θ +Dy
∂φ
∂ysin θ = 0 (6.3)
und
Dx∂φ
∂n= 0. (6.4)
Die Bedingung M = S = 0 gilt fur isolierte und undurchlassige Rander
und Symmetrieachsen. Situationen sind auch moglich, in denen M oderS, aber nicht beide, Null sind.
Die ableitungsenthaltenen Randbedingungen werden mit Hilfe des Inter-elementvektors
{
I(e)}
in FE-Analyse der Feldprobleme eingearbeitet. Die-
x
y
Γ1
Γ2
φ(x, y)φdie Ableitungen von enthalten.
Γ2Auf gelten Randbedingungen,
Γ1Auf ist bekannt.φ
θ
n
Abb. 6.1: Feldprobleme mit gemischten Randbedingungen.
6.1. ABLEITUNGSENTHALTENE RANDBEDINGUNGEN 83
ser Vektor ist nach Gl. (5.23) wie folgt definiert:
{
I(e)}
= −
∫
Γ
[N ]T(
Dx∂φ
∂xcos θ +Dy
∂φ
∂ysin θ
)
dΓ, (6.5)
wobei das Integral uber den Elementrand und gegen den Uhrzeigersinnauszuwerten ist. Fur das Dreieckselement ist das Integral in (6.5) die Sum-
me von drei Integralen, ein Integral fur jede Elementseite.Wir konnen den Vektor
{
I(e)}
in zwei Komponenten zerlegen
{
I(e)}
={
I(e)bc
}
+{
I(e)i
}
. (6.6)
Die erste Komponente steht fur die Elementseite mit vorgegebener Rand-
bedingung, Γbc, und daher
{
I(e)bc
}
= −
∫
Γbc
[N ]T(
Dx∂φ
∂xcos θ +Dy
∂φ
∂ysin θ
)
dΓ. (6.7)
Die zweite Komponente, der Vektor{
I(e)i
}
, beinhaltet Integrale von (6.5),
die uber die Elementseiten ohne vorgegebene Randbedingung auszuwerten
sind. Diese Integrale fuhren zu den interelementaren Bedingungen, dieerfullt werden mussen, bevor der Galerkin’sche Rest verschwindet.
Unter Berucksichtigung der Gl. (6.1) erhalten wir aus (6.7)
{
I(e)bc
}
=
∫
Γbc
[N ]T (Mφb − S)dΓ. (6.8)
φb ist durch folgende Elementgleichung gegeben
φ(e) = [N ]{
Φ(e)}
. (6.9)
Das Einsetzen von (6.9) in (6.8) ergibt
{
I(e)bc
}
=
∫
Γbc
[N ]T(
M [N ]{
Φ(e)}
− S)
dΓ, (6.10)
das man in zwei Integrale schreiben kann:
{
I(e)bc
}
=
(∫
Γbc
M [N ]T [N ]dΓ
)
{
Φ(e)}
−
∫
Γbc
S[N ]TdΓ. (6.11)
84KAPITEL 6. 2D-FELDPROBLEME MIT GEMISCHTEN RANDBEDINGUNGEN
������������������������������������������������������������������������������������������
������������������������������������������������������������������������������������������
Dx∂φ
∂x= Mφb − S
−Dx∂φ
∂x= Mφb − S
Dy∂φ
∂y= Mφb − S
−Dy∂φ
∂y= Mφb − S
R :1
R :4
R4
R4
R3
R3
R1 R2 R1 R2
R :3
R :2
Abb. 6.2: Randbedingungen fur eine rechteckige Region mit Aussparung.
Durch Randbedingung am Randelement entstehen also zwei Komponen-ten mit unterschiedlichen Auswirkungen. Die eine Komponente addiert
sich auf Elementsteifigkeitsmatrix[
k(e)]
, weil sie mit{
Φ(e)}
multipliziertwird. Die andere komponente addiert sich auf den Elementkraftevektor{
f (e)}
. Die zwei Komponenten von{
I(e)bc
}
in (6.11) werden wie folgt dar-
gestellt:{
I(e)bc
}
=[
k(e)M
]{
Φ(e)}
−{
f(e)S
}
, (6.12)
wobei die erste Komponente
[
k(e)M
]
=
∫
Γbc
M [N ]T [N ]dΓ (6.13)
und die zweite{
f(e)S
}
=
∫
Γbc
S[N ]TdΓ (6.14)
ist. Bevor wir mit der Auswertung dieser Elementintegrale fortfahren, ist
darauf hinzuweisen, dass die Randbedingung (6.1) verschiedene Formenannimmt, wenn wir die Außen- und Innenseite eines rechteckigen Gebietes
mit einem Loch betrachten. Die unterschiedlichen Formen fur
−
(
Dx∂φ
∂xcos θ +Dy
∂φ
∂ysin θ
)
(6.15)
sind im Bild 6.2 gegeben. Die mitM und S assoziierten Vorzeichen mussenanhand der Informationen aus dem Bild 6.2 bestimmt werden. Die art die-
ser Bestimmung wird in einem Anwendungsbeispiel gezeigt. Es sei darauf
6.1. ABLEITUNGSENTHALTENE RANDBEDINGUNGEN 85
hingewiesen, dass θ der Winkel zwieschen der x-Achse und der Außennor-
male (n im Bild 6.1) ist.
6.1.1 Auswertung von Elementintegralen
Die Integrale in (6.13) und (6.14) sind fur jedes zweidimensionales Ele-
ment gultig und konnen ausgewertet werden, sobald die Formfunktionenbekannt sind. Diese Integrale wollen wir fur das lineare Dreieckselement
und das bilineare Rechteckselement auswerten und beginnen wir mit derGl. (6.14) und dem Rechtseckselement, weil es die einfachste Kombination
ist.Nehmen wir an, dass S uber die ij-Seite eines Rechteckelementes vorgege-ben ist und die Dicke des Elementes einer Langeneinheit entspricht, somit
bekommen wir
∫
Γbc
S[N ]TdΓ =
∫ b
−b
S
Ni
Nj
Nk
Nm
dq, (6.16)
wobei die Formfunktionen im qr-Koordinatensystem durch die Gln. (3.16)
gegeben sind. Unter Berucksichtigung, dass fur die ij-Seite Nk = Nm = 0ist, erhalten wir durch Einsetzen von Formfunktionen und r = −a:
{
f(e)S
}
=
∫ b
−b
S
2b
(b− q)(b+ q)
00
dq =SLij
2
11
00
. (6.17)
Die Große S ist mit der Lange der ij-Seite, Lij = 2b, multipliziert und
dann zwischen den beiden Knoten auf der ij-Seite geteilt worden. Es blei-ben noch drei weitere Auswertungen ubrig, eine fur jede weitere Seite.Ergebnisse dieser Auswertungen sind fur die jk-, km- und im-Seiten:
{
f(e)S
}
=SLjk
2
0
11
0
,SLkm
2
0
01
1
undSLim
2
1
00
1
(6.18)
86KAPITEL 6. 2D-FELDPROBLEME MIT GEMISCHTEN RANDBEDINGUNGEN
Wenn S fur mehr als eine Seite eines Elementes vorgegeben ist, sind die
Werte von{
f(e)S
}
fur entsprechende Seiten aufeinander zu addieren. Die
Auswertung von (6.14) fur das Dreieckselement ergibt:
{
f(e)S
}
=SLij
2
11
0
SLjk
2
01
1
undSLik
2
10
1
(6.19)
wobei Lij, Ljk und Lik die Langen von ij-, jk- bzw. ik-Seiten sind und
nicht die Flachenkoordinaten. Das erste Ergebnis in (6.19) wurde wie folgtberechnet: Fur gegebene Seite ij haben wir
{
f(e)S
}
=
∫
Γbc
S[N ]TdΓ = Lij
∫ 1
0
S
Ni
Nj
Nk
dℓ2 (6.20)
Da entlang der ij-Seite Nk = 0 ist, erhalt man
{
f(e)S
}
= Lij
∫ 1
0
S
Ni
Nj
0
dℓ2 = Lij
∫ 1
0
S
ℓ1ℓ20
dℓ2, (6.21)
weil die Formfunktionen Ni und Nj sich auf der ij-Seite auf
Ni = L1 = ℓ1 und Nj = L2 = ℓ2 (6.22)
reduzieren. Die Integration findet entlang einer Linie statt, daher konnen
wir die Fakultatsformel (4.17) anwenden und direkt zu den Ergebnissenin (6.19) gelangen.
Die mit [k(e)M ] assoziierten Integralen werden auf gleicher Art und Weise
ausgewertet. Der Unterschied liegt in der Anzahl der Terme, die zu beruck-sichtigen sind. Fur das Rechteckselement fuhrt das Integral in (6.13) zu
[
k(e)M
]
=
∫
Γbc
M
N2i NiNj NiNk NiNm
NiNj N2j NjNk NjNm
NiNk NjNk N2k NkNm
NiNm NjNm NkNm N2m
dΓ (6.23)
6.1. ABLEITUNGSENTHALTENE RANDBEDINGUNGEN 87
Mit der Annahme, dass M auf der ij-Seite definiert ist, haben wir Nk =
Nm = 0 und (6.23) wird zu
[
k(e)M
]
=
∫ b
−b
M
N2i NiNj 0 0
NiNj N2j 0 0
0 0 0 00 0 0 0
dq. (6.24)
Die Auswertung der einzelnen Komponente unter Berucksichtigung vonr = −a ergibt
∫ b
−b
N2i dq =
∫ b
−b
(b− q)2
4b2dq =
2b
3=Lij
3(6.25)
∫ b
−b
NiNjdq =
∫ b
−b
(b− q)(b+ q)
4b2dq =
2b
6=Lij
6(6.26)
und∫ b
−b
N2j dq =
∫ b
−b
(b+ q)2
4b2dq =
2b
3=Lij
3. (6.27)
Durch Einsetzen der Integralergebnisse erhaltren wir
[
k(e)M
]
=MLij
6
2 1 0 0
1 2 0 00 0 0 0
0 0 0 0
(6.28)
Es gibt fur jede weitere Seite ein Ergebnis fur[
k(e)M
]
. Diese drei Ergebnisse
sind
[
k(e)M
]
=MLjk
6
0 0 0 0
0 2 1 00 1 2 0
0 0 0 0
(6.29)
[
k(e)M
]
=MLkm
6
0 0 0 0
0 0 0 00 0 2 1
0 0 1 2
(6.30)
88KAPITEL 6. 2D-FELDPROBLEME MIT GEMISCHTEN RANDBEDINGUNGEN
und
[
k(e)M
]
=MLim
6
2 0 0 10 0 0 00 0 0 0
1 0 0 2
. (6.31)
wobei Ljk, Lkm und Lim die Lange der entsprechenden Seiten darstellen.
Die Auswertung von (6.13) fur das Dreieckselement ergibt
[
k(e)M
]
=MLij
6
2 1 01 2 0
0 0 0
(6.32)
[
k(e)M
]
=MLjk
6
0 0 0
0 2 10 1 2
(6.33)
und[
k(e)M
]
=MLik
6
2 0 10 0 0
1 0 2
. (6.34)
6.2 Punktquellen und Punktsenken
Mit dem Konzept von Punktquelle und Punktsenke wird eine physikalischwichtige Situation beschrieben. Eine Punktquelle oder -senke existiert,
wenn Q in Dgl. (5.1) in einem sehr kleinen Bereich konzentriert ist. Bei-spiele aus dem Bereich der Linienquellen sind unterirdische Dampf- und
Heißwasser-Rohranlagen, sowie konduktive drahtformige Heizkorper in-nerhalb eines massiven Teils. In solchen Fallen ist die Schnittflache des Lei-tungsrohres oder des Heizelementes im Vergleich mit ihrem Umgebungs-
medium sehr klein und wird als einen Punkt auf einer Ebene betrachtet.Punktsenke treten bei Grundwasserproblemen auf, z.B. bei Auspumpen
von Wasser aus einem Brunnen.Unsere Diskussion hier ist auf zweidimensionalen Elementen basiert. Die
Prozedur lasst sich jedoch zur Behandlung von ein- oder dreidimensiona-lem Element einfach modifizieren. Betrachten wir das Dreieckselement im
Bild 6.3 mit einer Quelle Q∗ am Ort (X0, Y0). Da die Quelle an einem
6.2. PUNKTQUELLEN UND PUNKTSENKEN 89
X0 , Y0( )
Q*
i
y
x
j
k
Abb. 6.3: Ein Element mit einerPunktquelle oder Punktsenke, Q∗
am (X0, Y0).
Punkt fest liegt, ist Q uber das gesamte Element nicht konstant, son-
dern eine Funktion von x und y. Mit Hilfe von Impuls-Einheitsfunktionenδ(x−X0) und δ(y − Y0) kann man fur die Feldvariable Q(x, y) innerhalbdes Elementes schreiben
Q = Q∗δ(x−X0)δ(y − Y0) (6.35)
Das Integral{
f(e)Q
}
=
∫
A
Q[N ]TdA (6.36)
fuhrt zu
{
f(e)Q
}
= Q∗
∫
A
Ni
Nj
Nk
δ(x−X0)δ(y − Y0)dxdy. (6.37)
Das Integral einer Große, die mit einer Impulsfunktion multipliziert wurde,
ist dem Wert der selben Große am X0 und Y0 gleich. Daher gilt
{
f(e)Q
}
= Q∗
Ni(X0, Y0)
Nj(X0, Y0)Nk(X0, Y0)
. (6.38)
Das Verhaltnis, zu dem Q∗ jedem Knoten zugeteilt wird, beruht auf denrelativen Werten von Ni, Nj und Nk, die mit den Koordinaten der Quel-
le ermittelt worden sind. Da die Summe der Formfunktionen fur jedenPunkt innerhalb des Elementes den Wert 1 ergibt, haben wir nur die Q∗
zuzuweisen.
90KAPITEL 6. 2D-FELDPROBLEME MIT GEMISCHTEN RANDBEDINGUNGEN
Q*
j
α
k
i
Abb. 6.4: Eine Punktquelle aneinem Knoten.
Die beste Lage fur eine Quelle oder Senke ist an einem Knoten. DieseLage andert die Ergebnisse in (6.38). Nehmen wir an, dass die Quelle am
Knoten j (s. Bild 6.4) liegt, dann gilt Ni = Nk = 0 und
{
f(e)Q
}
= Q∗
0
10
. (6.39)
Der Betrag von Q∗ muss korrigiert werden, wenn die Quelle (oder Senke)mehr als einem Element angehort. Der Wert der Quelle wird zwischen den
Elementen geteilt, die den gemeinsamen Knoten an der Quelle haben. DieQuelle wird unter Berucksichtigung des Verhaltnisses α/360◦ zugeteilt,
wobei α der an der Quelle liegende Winkel des Elementes ist. Daher istdie korrekte Beziehung fur das Element (e) im Bild 6.4
{
f(e)Q
}
=αQ∗
360
01
0
. (6.40)
Es ist nicht notig den Winkel α fur jedes Element um die Quelle (oder
Senke) zu bestimmen. Wenn die Gleichungen durch direkte Steifigkeits-methode zusammengesetzt werden, summieren sich die Elementanteile auf
Q∗. Daher kann man einfach Q∗ auf die Reihe des Vektors {F} addieren,die der Knotennummer der Quelle entspricht. Eine Quelle hat einen posi-
tiven und eine Senke einen negativen Wert.
Beispiel 6.1
6.3. DREHFREIE STROMUNG IDEALER FLUSSIGKEIT 91
6.3 Drehfreie Stromung idealer Flussigkeit
Die drehfreie Stromung von idealen Flussigkeiten wurde intensiv studiert,
weil sie reichlich Informationen uber die Stromung um die Ecken, uberStauanlagen , durch Bauelementen und um die Tragflugel liefert. Idea-le drehfreie Stromung ist eine Approximation, die annimt, dass es keine
Reibung zwischen der Flussigkeit und den mit ihr in Kontakt stehendenFlachen gibt (ideal) und dass es keine Rotation oder Formanderung der
Flussigkeitspartikeln wahrend des Fließens gibt (drehfrei). Das Fließen desGrundwassers durch die Erde kann mit der Annahme der Rotationsfreiheit
sehr gut approximiert werden.
6.3.1 Stromung einer idealen Flussigkeit
Die zweidimensionale Stromung einer idealen Flussigkeit kann in Termen
von Stromfunktion ψ oder in Termen von Geschwindigkeitspotential φ for-muliert werden. Linien der konstanten ψ stehen senkrecht zu den Linien
der konstanten φ und die herrschenden Dgln. sind identisch und wie folgt
∂2ψ
∂x2+∂2ψ
∂y2= 0 und
∂2φ
∂x2+∂2φ
∂y2= 0. (6.41)
Die Randbedingungen fur ψ und φ sind jedoch nicht gleich und fuhren zu
unterschiedlichen Ergebnissen, die hier diskutiert werden.
Stromlinien-Formulierung
Eine Linie, bei der alle Punkten einen konstanten ψ-Wert besitzen, nennt
man eine Stromlinie. Die Stromungsrate, Qij, zwischen jedem Paar vonStromlinien ist gleich der Differenz in ihren ψ-Werten. Zwischen denStromlinien i und j ist die Stromungsrate
Qij = ψi − ψj. (6.42)
Senkrecht zu einer Stromlinie gibt es keine Stromung. Die Stromungsge-schwindigkeit an jedem Punkt ist ein Vektor, der tangential zur Stromlinie
an diesem Punkt lauft und seine Komponenten mittels
Vx =∂ψ
∂yund Vy =
∂ψ
∂x(6.43)
92KAPITEL 6. 2D-FELDPROBLEME MIT GEMISCHTEN RANDBEDINGUNGEN
4,5
cm
xV =5 cm/Sek
12 c
m
Abb. 6.5: Ebene Stromung umeinen Zylinder.
gerechnet werden. Die Annahme von einer idealen Flussigkeit impliziert,
dass die Flussigkeit nicht in ihre Umgebung dringt und sich auch nicht vonihren Kontakflachen trennt und damit keine freie Raume entstehen lasst.
Diese Bedingungen haben es zur Folge, dass die Stromungsgeschwindigkeitnormal zur Kontaktflache mit der Geschwindigkeit der Kontaktflache inder gleichen Richtung gleich ist (Duncan et al., 1970). Dies bedeutet, dass
es keine Stromung senkrecht zu einem festen Rand gibt, und die normaleGeschwindigkeitskomponente zum Rand gleich Null ist. Feste Rander so-
wie die Symmetrisachse parallel zur Fließrichtung stellen Stromlinien dar,weil die Geschwindigkeitskomponente und damit die Stromung senkrecht
zu denen Null ist.Die Randbedingungen fur Stromlinien-Formulierung werden bezuglich desim Bild 6.5 dargestellten Problems diskutiert. Eine horizontale und eine
vertikale Linie, die durch den Zylindermittelpunkt laufen, bilden fur die-ses Problem zwei Symmetrieachsen. Daher wird ein Viertel des Problems
betrachtet (s. Bild 6.6). Randbedingungen fur Stromlinien-Formulierungbestehen aus bekannten ψ auf einem Teil des Randes und aus Ableitun-
gen von ψ auf dem Restteil. Auf dem linken Rand ist eine konstante und
∂ψ
∂x= 0
∂ψ
∂x= 0
ψ = 30
ψ = 0
ψ = 0S.A
S.A
Abb. 6.6: Randbedingungen furStromlinien-Stromung um einenZylinder.
6.3. DREHFREIE STROMUNG IDEALER FLUSSIGKEIT 93
gleichmaßige Vx vorgegeben und damit ist Vy = 0. Mit Vy = 0 erhalt man
aus (6.43) fur diesen Rand ∂ψ/∂x = 0. Die gleichen Randbedingungengelten fur den rechten vertikalen Rand, weil er eine Symmetrieachse ist,
und die Stromlinien mussen um sie symmetrisch sein.Die horizontale Symmetrieachse, der Zylinderrand sowie der obere Randformen jeweils eine Stromlinie (s. Bild 6.6 und argumentiere warum). Der
Stromlinie am unteren Rand hat den Wert ψ = 0, weil auf dieser Liniekeine Stromung stattfindet. Die obere Stromlinie kann jeden vom Null
verschiedenen Wert haben; der entsprechende Wert ist hier 30, weil dieStromungsrate fur die Einheitstiefe 30 [cm3/sek] betragt (die Halfte der
Stromung uber den Gesamtquerschnitt).
Potential-Formulierung
In Potentialformulierung stehen die Geschwindigkeitskomponenten im Zu-
sammenhang mit φ:
Vx =∂φ
∂xund Vy =
∂φ
∂y. (6.44)
Bezogen auf die Stromung um das Zylinder im Bild 6.5 erhalt man fur den
linken Rand ∂φ/∂x = 5 und entlang der horizontalen Symmetrieachse unddem oberen Rand ∂φ/∂y = 0, weil dort keine vertikalen Geschwindig-
keitskomponenten vorhanden sind. Da die Geschwindigkeit normal zumZylinderrand Null ist, gilt am Zylinderrand ∂φ/∂n = 0. Dass die Potenti-
allilien senkrecht zu den Stromlinien stehen, nutzen wir zur Bestimmung
0∂φ
∂n= 0
∂φ
∂x= 5
5∂φ
∂y= 0
5∂φ
∂y= 0
φ = 50
S.A
S.A
Abb. 6.7: Randbedingungen furPotential-Stromung um einen Zy-linder.
94KAPITEL 6. 2D-FELDPROBLEME MIT GEMISCHTEN RANDBEDINGUNGEN
der letzten Randbedingung auf der rechten Seite. Die rechte vertikale Kan-
te muss eine Potentiallinie sein, weil sie eine Symmetrieachse ist und alleStromlinien senkrecht dazu verlaufen. Dieser Kante wird den beliebigen
Wert φ = 50 zugewiesen. Der zugewiesene Wert hat keinen Einfluss aufdas Ergebnis, weil die Komponenten der Geschwindigkeit allein von denGradienten von φ abhangen und nicht von seinem absoluten Wert. Diese
Randbedingungen sind in Bild 6.7 dargestellt.Der quantitative Wert der Randbedingung entlang der linken Kante wird
durch den Vergleich zwischen der aktuellen und der theortischen Randbe-dingung bestimmt. Die aktuelle Randbedingung ist
∂φ
∂x= 5 (6.45)
und die theoretische Randbedingung nach Gl. (6.1) und Bild 6.2 lautet
−Dx∂φ
∂x= S oder
∂φ
∂x= −S, (6.46)
wenn Dx = 1 ist. Der Vergleich zwischen (6.46) und (6.45) ergibt S = −5.Die Randbedingung φ = 50 auf der rechten Seite im Bild 6.7 ist ein
Sonderfall. Wenn der umstromte Korper eine irregulare Form hat, mussdas FE-Netz soweit in Stromabwarts gestreckt sein, dass die Stromungeine konstante Geschwindigkeit erreicht (s. Bild 6.8). An solchen Stellen
sind die konstanten Potentiallinien zur Fließrichtung senkrecht und φ darfein beliebiger Wert zugewiesen werden.
0∂φ
∂n= 0
5∂φ
∂y= 0
5∂φ
∂y= 0
∂φ
∂x= Vx φ=konst.
Abb. 6.8: Randbedingungen fur Potential-Stromung um eine irregulare Form.
6.3. DREHFREIE STROMUNG IDEALER FLUSSIGKEIT 95
6.3.2 Grundwasserstromung
Es gibt zwei wichtige Grundwasserprobleme, die mit der zweidimensio-
nalen Feldgleichung beschrieben werden. Das erste Problem ist das Ver-sickern vom Wasser unter einem Staudamm und wird mit folgender Glei-
chung
Dx∂2φ
∂x2+Dy
∂2φ
∂y2= 0 (6.47)
beschrieben. Dx und Dy sind in dieser Gleichung Durchlassigkeitskoef-fizienten [m/Tag] und φ steht fur die Wassersaulenhohe [m], die vondem Boden des Wasserspeichers gerechnet wird. Die Randbedingungen
bestehen in der Regel aus bekannten Werten von φ innerhalb des Wassersund verschwindendem Durchsickern auf anderen Rander (s. Bild 6.9). Ei-
ne undurchlassige vertikale Wand unter der Staumauer wird durch einenschmalen Spalt modeliert. Die FE-Methode erzeugt automatisch eine un-
durchlassige Randbedingung, ∂φ/∂n = 0, auf beiden Seiten des Spaltes,wenn dort keine expliziten Randbedingung vorgegeben wird.
Das zweite Grundwasserproblem besteht in der Bestimmung von Absen-
kung des Wasserspiegels wahrend des Pumpens im Boden rund um einenBrunnen und fern davon. Die folgende Gleichung ist die herrschende Feld-
gleichung fur einen beschrankten Grundwasserspeicher mit einem Zapf-
Dx
Dy
φ = bφ = a
0∂φ
∂n= 0
0∂φ
∂n= 0
∂φ
∂n
ba
Abb. 6.9: Randbedingungen fur das Versickern vom Wasser unter einem Staudamm.
96KAPITEL 6. 2D-FELDPROBLEME MIT GEMISCHTEN RANDBEDINGUNGEN
brunnen
Dx∂2φ
∂x2+Dy
∂2φ
∂y2+Q = 0, (6.48)
wobei Dx, Dy und φ die gleichen Bedeutungen wie bei der Gl. (6.47) ha-
ben. Der Term Q steht fur eine Punktsenke, nammlich den Brunnen, undmuss mittels des diskutierten Konzeptes in 6.2 ausgewertet werden. Das
beste Ergebnis wird naturlich dann erreicht, wenn die Senke sich direktauf einem Knoten befindet. Die mit der Gl. (6.48) assoziierten Randbe-
dingungen bestehen aus bekannten φ-Werten auf einen Teil des Randesund aus dem Versickern des Wassers in den Speicher entlang des restlichenRandes. Das Versickern ist durch folgende ableitungsenthaltene Randbe-
dingung definiert
−
(
Dx∂φ
∂xcos θ +Dy
∂φ
∂ysin θ
)
= S (6.49)
Die Geschwindigkeitskomponenten des Wassers werden mittels Darcy-Gleichungen gerechnet
Vx = −Dx∂φ
∂x
Vy = −Dy∂φ
∂y(6.50)
6.4 Warmeubertragung durch Leitung und Strah-
lung
Eine der fruheren Anwendungen von FE-Methode außerhalb der Struk-turmechanik ist in den Problemen der Warmeubertragung durch Leitungund Strahlung gewesen (Visser, 1965). Losung von Warmeubertragungs-
problemen mittels FE-Methode ist besonders bei den jenigen beliebt, diesich mit der Analyse von thermischen Spannungen beschaftigen, weil die
Losung des Warmeubertragungsproblems als Eingabe fur die Spannungs-analyse dient und das gleiche FE-Netz fur beide Aufgaben verwendet wer-
den kann.Wir diskutieren in diesem Abschnitt die Losung von vier unterschiedli-
chen Warmeubertragungsproblemen, wovon zwei die Warmeubertragung
6.4. WARMEUBERTRAGUNG DURCH LEITUNG UND STRAHLUNG 97
bei einer Schaufel behandeln. Bei dem dritten geht es um eine Verbund-
mauer, und das vierte ist die Analyse eines klassischen zweidimensionalenProblem mit konvektiven Randbedingungen.
6.4.1 Eindimensionale Schaufel
Die herrschende Dgl. fur stationare Warmeubertragung von einer eindi-
mensionalen Schaufel ist
kAd2φ
dx2− hPφ + hPφf = 0, (6.51)
wobei k die thermische Leitfahigkeit, h den Konvektionskoeffizient, A dieQuerschnittsflache, P der Querschnittsumfang von Schaufel und φ die
Temperatur ist. Alle Punkte auf einem Querschnitt mit einem bestimmtenKoordinate x haben den gleichen Temperaturwert. Die mit der Gl. (6.51)
assoziierte Randbedingung ist normalerweise eine bestimmte Temperaturan der Stelle x = 0
φ(0) = φ0 (6.52)
und der konvektive Warmeverlust am freien Ende der Schaufel
−kAdφ
dx= hA(φb − φf) an x = H, (6.53)
wobei φb die Temperatur am freien Schaufelende und vorerst unbekannt
ist. Der Konvektionskoeffizient in (6.53) kann der gleiche sein wie in (6.51)oder auch nicht.Die herrschende Dgl. (6.51) hat die allgemeine Form
Dd2x
dx2−Gφ+Q = 0, (6.54)
bei der D = kA, G = hP und Q = hPφf ist. Der Elementbeitrag zu denGalerkin’schen Restgleichung {R(e)} ist
{
R(e)}
= −
∫ xj
xi
[N ]T(
Dd2φ
dx2−Gφ+Q
)
dx
= −
∫ xj
xi
[N ]T(
Dd2φ
dx2+Q
)
dx+
∫ xj
xi
G[N ]Tφdx (6.55)
98KAPITEL 6. 2D-FELDPROBLEME MIT GEMISCHTEN RANDBEDINGUNGEN
Das erste Integral in (6.55) haben wir in den vorigen Kapiteln studiert.
Es ist gleichwertig mit {I(e)} + [k(e)]{Φ} − {f (e)}. Das zweite Integral in(6.55) ist eine neue Große, die auszuwerten ist. Mit Berucksichtigung von
φ(e) = [N ]{Φ(e)} und dem Einsatz ins Integral erhalten wir∫ xj
xi
G[N ]Tφdx =
(∫ xj
xi
G[N ]T [N ]dx
)
{
Φ(e)}
(6.56)
Da dieses Integral mit {Φ(e)} multipliziert wird, ist es einen Teil von Ele-
mentsteifigkeitsmatrix. Fuhren wir die Definition[
k(e)G
]
=
∫ xj
xi
G[N ]T [N ]dx, (6.57)
ein, dann erhalten wir{
R(e)}
={
I(e)}
+([
k(e)D
]
+[
k(e)G
]){
Φ(e)}
−{
f(e)Q
}
. (6.58)
wobei [k(e)D ] in (2.57) und {f
(e)Q } in (2.58) definiert wurde.
Das Integral in (6.57) wird am einfachsten in s- oder ℓ1, ℓ2-Koordinatenausgewertet. Das Auswertungsergebnis ist (vom Leser zu zeigen)
[
k(e)G
]
=GL
6
[
2 11 2
]
. (6.59)
Die Formulierung mittels Interelementvektor {I(e)} beinhaltet die ablei-tungsenthaltene Randbedingung, die durch (6.53) definiert wurde. Dieser
Vektor ist nach (2.56)
{
I(e)}
=
Ddφdx
∣
∣
∣
x=Xi
−Ddφdx
∣
∣
∣
x=Xj
(6.60)
und kann in zwei Komponenten aufgeteilt werden:
{
I(e)}
=
{
Ddφdx
∣
∣
∣
x=Xi
0
}
+
{
0
−Ddφdx
∣
∣
∣
x=Xj
}
, (6.61)
dass mit dem Ausdruck{
I(e)}
={
I(e)i
}
+{
I(e)b
}
(6.62)
6.4. WARMEUBERTRAGUNG DURCH LEITUNG UND STRAHLUNG 99
gleichwertig ist, wenn {I(e)i } die interelementare Anforderung und {I
(e)b }
mit der Randbedingung assoziiert ist. Der vom Null verschiedene Term in
{I(e)b } ist allerdings gleich der linken Seite von (6.53). Daher
{
I(e)b
}
=
{
0hA(φb − φf)
}
=
{
0hAΦj
}
−
{
0hAφf
}
, (6.63)
weil φb und Φj dasselbe sind. Die Gl. (6.63) ist mit
{
I(e)b
}
=
[
0 0
0 hA
]{
Φi
Φj
}
−
{
0
hAφf
}
=[
k(e)M
]{
Φ(e)}
−{
f(e)S
}
(6.64)
gleichwertig, wobei gilt
[
k(e)M
]
=
[
0 00 hA
]
,{
f(e)S
}
=
{
0hAφf
}
. (6.65)
Die komplette Restgleichung erhalt man durch Einsetzen von {I(e)} in(6.58):{
R(e)}
={
I(e)i
}
+([
k(e)D
]
+[
k(e)G
]
+[
k(e)M
]){
Φ(e)}
−{
f(e)Q
}
−{
f(e)S
}
(6.66)
Mit dem Verzicht auf interelementare Forderung {I(e)i } erhalt man
{
R(e)}
=[
k(e)]{
Φ(e)}
−{
f (e)}
(6.67)
Ein Beitrag von [k(e)M ] zu [k(e)] gibt es nur fur das letzte Element und nur
dann, wenn fur das Schaufelende h 6= 0 ist. Fur ein isoliertes Ende gilt
[k(e)M ] = 0.
Beispiel 6.2
6.4.2 Verbundwand
Die herrschende Dgl. fur die Warmeubertragung durch eine Verbundwandist
kAd2φ
dx2= 0, (6.68)
100KAPITEL 6. 2D-FELDPROBLEME MIT GEMISCHTEN RANDBEDINGUNGEN
wobei entweder φ auf einer Seite oder auf beiden Seiten der Wand bekannt
ist, oder der konvektive Warmeverlust von einer Seite oder von beidenSeiten bekannt ist. Die konvektive Randbedingungen sind
kAdφ
dx= hA(φb − φf) am x = 0 (6.69)
und
−kAdφ
dx= hA(φb − φf) am x = H (6.70)
Die Elementsteifigkeitsmatrix ist durch
[
k(e)]
=kA
L
[
1 −1−1 1
]
+
[
hAi 00 0
]
+
[
0 00 hAj
]
=[
k(e)D
]
+[
k(e)Mi
]
+[
k(e)Mj
]
(6.71)
gegeben, wobei die zweite Matrix, [k(e)Mi
], aus der konvektiven Randbedin-
gung am Knoten i resultiert und die dritte Matrix, [k(e)Mj
], aus der konvek-tiveb Randbedingung am Konten j resultiert und mit der Matrix in (6.65)identisch ist. Der Elementkraftevektor ergibt sich als
{
f (e)}
=
{
hAiφf
0
}
+
{
0hAiφf
}
, (6.72)
wobei der erste Vektor von (6.69) und der zweite von (6.70) resultiert.Wenn die Temperatur auf der beiden Seiten der Wand gegeben ist, werden
die Großen {f (e)}, [k(e)Mi
] und [k(e)Mj
] vernachlassigt.
Beispiel 6.3
6.4.3 Zweidimensionale Schaufel
Eine zweidimensionale Schaufel ist eine dunne Metalplatte, die an ei-
nem Heißwasser- oder Heißdampfrohr, wie im Bild 6.10, appliziert wurde.Warme wird durch Leitung vom Rohr in die Schaufel und von der Schaufel
durch Strahlung (Konvektion) zum Umgebungsmedium (oft Luft) ubert-ragen. Der konvektive Warmeverlust geschieht uber beide Seiten und die
Kante der Schaufel. Da die Kante im Vergleich mit den Flachen sehr klein
6.4. WARMEUBERTRAGUNG DURCH LEITUNG UND STRAHLUNG 101
ist, wird der Warmeverlust uber die Kante vernachlassigt.
Die herrschende Dgl. fur die zweidimensionale Schaufel ist durch
kxt∂2φ
∂x2+ kyt
∂2φ
∂y2− 2hφ+ 2hφf = 0 (6.73)
gegeben, wobei kx und ky die thermischen Leitfahigkeitskoeffizienten in x-
und y-Richtung sind, t die Dicke der Schaufelplatte angibt, h fur die Kon-vektionskoeffizient steht und φf die Umgebungstemperatur von Schaufelist. Die Randbedingungen sind: Entlang des Rohrrandes gilt
φ(Γ) = φs (6.74)
und entlang der Schaufelkante gilt
kxt∂φ
∂xcos θ + kyt
∂φ
∂ysin θ = 0. (6.75)
Diese letzte Randbedingung ist eine isolierte Randbedingung. DieWarmeubertragung in dieser Schaufel ist ein zweidimensionales Problem,
weil die kleine Dicke der Schaufel lasst in Dickenrichtung keine Tempera-turgradienten entstehen.Die herrschende Dgl. (6.73) hat eine Form identisch mit der Dgl. (5.1),
wenn es gilt:
Dx = tkx, Dy = tky, G = 2h, Q = 2hφf . (6.76)
Die Elementmatrizen, die im Kapitel 5 entwickelt wurden, sind hier ohne
Modifikation anwendbar (sehen Sie dazu bitte nochmal das Beispiel imAbschnitt 5.3).
φ(Γ) bekannt
kxt∂φ
∂xcos θ + kyt
∂φ
∂ysin θ = 0
t
x
y
z Abb. 6.10: ZweidimensionaleSchaufel und ihre Randbedingun-gen.
102KAPITEL 6. 2D-FELDPROBLEME MIT GEMISCHTEN RANDBEDINGUNGEN
konvektiveWärmeübertragung
vorgegebenTemperatur
Γ1
Γ3
Γ2
Wärmefluss
Oberflächedurch
Abb. 6.11: Randbedingungsar-ten fur klassische Warmeubertra-gungsprobleme.
6.4.4 Lange zweidimensionale Korper
Eine andere Form von zweidimensionaler Warmeubertragung ist der Fall,in dem ein langer Korper mit konstantem Querschnitt entlang seiner ge-
samten Lange die gleichen Randbedingungen hat. In dieser Situation istder Temperaturgradient in der Langsrichtung des Korpers (z-Koorinate)
gleich Null und die herrschende Dgl. ist durch
kx∂2φ
∂x2+ ky
∂2φ
∂y2+Q = 0 (6.77)
gegeben, wobei kx und ky die thermischen Leitfahigkeitskoeffizienten sind
und Q eine innere Warmequelle oder -senke darstellt. Die innere Warme-quelle oder -senke muss entlang der gesamten Korperlange (z-Richtung)
existieren, um die Zweidimensionalitat gewahrleistet zu sein. Das Problemist damit zweidimensional, weil die Temperatur sich nur in xy-Ebene va-
riieren kann.Die Dgl. (6.77) ist mit folgenden Beziehungen bereits in der Dgl. (5.1)
enthalten
Dx = kx, Dy = ky, G = 0, Q = Q. (6.78)
und die im Kapitel 5.3 und 5.4 entwickelten Elementmatrizen sind ohneModifikation anwendbar.
Die Warmeubertragung in einem langen Korper unterscheidet sich von derin einer dunnen Schaufel in den Randbedingung. Die moglichen Randbe-
dingungen sind, vorgegebene Temperaturwerte, konvektive Warmeubert-ragung und Warmefluss durch Oberflachen (s. Bild 6.11). Die letzten
zwei Randbedingungen wurden in ihren allgemeinen Termen bereits im
6.4. WARMEUBERTRAGUNG DURCH LEITUNG UND STRAHLUNG 103
qc
∂φ
∂n
n
qc
∂φ
∂n
n
qh = hA(φb − φf) qh = hA(φf − φb)
(b)(a)
Abb. 6.12: Konvektive Warmeubertragung auf der Oberflache.
Abschnitt 5 diskutiert. Nun wollen wir diese Bedingungen in Bezug auf
Warmeubertragung betrachten und nehmen dabei an, dass kx = ky = kist. Diese Annahme andert die ableitungsenthaltene Randbedingung (6.1)zu
k∂φ
∂n= −Mφb + S (6.79)
Konvektive Randbedingung
Das Thema ist hier die Bestimmug von Parametern, M und S in (6.79)und ihre Vorzeichen fur konvektive Randbedingung. Es sind zwei Situa-
tionen zu beachten: Warme, die den Korper verlasst und Warme, die inden Korper fließt. Betrachte die Konvektionsbedingung in Bild 6.12a, wo-
bei der Korper Warme verliert. In diesem Fall ist der Temperaturgradient∂φ/∂n negativ und die zur Oberflache geleitete Warme ist:
qc = −kA∂φ
∂n. (6.80)
Das negative Vorzeichen bringt, dass die Warmemenge positiv bleibt. Die
zur Oberflache geleitete Warme muss der Warmemenge entsprechen, diedie Oberflache durch Konvektion verliert; daher
−kA∂φ
∂n= hA(φb − φf). (6.81)
104KAPITEL 6. 2D-FELDPROBLEME MIT GEMISCHTEN RANDBEDINGUNGEN
Das Umschreiben von (6.81) in die Form von (6.79) ergibt
k∂φ
∂n= −hφb + hφf (6.82)
und daraus schließen wir:
M = h und S = hφf (6.83)
Nun betrachten wir den Fall, wenn die Warme in den Korper fließt (s.
Bild 6.12b). Der Temperaturgradient ∂φ/∂n ist dann positiv und dahergilt
qc = kA∂φ
∂n. (6.84)
Ein negatives Vorzeichen ist hier nicht mehr notig, weil der gradient positivund damit auch die Warmemenge qc positiv ist. Durch Gleichsetzen von
dieser Warme mit der, die auf die Oberflache antrifft, erhalt man
kA∂φ
∂n= hA(φf − φb). (6.85)
und
k∂φ
∂n= −hφb + hφf (6.86)
Daraus schließen wir nochmal
M = h und S = hφf . (6.87)
Diese Analyse zeigt, dass M = h und S = hφf ist, unabhangig davon,
ob der Korper durch Konvektion Warme verliert oder gewinnt. Diese Er-gebnisse gelten sowohl fur die außere als auch fur innere Flachen wie bei
einem Schornstein. Eine positive Normale, n, zeigt bei diesen Analysenstets von der Flache nach außen.
Warmefluss in den Korper
Bild 6.13 zeigt schematisch die Situation, in der Warme auf einen Teil des
Randes aufgebracht wird. Diese Warme muss von dem Rand weggeleitetwerden. Daher hat man
qc = kA∂φ
∂n. (6.88)
6.4. WARMEUBERTRAGUNG DURCH LEITUNG UND STRAHLUNG 105
qc
q*A
∂φ
∂n
n
Abb. 6.13: Warmefluss durcheine Oberflache.
Die auf die Flache aufgebrachte Warme ist gleich q∗A und durch Gleich-setzen der Warmeflusse erhalt man
kA∂φ
∂n= q∗A. (6.89)
Das Gleichsetzen von (6.89) und (6.79) ergibt
M = 0 und S = q∗, (6.90)
wobei q∗ den Warmefluss pro Einheitsflache darstellt. Wenn Warme vom
Korper abgefuhrt wird, dann gilt S = −q∗. Die Schlussfolgerung aus die-ser Analyse ist, dass Warmefluss auf dem Rand mit einem Vorzeichen
assoziiert ist. Warme gilt als positiv, wenn sie in den Korper fließt undals negative, wenn sie vom Korper abgefuhrt wird. Wenn am Rand einWarmefluss statfindet, ist der Koeffizient M gleich Null.
Abschließende Bemerkungen
Konvektive Warmeubertragung und Warmefluss-Randbedingung werden
durch Elementmatrizen, die im Kapitel 5 entwickelt wurden, in die FE-Analyse eingebracht. Im Falle der konvektiven Warmeubertragung gilt:
M = h und S = hφf und im Falle der Warmefluss-Randbedingung:M = 0und S = ±q∗. Warme gilt sowohl fur q∗ als auch fur Q als positiv, wenn
sie in den Korper hineingeht.
106KAPITEL 6. 2D-FELDPROBLEME MIT GEMISCHTEN RANDBEDINGUNGEN
6.5 Akustische Vibrationen
Wenn in der stationaren Feldgleichung (5.1) der Wert von Q gleich Nullund der von G negativ ist, erhalt man die Helmholz’sche Gleichung
Dx∂2φ
∂x2+Dy
∂2φ
∂y2+Gφ = 0 (6.91)
Physikalische Probleme, die mit dieser Dgl. beschrieben werden, sind die
Wellenausbreitung in flachen Gewasser und akustischen Vibrationen ingeschlossenen Raumen.
Die Losung von Helmholz’scher Gleichung bedingt die Losung eines Eigen-wertproblems, weil die Randbedingungen der Art sind, dass der globale
Kraftevektor, {F}, gleich Null ist. Die Losung der Helmholz’schen Glei-chung ist das Thema dieses Abschnittes und die akustische Vibration stellt
das Problemgebiet dar, mit dessen Hilfe die Losungstechniken diskutiertwerden. Wir fangen mit eindimensionalem Probelm an, weil die Berech-nungen ohne Rechner durchgefuhrt werden konnen. Die naturlichen Fre-
quenzen und die Schwingungsformen eines zweidimensionalen Problemsrunden diesen Abschnitt ab.
6.5.1 Eindimensionale Vibrationen
Die herrschende Dgl. fur ein Druckfeld, das mit den akustischen Vibratio-nen in ein einem zweidimensionalen Raum mit festen Randern assoziiert,
ist durch∂2φ
∂x2+∂2φ
∂y2+w2
c2φ = 0 (6.92)
gegeben, wobei φ der Druckunterschied zur Umgebung, w die Wellenfre-quenz und c die Wellengeschwindigkeit im betrachteten Medium ist. Die
Randbedingung, die an jedem Rand erfullt sein muss, ist
∂φ
∂n= 0. (6.93)
Das eindimensionale Anlogon zu (6.92) lautet
d2φ
dx2+w2
c2φ = 0 (6.94)
6.5. AKUSTISCHE VIBRATIONEN 107
mit dφ/dx = 0 an jedem Ende. Die Gl. (6.94) hat die allgemeine Form
Dd2φ
dx2−Gφ = 0 (6.95)
mit D = 1 und G = −w2/c2. Die Elementsteifigkeitsmatrix fur den ersten
Term, d2φ/dx2, ist durch die Gl. (2.57) gegeben und der Beitrag vom Term−Gφ dazu durch die Gl. (6.59). Die komplette Elemetsteifigkeitsmatrix fur
(6.95) ergibt sich zu
[
k(e)]
=1
L
[
1 −1−1 1
]
−w2L
6c2
[
2 11 2
]
. (6.96)
Der Elementkraftevektor, {f (e)}, beinhaltet nur Nullen, weil es keinen
Quellenterm in (6.95) gibt und die Randbedingung dφ/dx = 0 an beidenEnden keinen vom Null unterscheidenen Term in {f (e)} erzeugt.
Das zu untersuchende Problem ist ein geschlossenes Rohr wie im Bild ??
dargestellt. Das FE-Modell hat zwei Elemente der gleichen Lange H/2.
Einsetzen von H/2 fur L in (6.96) und Multiplizieren des Gesamten mitH/2 ergibt
[
k(e)]
=
[
1 −1
−1 1
]
− Z
[
2 1
1 2
]
, (6.97)
wobei gilt:
Z =w2H2
24c2(6.98)
Das Kombinieren der zwei Elementsteifigkeitsmatrizen fuhrt zu einem Sy-
stem von drei Gleichungen:
1 −1 0
−1 2 −10 −1 1
Φ1
Φ2
Φ3
− Z
2 1 0
1 4 10 1 2
Φ1
Φ2
Φ3
=
0
00
(6.99)
oder in kompakter Form
([KD] − Z [KG]) {Φ} = {0} (6.100)
Beide Matrizen, [KD] und [KG], sind symmetrisch. Wahrend [KD] einepositiv-definierte Matrix ist, ist [KG] eine semidifinierte Matrix (Determi-
nante gleich Null). Nach Eigenwert-Theorie sind alle Eigenwerte, Zi, die
108KAPITEL 6. 2D-FELDPROBLEME MIT GEMISCHTEN RANDBEDINGUNGEN
die Gl. (6.99) erfullen, reele, positive und getrennte Zahlenwerte und die
entsprechenden Eigenvektoren, {Φ}i, sind voneinander unabhangig.Die Eigenwerte, Zi, sind die Werte von Z, mit denen die Determinante von
(6.99) Null wird. Die Kombination der beiden Matrizen in (6.99) ergibt
(1 − 2Z) −(1 + Z) 0−(1 + Z) (2 − 4Z) −(1 + Z)
0 −(1 + Z) (1 − 2Z)
Φ1
Φ2
Φ3
=
00
0
(6.101)
Die Determinante der Koeffizientenmatrix ist
2(1 − 2Z)[
(1 − 2Z)2 − (1 + Z)2]
= 0, (6.102)
die die folgenden Nullstellen (Eigenwerte) hat:
Z1 = 0, Z2 =1
2und Z3 = 2 (6.103)
Es gibt zu jeder Nullstelle in (6.103) einen assoziierten Eigenvektor {Φ}i.
Es ist unmoglich die drei Komponenten von {Φ}i eindeutig zu bestim-men, weil das resultierende Gleichungssytem homogen ist. Das ubliche
Verfahren hier ist die Zuweisung eines beliebigen Wertes zu einem derKomponenten und die Berechnung der weiteren zwei in Termen des er-
sten.Mit Z1 = 0 liefert die Gl.(6.101) das folgende Gleichungssystem
Φ1 − Φ2 = 0
−Φ1 + 2Φ2 − Φ3 = 0 (6.104)
−Φ2 + Φ3 = 0
Aus diesem Gleichungssystem resultiert: Φ1 = Φ2 = Φ3 und der ersteEigenvektor ergibt sich mit dem beliebigen Wert Φ1 = 1 zu
{Φ}T1 = [ 1 1 1 ] (6.105)
Mit Z2 = 12 in (6.101) erhalt man das Gleichungssystem
−3
2Φ2 = 0
−3
2Φ1 −
3
2Φ3 = 0 (6.106)
−3
2Φ2 = 0
6.5. AKUSTISCHE VIBRATIONEN 109
Daraus resultiert: Φ1 = 1, Φ2 = 0 und Φ3 = −1. Der zweite Eigenvektor
ergibt sich mit dem beliebigen Wert Φ1 = 1 zu
{Φ}T2 = [ 1 0 −1 ] (6.107)
Den dritten Eigenvektor bekommt man mit Z3 = 2 und Φ1 = 1 als
{Φ}T3 = [ 1 −1 1 ] (6.108)
Die theoretischen Werte fur die naturlichen Frequenzen, wn, sind durch
wn =nπc
H(6.109)
gegeben. Die rechnerischen Werte fur wn erhalt man durch Einsetzen der
Nullstellenwerte, Z1, Z2 und Z3, in die Gl. (6.98) und Losen fur wn. Diegerechneten Werte fur wn sind
w1 = 0, w2 =3, 464c
Hund w3 =
6, 928c
H,
die sich sehr gut mit den theoretischen Werten
w1 = 0, w2 =3, 142c
Hund w3 =
6, 283c
H,
vergleichen lassen, wenn das FE-Netz nur aus zwei Elementen besteht.Die theoretischen Schwingungsformen haben die allgemeine Form P =
cos(nπx/H). Die theoretischen Schwingungsformen und die berechnentenEigenvektoren sind in Bild ?? dargestellt. Fur die erste Schnigungsformsind die theoretischen und gerechneten Formen deckungsgleich.
6.5.2 Zweidimensionale Vibrationen
Die zweidimensionale akustische Schwingungsgleichung und ihre Randbe-dingungen sind durch (6.92) und (6.93) gegeben. Fur diese Gleichung sind
die Elemetmatrizen fur das Dreieckselement durch (5.37) und (5.39) undfurs Rechteckselement durch (5.50) und (5.56) gegeben.
Ein rechteckiger Raum mit der Lange 20 m und der Breite 10 m wirddurch vier Dreieckselemente diskretisiert (s. Bild ??). Mit der Definition
Z = w2/c2 erhalt man das globale Gleichungssystem als
([KD] − Z[KG]){Φ} = {0}, (6.110)
110KAPITEL 6. 2D-FELDPROBLEME MIT GEMISCHTEN RANDBEDINGUNGEN
wobei
[KD] =1
8
10 3 0 −3 −10
3 10 −3 0 −100 −3 10 3 −10
−3 0 3 10 −10
−10 −10 −10 −10 40
(6.111)
[KG] =1
6
100 25 0 25 5025 100 25 0 500 25 100 25 50
25 0 25 100 5050 50 50 50 200
(6.112)
und
{Φ}T = [ Φ1 Φ2 Φ3 Φ4 Φ5 ]. (6.113)
Handberechnungen der Z-Werte, die den Wert von Determinante, |([KD]−
Z[KG])|, gleich Null machen, ist nicht gerade sehr vernunftig; ein Rechen-programm muss hier eingesetzt werden. Diskussion uber direkten und it-
terativen Methoden fur die Auswertung der Eigenwerte gibt es in vielenTextbuchern. Bathe und Wilson (1976) presentieren eine ausfuhrliche Dis-
kussion uber Eigenwertberechnung in Bezug auf FE-Probleme. Fur die Gl.(??) sind die Eigenwerte und Eigenvektoren wie folgt:
Z1 = 0 {Φ}T1 = [ 1 1 1 1 1 ]
Z2 = 0, 030 {Φ}T2 = [ 1 −1 −1 1 0 ]
Z3 = 0, 120 {Φ}T3 = [ 1 1 −1 −1 0 ] (6.114)
Z4 = 0, 150 {Φ}T4 = [ 1 −1 1 −1 0 ]
Z5 = 0, 450 {Φ}T5 = [ −0, 5 −0, 5 −0, 5 −0, 5 1 ]
Die berechneten Eigenwerte lassen sich gut mit den theoretischen Ergeb-nissen aus Z = w2/c2 vergleichen. Die theoretischen Werte sind:
Z1 = 0 , Z2 = 0, 0247 , Z3 = 0, 0987 , Z4 = 0, 123 , Z5 = 0, 395
Jede Schwingungsform kann wie im eindimensionalen Fall grafisch darge-
stellt werden. Die Schwingungsformen, die {Φ}2 und {Φ}4 entsprechen,sind im Bild ?? dargestellt.