2.1. Differentialgleichungen der linearen Elastostatik 43
2. Grundlagen der FEM für die ElastomechanikIm folgenden Abschnitt wollen wir wichtige Grundgleichungen der FEM zur Ent-wicklung von finiten Elementen für Probleme der Elastomechanik bereitstellen. Dabeibeschränken wir uns hier auf linear-elastische Probleme der Statik und Dynamik, leitenaber die Grundgleichungen der FEM für den allgemeinen 3D-Fall ab. Dabei erfolgt dieHerleitung der Elementsteifigkeitsmatrix, der Elementmassenmatrix und des Element-kraftvektors für ein beliebiges finites Element, dessen Verformungszustand durch zu-nächst nicht näher spezifizierte Ansatzfunktionen approximiert wird. Auf die Assem-blierung des Elementgleichungssystems zum Gesamtgleichungssystem einer Vernet-zung, die Einarbeitung von Randbedingungen und die Lösung des Gleichungssystemswird nur kurz eingegangen.1
Von den Grundgleichungen ausgehend werden dann im Abschnitt 3 für einfachefinite Elemente - die sogenannte Simplexelemente2 - die Ableitung der Elementmatri-zen und der Elementkraftvektoren exemplarisch dargestellt.
2.1. Differentialgleichungen der linearen Elastostatik
In diesem Abschnitt werden einleitend die wichtigsten Gleichungen der linearen Ela-stostatik in Matrizenschreibweise angegeben.3 Probleme der Elastomechanik lassensich durch die folgenden drei Sätze von Gleichungen beschreiben:
♦ Zusammenhang zwischen Verrzerrungen und Verschiebungen♦ Gleichgewichtsbedingungen♦ Konstitutive Gleichungen (Materialgesetz)
Dazu kommen noch Gleichungen für die Randbedingungen in den Verschiebungenund in den Spannungen.
Zusammenhang zwischen Verzerrungen und Verschiebungen4
Der Zusammenhang zwischen Verzerrungen (Dehnungen) und Verschiebungen lautetin Matrizenschreibweise
uD=ε in V (2.1-1)
Dabei ist e der Vektor der Verzerrungen, D ist eine Differentiationsmatrix, die partielleAbleitungen nach den globalen Koordinaten enthält, und u ist der Verschiebungsvek-tor. Für ein-, zwei- und dreidimensionale Probleme wird die Gleichung (2.1-1) nach-folgend spezifiziert.
Eindimensionaler Fall (1D)
[ ] [ ]uxx
∂∂
==ε= Duε (2.1-2)
1 Siehe dazu auch die Vorlesung Mathematische und Numerische Methoden der Mechanik – Teil II – Computer-Numerik.2 Finite Elemente mit den einfachsten möglichen Ansatzfunktionen.3 Siehe dazu die auch die Vorlesungen Kontinuumsmechanik/Elastizitätstheorie und Mathematische und Nume-rische Methoden der Mechanik - Teil I.4 Manchmal auch als kinematische Verträglichkeitsbedingungen bezeichnet.
44 2. Grundlagen der FEM für die Elastomechanik
Zweidimensionaler Fall (2D)5
∂∂
∂∂
∂∂
∂∂
==
γεε
=vu
xy
y0
0x
xy
y
x
Duε (2.1-3)
Dreidimensionaler Fall (3D)
∂∂
∂∂
∂∂
∂∂∂∂
∂∂
∂∂
∂∂
∂∂
==
γγγεεε
=wvu
x0
z
yz0
0xy
z00
0y
0
00x
zx
yz
xy
z
y
x
Duε (2.1-4)
Verschiebungsrandbedingungen6
Wenn auf einem Teil des Randes Ou Verschiebungen u vorgeschrieben sind, mußgelten
uu = auf uO . (2.1-5)
Gleichgewichtsbedingungen7
Die Gleichgewichtsbedingungen lauten
0T =+ pD σ in V, (2.1-6)
wobei p den Vektor der Volumenlasten bezeichnet und s der Spannungsvektor ist, derim 1D-, 2D- und 3D-Fall folgende Komponenten hat.
Eindimensionaler Fall (1D)
[ ]xσ=?σ
5 Als 2D-Modell wird hier stets ein Scheibenmodell mit einem ebenen Spannungszustand vorausgesetzt. Auf an-dere, oft auch mit 2D bezeichnete Modelle (z.B. ebener Verformungszustand, rotationssymmetrische Problemeund Plattenprobleme) wird später noch eingegangen werden.6 Manchmal auch als kinematische Randbedingungen bezeichnet.7 Manchmal auch als statische Verträglichkeitsbedingungen bezeichnet.
2.1. Differentialgleichungen der linearen Elastostatik 45
Zweidimensionaler Fall (2D)
[ ]xyyx ss τ=?σ
Dreidimensionaler Fall (3D)
[ ]zxyzxyzyx τττσσσ=Τσ
Spannungsrandbedingungen8
Wenn auf einem Teil der Oberfläche Oq Flächenlasten wirken, dann muß auf jedemdifferentiellen Teil dieser Oberfläche Gleichgewicht mit den Spannungen an dieserStelle herrschen, d.h. es gilt
qsTq == Tn auf qO . (2.1-7)
Die Matrix Tn transformiert den Spannungsvektor σ auf die Randspannungen q (eng-lisch: tractions). Sie enthält die Richtungskosinus der Randnormale. Für den 3D-Falllautet diese Matrix beispielsweise
=
nxny
nznx
nzny
nz
ny
nxTn
cc00cc
c0c
c000c000c
T , (2.1-8)
wobei die cni =cos(n, xi), mit x1=x, x2=y, x3=z, die Randnormalen sind. Der Vektor qenthält die vorgegebenen (bekannten) Flächenlastkomponenten bezogen auf das glo-bale kartesische Bezugssystem.9
Konstitutive Gleichungen10
Wir wollen hier zunächst nur Materialien betrachten, die durch das klassische Hooke-sche Gesetz beschrieben werden können. Die in die Materialmatrix E eingehendenMaterialparameter müssen aus Experimenten bestimmt werden.
Ees = (2.1-9)
Dabei bezeichnet E die Hookesche Matrix, die bekanntlich im isotropen Fall durchzwei unabhängige Stoffgrößen, den Elastizitätsmodul11 E und die Querkontraktions-zahl12 ? beschrieben ist13.
Eindimensionaler Fall (1D)
xx Eε=σ (2.1-10)
8 Manchmal auch als statische Randbedingungen bezeichnet.9 Im Rahmen dieser Vorlesung werden wir als globales Bezugssystem stets ein kartesisches Koordinatensystemverwenden.10 Materialgesetz11 In der englischsprachigen Literatur meist als Youngscher Modul bezeichnet.12 Manchmal auch als Poissonsche Zahl bezeichnet.13 Teilweise werden statt E und ν auch die Lameschen Konstanten ? und µ verwendet:
)21)(1(
Eund
)1(2E
Gν−ν+
ν=λ
ν+==µ
46 2. Grundlagen der FEM für die Elastomechanik
Zweidimensionaler Fall (2D)Für den 2D-Fall gilt allgemein:
τεε
=
τσσ
xy
y
x
33
2322
131211
xy
y
x
E.symmEEEEE
(2.1-11)
Für den allgemeinen anisotropen Fall müssen also 6 Materialkenngrößen bestimmtwerden. Wenn das Material isotrop ist, gilt für die Eij
22211 1E
EEν−
== , 212 ?1E?
E−
= , GE33 = , 0EE 2313 == (2.1-12)
Im isotropen Fall ist das Materialgesetz durch zwei Konstanten, z.B. E und ? eindeutigbeschrieben, der Schubmodul G ist mit E und ? durch folgende Gleichung verknüpft:
)1(2E
Gν+
= (2.1-13)
Ist das Material orthotrop, hat es bezogen auf zwei senkrecht zueinander stehendeRichtungen, die durch die Indizes 1 und 2 gekennzeichnet sind14, unterschiedliche Ei-genschaften. Man erhält für diesen Fall
0EE,mE
E
,n)n1(
EE,
n1E
E,n1E
E
23131
33
221
1222
21
121122
21
111
===
ν−=
ν−ν
=ν−
= (2.1-14)
mit
12
1
2
1
GE
m,EE
n == (2.1-15)
Es ergeben sich also insgesamt 4 unabhängige Materialkonstanten E1, E2, ?21 und G12.Aus Symmetriegründen muß außerdem die Beziehung
1
12
2
21
EEν
=ν
(2.1-16)
gelten, wodurch 12ν durch 21ν ausgedrückt werden kann.
Dreidimensionaler Fall (3D)Im beliebig anisotropen Fall hat die symmetrische Elastizitätsmatrix Eij das Format6×6, und man benötigt 21 Materialkonstanten. Im isotropen Fall reichen wieder diezwei Größen E und ? aus, und es ergibt sich 14 Es wird hier angenommen, daß das globale Koordinatensystem mit den lokalen Materialrichtungen überein-stimmt. Anderenfalls ist eine Transformation der Hookeschen Matrix erforderlich. Wenn sich die Spannungenüber Tss = auf eine neues Koordinatensystem x transformieren, wobei T eine Orthogonalmatrix ist mit
T1 TT =− , ergibt sich für die Transformation der Hookeschen Matrix ETTE T= .
2.1. Differentialgleichungen der linearen Elastostatik 47
ν−ν−
ν−ν−
ν−ν−
−νν−νν
−νν
ν−ν+ν−
)1(221
0)1(2
21
00)1(2
210001
0001
1
00011
1
)21)(1()1(E
symm.
=E (2.1-17)
Wenn das Materialgesetz orthotrop bezüglich dreier senkrecht aufeinander stehenderKoordinatenrichtungen x1, x2 und x3 ist, kann die Hookesche Matrix15 in folgenderForm geschrieben werden:
γγγεεε
νν−νν+ννν−νν+ννν+ννν−
=
σσσσσσ
31
23
12
33
22
11
31
23
12
21123
132123231132
2312131321312132231
31
23
12
33
22
11
DGsymm.0DG00DG000)1(E000)(E)1(E000)(E)(E)1(E
D1
(2.1-18)
wobei hier D die folgende Abkürzung ist:
312312311332232112 21D ννν−νν−νν−νν−= (2.1-19)
Auch hier gelten wieder die Symmetriebedingungen
j
ji
i
ij
EE
ν=
ν , (2.1-20)
wobei ?ij die Querkontraktion in Richtung xj infolge einer Dehnung in Richtung xi dar-stellt. Es sind in diesem Fall also 9 unabhängige Materialkonstanten erforderlich undzwar drei Elastizitätsmoduln E1, E2, E3 in den drei Koordinatenrichtungen x1, x2 undx3, drei Querkontraktionszahlen ?12, ?23, ?13 und drei Schubmoduln G12, G23, G31. Na-türlich können diese Größen über die Symmetriebedingungen wieder durch andereGrößen ausgedrückt werden.
15 Die Ableitung erfolgt durch Aufstellen des Zusammenhanges zwischen den Dehnungen und den Spannungen inder Form
3131
31333
222
23111
1333
2323
23333
32222
111
1222
1212
12333
31222
21111
11
G1
E1
E1
E1
G1
E1
E1
E1
G1
E1
E1
E1
σ=γσ+σν−σν−=ε
σ=γσν−σ+σν−=ε
σ=γσν−σν−σ=ε
und die Inversion dieses Zusammenhanges.
48 2. Grundlagen der FEM für die Elastomechanik
2.2. Variationsformulierungen
Grundsätzlich besteht die Möglichkeit, eine Randwertaufgabe der Elastostatik mit Hil-fe von Differentialgleichungen oder mit Hilfe einer Variationsformulierung zu be-schreiben. Die Methode der finiten Elemente basiert üblicherweise auf einer Variati-onsformulierung des Problems. Falls keine Variationsformulierung bekannt ist, kannbeispielweise unter Verwendung der Methode der gewichteten Residuen (MWR) auchdirekt von den Differentialgleichungen ausgegangen werden.16 Wir betrachten hier zu-nächst nur die beiden klassischen Einfeldformulierungen und zwar:
1. Das Prinzip der virtuellen Verschiebungen (bzw. das Prinzip vom Minimum desElastischen Gesamtpotentials).
2. Das Prinzip der virtuellen Spannungen (Prinzip vom Minimum des Ergän-zungspotentials).17
1. Das Prinzip der virtuellen VerschiebungenFür einen Spannungszustand, der die statischen Verträglichkeits- und Randbedingun-gen (2.1-6) und (2.1-7) erfüllt, ist die bei einer virtuellen Verschiebung geleistete Ar-beit der inneren Kräfte IWδ gleich der virtuellen Arbeit der äußeren Kräfte AWδ .
0WW AI =δ−δ (2.2-1)
Das Prinzip kann in folgender Form geschrieben werden:
0dOdVdVqO
T
V
T
V
T =δ−δ−δ ∫∫∫ qupuse (2.2-2)
mit den Nebenbedingungen (2.1-1) und (2.1-5)
Due = in Vuu = auf O u
Damit ist (2.2-2) den Grundgleichungen (2.1-6) und (2.1-7) identisch. Zum Beweiswird (2.2-2) unter Verwendung des Gaußschen Integralsatzes umgeformt, und man er-hält:
( ) 0dOdO)(dVuq O
T
O
T
V
TT =δ+−δ++δ− ∫∫∫ quqqupsDu (2.2-3)
Wegen (2.1-5) müssen die Verschiebungsrandbedingungen von u erfüllt werden, sodaß der unterstrichene Term in (2.2-3) verschwindet ( 0u =δ ), und man erhält nachdem Fundamentallemma der Variationsrechnung als Eulersche Gleichungen des Va-riationsprinzips die Gleichgewichtsbedingungen (2.1-6) und außerdem die Spannungs-randbedingungen (2.1-7). Wird vorausgesetzt, daß sich die Spannungen und Belastun
16 Siehe dazu die Vorlesung Mathematische und Numerische Methoden der Mechanik - Teil I17 Eine Diskussion zu erweiterten Prinzipien und deren Anwendung in der FEM findet sich in dem Vorlesungss-kript U. Gabbert: Variationsprinzipe für die Methode der finiten Elemente.
2.2. Variationsformulierungen 49
gen aus Potentialen ableiten lassen,18 kann das Prinzip der virtuellen Verschiebungendurch die Variation des Gesamtpotentials (Elastisches Potential) p abgeleitet werden.
.MindOdVdV21
qO
T
V
T
V
T →−−=π ∫∫∫ qupuse (2.2-4)
Setzt man in Gleichung (2.2-4 die Gleichungen (2.1-1) und (2.1-9) ein, ergibt sich
( ) ( ) .MindOdVdV21
qO
T
V
T
V
T →−−=π ∫∫∫ qupuDuEDu (2.2-5)
0=πδmit den Nebenbedingungen Due = in V
uu = auf O u
Die Formulierung (2.2-5) besagt, daß für die exakte Lösung u, die die kinematischenRandbedingungen erfüllt, das Elastische Potential einen minimalen Wert annimmt.Diese Bedingung stellt in der Elastomechanik die übliche Ausgangsgleichung für dieAnwendung der Methode der finiten Elemente dar. Wir werden im nächsten Abschnittnur noch die Anfangsverzerrungen e0 in Gleichung (2.1-1) ergänzen, wie sie beispiels-weise durch Temperaturänderungen in der Form e0=a th? T entstehen. Dadurch erweitertsich Gleichung (2.1-1) zu e=Du-e0.
Über die Art des Extremwertes von p gibt die zweite Variation Aufschluß. Diese ist beipositiv definierter Matrix E (was für Hookesches Material stets gilt) größer als Null, sodaß der Extremwert von p ein absolutes Minimum ist.19
18 Die Existenz von Potentialen besagt, daß die Arbeit der inneren und äußeren Kräfte nicht vom Weg (d.h. derBelastungsgeschichte) abhängen, sondern nur vom Anfangs- und Endpunkt. Bezeichnet man beispielsweise dieauf das Volumen bezogene Arbeit der inneren Kräfte mit
Eeee T21
I )(W = ,
dann folgt aus der Wegunabhängigkeit auch die Existenz des totalen Differentials
eseEeee
ddd)W
(dW TTTII ==
∂∂
= ,
mit sEee
==∂
∂ Iw .
Siehe dazu beispielsweise G. Mehlhorn (Hrsg.): Der Ingenieurbau - Rechnerorientierte Baumechanik.Ernst&Sohn Verlag 1995, S.153ff.19 Aus (2.2-5) folgt
∫ ∫ ∫ δ−δ−δ=δπV V O
TTT
q
dOdVdV)()( qupuDuEuD ,
und daraus 0dVdV)()(V
T
V
T2 ≥δδ=δδ=πδ ∫∫ eEeuDEuD ,
wenn E eine positiv definite Matrix ist.
50 2. Grundlagen der FEM für die Elastomechanik
2. Das Prinzip der virtuellen SpannungenFür einen Verzerrungszustand, der die kinematischen Verträglichkeits- und Randbe-dingungen (2.1-1) und (2.1-5) erfüllt, ist die bei einer virtuellen Spannungsänderunggeleistete Arbeit der inneren Kräfte gleich der virtuellen Arbeit der äußeren Kräfte(2.2-1).
Wenn man beachtet, daß die fest vorgegebenen Größen (Volumen- und Oberflächen-kräfte) nicht variiert werden, kann das Prinzip der virtuellen Spannungen folgender-maßen geschrieben werden:
0dOdVuO
T
V
T =δ−δ ∫∫ uqes (2.2-6)
mit den Nebenbedingungen (2.1-6) und (2.1-7)
0psD =+T in V qq = auf O q
Auch hier wird das Materialgesetz a priori als gegeben vorausgesetzt. Damit ist (2.2-6)den Grundgleichungen (2.1-1) und (2.1-5) äquivalent. Umformen von (2.2-6) unterVerwendung des GAUßschen Integralsatzes liefert
0dV)(dOdO)(dV)(V
TT
O
T
O
T
V
T
qu
=δ−δ+−δ−−δ ∫∫∫∫ usDuquuqDues
Wegen (2.1-6) und (2.1-7) verschwinden die unterstrichenen Ausdrücke und man er-hält nach dem Fundamentallemma der Variationsrechnung als Eulersche Gleichungendes Variationsprinzips die kinematischen Gleichgewichtsbedingungen und die kine-matischen Randbedingungen. Es sei wieder vorausgesetzt, daß sich die Spannungenund Belastungen aus Potentialen ableiten lassen. Dann kann das Prinzip der virtuellenSpannungen auch durch die Variation des konjugierten Gesamtpotentials (auch Ergän-zungspotential) p* ausgedrückt werden.
.MindOdV21
*uO
T
V
T →−=π ∫∫ uqes (2.2-7)
0* =πδ
mit den Nebenbedingungen 0psD =+T in V qq = auf O q
Über die Art des Extremwertes gibt die zweite Variation Aufschluß. Für den Gültig-keitsbereich des Hookeschen Gesetzes (E - positiv definit) ist die zweite Variationgrößer als Null, so daß der Extremwert von p* ein absolutes Minimum ist.
2.2. Variationsformulierungen 51
Bei der Betrachtung der beiden genannten Variationsprinzipien wird ihr dualer Zu-sammenhang deutlich: Die Zwangsbedingungen des einen Prinzips stellen die natürli-chen Bedingungen des anderen dar und umgekehrt.
In der FEM wird das Prinzip vom Minimum des Elastischen Potentials bevorzugt, dasich problemlos Näherungsfunktionen für die Verschiebungen konstruieren lassen, diedie Nebenbedingungen des Prinzips erfüllen. Es ist für den allgemeinen Fall schwieri-ger, Spannungsansätze zu finden, die sämtliche Gleichgewichtsbedingen erfüllen. Au-ßerdem ist dann die Bestimmung der Verschiebungen aus den Spannungen (Integrati-on) ein weiteres Problem.
52 2. Grundlagen der FEM für die Elastomechanik
2.3. Die Verschiebungsgrößenmethode
Das Prinzip vom Minimum des Elastischen Potentials
Die Ausgangsbasis für die Ableitung von finiten Elementen nach der Verschiebungs-größenmethode ist das Prinzip vom Minimum des Elastischen Potentials in der Formder Gleichung (2.2-4). Wir wollen, wie bereits erwähnt, Anfangsverzerrungen ergän-zen und erhalten20 aus
.MindOdVdV21
qO
T
V
T
V
T →−−=π ∫∫∫ qupuse (2.3-1)
mit den Nebenbedingungen
0eDue −= in V, (2.3-2)u u= auf O u , (2.3-3)
und dem Zusammenhang zwischen Spannungen und Verzerrungen
)( 0eDuEEes −== , (2.3-4)
die gewünschte Ausgangsgleichung. Wenn außerdem noch Einzelkräfte berücksichtigtwerden sollen,21 ist in Gleichung (2.3-1) die Arbeit der äußeren Lasten um einen Anteilzu erweitern, der sich aus der Kraft multipliziert mit der Verschiebung in der Kraf-trichtung ergibt. Wenn wir generell nur Kräfte in Richtung der Achsen des globalenkartesischen Koordinatensystems zulassen, was keine Einschränkung darstellt, da jedeKraft in diese Richtungen zerlegt werden kann, ergibt sich ein zusätzlicher Term derForm
∑=
Fn
1kk
Tk Fu (2.3-5)
Dabei bezeichnet k den Ort des Kraftangriffs, uk und kF sind der vollständige Ver-schiebungs- bzw. der gegebene Kraftvektor an dieser Stelle, und nF bezeichnet die An-zahl der Punkte, an denen Einzelkräfte angreifen. Damit ergibt sich aus (2.3-1)
.MindOdVdV21 F
q
n
1kk
Tk
O
T
V
T
V
T →−−−=π ∑∫∫∫=
Fuqupuse (2.3-6)
In der Verschiebungsgrößenmethode sind die Verschiebung die primären Unbekann-ten. Den Verschiebungszustand in einem Bauteil oder in einem beliebigen Körperwollen wir durch elementweise definierte Näherungsfunktionen approximieren. DieseFunktionen sind nur in jeweils einem Element definiert, außerhalb des Elementes sindsie Null. Unter Beachtung dieser Tatsache kann das elastische Gesamtpotential p als
20 Über die Beschreibung der Anfangsverzerrungen lassen sich hier auch problemlos Anfangsspannungen berück-sichtigen, so daß diese nicht gesondert behandelt werden.21 Verschieben und Einzelkräfte sind hier als verallgemeinerte Größen aufzufassen, d.h. es kann sich in be-stimmten Fällen (z.B. bei Platten oder Schalen) auch um Verdrehungen bzw. Momente handeln.
2.3. Die Verschiebungsgrößenmethode 53
eine Summe von elastischen Potentialen p(e) aller n(e) finiten Elemente eines Netzes ge-schrieben werden.
∑=
π=π)e(n
1e
)e( (2.3-7)
Damit kann auch die Minimierung von p, d.h. die Ableitung von p nach den Unbe-kannten des Verschiebungsansatzes – den Knotenverschiebungen - für jedes Elementeinzeln ausgeführt werden.
)e()e(n
1e
)e( n1e,00)e(
L==δπ→=δπ=δπ ∑=
(2.3-8)
Die Ergebnisse der elementweisen Minimierung von π müssen dann entsprechendGleichung (2.3-8) nur noch „addiert“ werden. Aus dieser Addition läßt sich leicht derAlgorithmus zum Aufbau des Gesamtgleichungssystems aus den Elementbeiträgen(Elementsteifigkeitsmatrizen und –lastvektoren) herleiten (vergleiche dazu die Ausfüh-rungen im Abschnitt 1.4).
Üblicherweise verwendet man in jedem finiten Element eines Elementnetzes den glei-chen Typ von Näherungsfunktionen. Das muß aber nicht immer der Fall sein, wie wirspäter (z.B. bei p-Elementen oder der Verknüpfung unterschiedlicher Elementypen ineinem Netz mit Hilfe von Übergangselemente, der Methode der LAGRANGEschenMultiplikatoren oder der Penalty-Methode) noch sehen werden.
Wir wollen uns nachfolgend typischen Näherungsfunktionen, die wir auch als Ansatz-funktionen (oder im englischen Sprachgebrauch shape functions) bezeichnen, genaueransehen.
Verschiebungsansatzfunktionen (shape functions)Für ein beliebiges finites Element führen wir folgenden Näherungsansatz für die ge-suchte Feldgröße u ein
axu )(P= (2.3-9)
Die Matrix P(x) enthält dabei die Näherungsfunktionen, die im dreidimensionalen Fallvon den Koordinaten xT =[x1, x2, x3] abhängig sind. Die Matrix P(x) enthält meist dieTerme von möglichst vollständigen Polynomen. Wir werden später noch sehen, daßdurch das höchste vollständige Polynom in dem Näherungsansatz die Konvergenzge-schwindigkeit der Lösung bestimmt wird. Natürlich können prinzipiell auch trigono-metrische Funktionen oder andere Funktionensysteme zur Konstruktion geeigneterAnsatzfunktionen verwendet werden. Der Vektor a enthält die Unbekannten des An-satzes, die keine physikalische Bedeutung besitzen. Wir werden nachfolgend zeigen,daß es vorteilhaft ist, diese Unbekannten äquivalent durch andere Unbekannte zu er-setzen, die eine physikalische Bedeutung haben – z.B. als Verschiebungen an denKnotenpunkten eines Elementes.
54 2. Grundlagen der FEM für die Elastomechanik
Die Beispiele für Ansatzfunktionen des Typs (2.3-9) zeigen den typischen Aufbau derMatrix P, die üblicherweise aus Polynomansätzen (daher die Bezeichnung P) gebildetwerden.22
Nachfolgend werden wir die Unbekannten ai äquivalent durch die Knotenverschiebun-gen ersetzt23 Dazu müssen wir den Ansatz (2.3-9) für jeden Knoten eines finiten Ele-mentes anschreiben, wobei sich dann jeweils die an dem Knoten vorhandenen Knoten-verschiebungen ergeben. Fassen wir diese Knotenverschiebungen zu einem Vektor
[ ]nn2111T wvuwvu K=v zusammen, der sämtliche Knotenunbekannten
enthält, ergibt sich
Aav = (2.3-10)
Die Matrix A enthält die Koordinaten der Knotenpunkte, denn diese wurden ja in Peingesetzt, um die quadratische Matrix A zu erhalten. Wir können jetzt aus (2.3-10) dieprimären Unbekannten a durch die Unbekannten v ausdrücken und erhalten
vAa 1−= (2.3-11)
Wir setzen (2.3-11) in den ursprünglichen Ansatz (2.3-9) ein und erhalten
vAxPxu 1−= )()( (2.3-12)
Führen wir noch die Abkürzung1)()( −= AxPxN (2.3-13)
ein, wobei die Matrix N(x) als die Matrix der Formfunktionen bezeichnet werden soll,erhalten wir schließlich die neue Form des Verschiebungsansatzes
vxNxu )()( = , (2.3-14)
in die als Unbekannte die Knotenverschiebungen eingehen, die im Vektor v zusam-mengefaßt sind. In der Fußnote ist beispielhaft für ein Dreieckselement mit drei Kno-
22 Im eindimensionalen Fall (1D) sieht die Funktion bei einem linearen Ansatz im Element folgendermaßen aus:
[ ] aP )x(aa
x1xaau2
121 =
=+=
P(x) ist also eine Matrix, die die gewählten Polynomansätze enthält. Für einen quadratischen Ansatz ergibt sichentsprechend
[ ] aP )x(aaa
xx1xaxaau
3
222
321
1
=
=++= ,
und für ein Dreieckselement mit drei Knoten erhält man
aPu )y,x(a
a
yx1000000yx1
vu
6
1
=
=
= M
Deutlich erkennt man, daß jetzt in jeder der beiden Koordinatenrichtungen x1 und x2 ein linearer Näherungsan-satz eingeführt wurde.23 Das ist jedoch nicht zwangsläufig erforderlich, wie wir später noch sehen werden.
2.3. Die Verschiebungsgrößenmethode 55
tenpunkten (ein finites Element zur Scheibenberechnung, siehe auch Abschnitt 3) dar-gestellt, wie man von (2.3-1) zu (2.3-14) gelangt.24
Anforderungen an die FormfunktionenUm die Konvergenz der Finite-Element-Methode zu sichern, muß bei einer unendlichfeinen Zerlegung des Lösungsgebietes in finite Elemente die Näherungslösung gegendie exakte Lösung des Modellproblems (mathematisches Modell) konvergieren. Umdas zu garantieren, müssen die Formfunktionen eine Reihe von Bedingungen erfüllen.Bei der Behandlung des Ritzschen Verfahrens25 hatten wir schon über Bedingungen fürdie Näherungsfunktionen bei diesem Verfahren diskutiert. Wenn man die FEM alsRitz-Verfahren mit stückweise (elementweise) definierten Näherungsfunktionen auf-faßt, kann man diese Bedingungen auf die FEM übertragen. Wir werden im Abschnitt6 noch einmal auf dieses Problem eingehen und halten hier zunächst nur einige wichti-ge Bedingungen fest.
1. Die Fomfunktionen müssen zwischen benachbarten Elementen stetig sein.Die Funktionen in N(x) müssen so gewählt werden, daß die Kontinuität der Verschie-bungen und ihrer Ableitungen über die Elementgrenzen hinweg garantiert ist. DerGrad dieser Stetigkeitsforderung, d.h. bis zu welcher Ableitung die Bedingung erfülltsein muß, hängt von der höchsten Ableitung ab, die in der Variationsformulierungauftritt. Wir sprechen von einer 0C - Stetigkeit, wenn nur die Verschiebungen selbstzwischen benachbarten Elementränder stetig sind und diese auf dem Gebietsrand dieVerschiebungsrandbedingungen erfüllen. Diese Bedingungen müssen erfüllt sein,wenn im Variationsfunktional maximal erste Ableitungen der Verschiebungen auftre-ten. Von 1C -Stetigkeit sprechen wir, wenn sowohl die Verschiebungen als auch ihre 24 An dem nachfolgend skizzierten Dreieckselement mit 3 Knoten, soll die Herleitung von G aus N erläutert wer-den.
x3
y, v
12
3
u2
v2
u3
v3
u1
v1
x1 x2
x, u
y2
y3
y1
Der Auswertung des Polynomansatzes für das lineare Dreiecksele-ment (siehe Fußnote 9) lautet
aAv ⋅===
5
5
4
3
2
1
33
22
33
11
22
11
3
3
2
1
1
aaaaaa
yx1000
yx1
000yx1000
000yx1000
xx1000yx1
vuvuvu
2
Die Inversion von A und Multiplikation mit P liefert N=PA-1 :
[ ]321321
321
N0N0N00N0N0N1 NNNPAN =
=−=
mit
( ) ( ) ( )[ ] ( ) ( ) ( )[ ] ( ) ( ) ( )[ ]A21
2xxxyxyN,A21
yxxyxyN,A21
yxxyxyN 2112123133123223231 31 ++=++=++=
und ( ) kllklk yxyxxy −= , lklk xx)x( −= , lklk yy)y( −= , ( ) ( ) ( )123123 xyxyxyA2 ++= ,
wobei A der Flächeninhalt des Dreiecks ist.25 Siehe Teil I des Vorlesungsskript Mathematische und Numerische Methoden der Mechanik.
56 2. Grundlagen der FEM für die Elastomechanik
sämtlichen ersten Ableitungen zwischen benachbarten Elementränder stetig sind. DieseBedingung muß erfüllt sein, wenn im Variationsfunktional maximal Ableitungen derVerschiebungen bis zur zweiten Ordnung auftreten. Das ist beispielsweise bei Balken-,Platten- und Schalenmodellen der Fall. Wir sprechen allgemein von einer nC -Stetigkeit, wenn die Verschiebungen und sämtliche Ableitungen bis zur n-ten Ordnungzwischen benachbarten finiten Elementen stetig sind. Diese Bedingung muß dann er-füllt sein, wenn im Variationsfunktional maximal Verschiebungsableitungen bis zurOrdnung n+1 auftreten.
2. Starrkörperverschiebungen dürfen keine Verzerrungen hervorrufenBei Starrkörperverschiebungen und/oder Verdrehungen eines finiten Elementes dürfenim Element keine durch den Näherungsansatz bedingten Verzerrungen und damitSpannungen auftreten. Diese Bedingung klingt zwar selbstverständlich, ist aber in be-stimmten Fällen (z.B. bei Schalenelementen) nicht immer leicht zu erfüllen. Sie istoftmals nur für kleine Verschiebungen näherungsweise erfüllt, was zur Konsequenzhat, daß derartige finite Element nicht für die Berechnung großer Verformungen be-nutzt werden dürfen. Wir wollen an einem einfachen Fall zeigen, welche Konsequenzsich beispielsweise aus dieser Forderung für die Ansatzfunktionen ergibt. Der Ver-schiebungsansatz für ein finites Element lautet entsprechend Gleichung (2.3-14)
∑=
==n
1LLL )()( uxNvxNu , (2.3-15)
wobei n die Anzahl der Knotenpunkte des Elementes ist. Wenn das Element währendder Rechnung sowohl eine Deformation erleidet, die durch den Knotenverschiebungs-vektor [ ]T
nTL
T1
T uuuv LL= beschrieben wird, als auch einer vorgegebenen
Starrkörperverschiebung, z.B. [ ]wvuT ∆∆∆=∆u unterliegt, so ergibt sich für dieVerzerrungen e an einem beliebigen Punkt im Element
uNDuND
uuNDDue
∆⋅
+
=
∆+==
∑∑
∑
==
=n
1LL
n
1LLL
L
n
1LL
)(
)( (2.3-16)
Der erste Summand liefert die sich aus der Deformation ergebenden Verzerrungen.Der zweite Summand rührt aus der Starrkörperverschiebung her und muß Null erge-
ben. Das ist offensichtlich dann der Fall, wenn ∑=
n
1LLN konstant ist, da dann die Anwen-
dung der Differentiationsmatrix D auf diesen Term Null ergibt. Die Summe muß ge-nau gleich 1 sein, da sich dann die Starrkörperverschiebung auch richtig aus dem Ge-samtverschiebungszustand berechnen läßt.
uu
uNuNuuNu
∆+=
∆⋅
+
=
∆+= ∑∑∑===
n
1LL
n
1LLLL
n
1LLgesamt )()( (2.3-17)
2.3. Die Verschiebungsgrößenmethode 57
Daraus folgt als Bedingung für die Formfunktionen:
IN =∑=
n
1LL , (2.3-18)
wobei I die Einheitsmatrix ist. Die Summe aller Formfunktionen muß also genau denWert 1 ergeben, was man z.B. bei dem Ergebnis in Fußnote 24 nachprüfen kann.
3. Die finiten Elemente müssen mindestens einen konstanten Spannungszustand exakt approximieren.Wenn die Knotenverschiebungen einem konstanten Verzerrungszustand entsprechen,muß die Näherungsfunktion für die Verschiebungen diesen auch exakt liefern. Damitschließt dieses Kriterium das vorhergehende mit ein, da ja Starrkörperbewegungen nureinen Sonderfall konstanter Verzerrungen mit dem Wert Null sind.
Auf weitere Probleme bezüglich der Wahl der Näherungsfunktionen und ihres Einflus-ses auf die Genauigkeit und Konvergenzgeschwindigkeit der Lösung wird in Abschnitt7 eingegangen.
58 2. Grundlagen der FEM für die Elastomechanik
2.4. Steifigkeits-, Massen- und Belastungsmatrizen
Alle nachfolgenden Betrachtungen gelten für ein beliebiges finites Einzelelement, d.h.wenn wir das elastische Gesamtpotential π betrachten, berechnen wir nachfolgend nurden Anteil π(e) aus Gleichung (2.3-7). Wir werden im Abschnitt 2.5. noch sehen, wiesich aus diesen Elementbeiträgen das Gesamtsystem aufbauen läßt.
Elementsteifigkeitsmatrix und -kraftvektorAls Ausgangsbasis benutzen wir das Prinzip vom Minimum des Elastischen Potentialsentsprechend Gleichung (2.3-6). Setzen wir in den Zusammenhang zwischen Verzer-rungen und Verschiebungen Gleichung (2.3-2) 0eDue −= den Verschiebungsansatz(2.3-14) ein, ergibt sich
0)( evNDe −= (2.4-1)
mit der Abkürzung
)(NDB = (2.4-2)
erhalten wir daraus
0eBve −= (2.4-3)
Für die Spannungen erhält man dann aus (2.3-4)
0eEEBvEes −== (2.4-4)
Setzen wir die Gleichungen (2.4-1) und (2.4-4) in die Gleichung für das ElastischePotential (2.3-6) ein, erhalten wir
k
n
1k
Tk
O
T
V
T
V0
T0
F
q
])([dO)(
dV)(dV)()(21
FvxNqNv
pNveBvEeBv
∑∫
∫∫
=−−
−−−=π
Durch Ausmultiplizieren ergibt sich daraus
( ) .MindOdVdV
21
dV21
dV21
dV21
F
q
n
1kkk
TT
O
TT
V
TT
V0
T0
V
T0
V0
TT
V
TT
→−−−−
−⋅−=π
∑∫∫∫
∫∫∫
=FxNvqNvpNveEe
EBveeEBvEBvBv
Die Knotenverschiebungsvektoren können ausgeklammert werden, und es ergibt sich
∫ εε−−=πV
oTo
TT dV21
21
EfvKvv
mit den folgenden Abkürzungen.
2.4. Steifigkeits-, Massen- und Belastungsmatrizen 59
Elementsteifigkeitsmatrix
dVV
TEBBK ∫= (2.4-5)
Elementkraftvektor
Fqp0fffff +++= ε (2.4-6)
Kraftvektor infolge von Anfangsverzerrungen
∫=εV
0T dV
0eEBf (2.4-7)
Kraftvektor infolge von Volumenlasten
∫=V
Tp dVpNf (2.4-8)
Kraftvektor infolge von Oberflächenlasten
dOqO
Tq ∫= qNf (2.4-9)
Kraftvektor infolge von Einzelkräften
∑=
=Fn
1kkk
TF )( FxNf (2.4-10)
Die Minimierung des Elastischen Potentials bezüglich der Knotenverschiebungen vliefert für ein Einzelelement26
0fvKv
=−⋅=∂
π∂=δπ (2.4-11)
und wir erhalten die sogenannte Elementsteifigkeitsbeziehung
fKv = , (2.4-12)
die den Zusammen zwischen den an den Knoten angreifenden Kräften f und den Kno-tenverschiebungen v vermittelt.
ElementmassenmatrixFür die Lösung dynamischer Probleme mit Hilfe der FEM benutzen wir als Ausgangs-basis statt des Prinzips vom Minimum des Elastischen Potentials das HamiltonschePrinzip
stationärLdtt
→∫ . (2.4-13)
26 Zur Erinnerung vermerken wir nochmals, daß wir hier nur ein Element betrachten und daher genaugenommenimmer einen Index (e), z.B. p(e) – d.h. bezogen auf das Element (e) – vermerken müßten. Da hier aber alle Be-trachtungen stets auf ein Element bezogen sind, wird der Index (e) weggelassen.
60 2. Grundlagen der FEM für die Elastomechanik
Dabei bezeichnet L die sogenannte Lagrangesche Funktion, die sich folgendermaßenberechnen läßt
π−= TL . (2.4-14)
p ist das schon bekannte Elastische Potential, und T ist die kinetische Energie
∫ρ=V
T dV21
T uu && , (2.4-15)
mit der Dichte ρ, und dtdu
u =& bezeichnet die Geschwindigkeiten. Wenn man bei dy-
namischen Problemen die gleichen Verschiebungsfunktionen wie bei statischen Pro-blemen benutzt und beachtet, daß die Unbekannten des Näherungsansatzes – die Kno-tenverschiebungen – jetzt zeitlich veränderlich sind, erhält man
)t()( vxNu = (2.4-16)
)t()( vxNu && ⋅= (2.4-17)
Einsetzen von (2.4-17) in den Ausdruck für die kinetische Energie (2.4-15) liefert
vxNxNv T &&
ρ= ∫
V
T dV)()(21
T . (2.4-18)
Führt man als Abkürzung die
Massenmatrix
∫ρ=V
T dV)()( xNxNM (2.4-19)
ein, erhält man für die Lagrangesche Funktion
fvKvvvMv TTT
21
21
TL −−=π−= && (2.4-20)
Aus den Lagrangeschen Bewegungsgleichungen (Eulersche Gleichung des Hamilton-schen Prinzips)
0vv
=∂∂
−
∂∂
L
L
dtd
& (2.4-21)
ergibt sich schließlich
0fKvvM =−+&& . (2.4-22)
Wird noch eine geschwindigkeitsproportionale Dämpfung eingeführt, erhält man
fKvvCvM =++ &&& , (2.4-23)
2.4. Steifigkeits-, Massen- und Belastungsmatrizen 61
wobei C die Dämpfungsmatrix bezeichnet, die häufig mittels der sogenannten Be-quemlichkeitshypothese27 aus
KMC β+⋅α= (2.4-24)
bestimmt wird28. Das Differentialgleichungssystem (2.4-23) kann direkt numerisch in-tegriert werden, z.B. mittels der Zentralen Differenzenmethode, dem Newmark-Verfahren oder einem anderen Verfahren29. Wenn man zuerst das Eigenwertproblemlöst, kann das System auf Modalkoordinaten transformiert werden. Das Eigenwertpro-blem erhält man, wenn man für die Lösung des homogenen ungedämpften Differenti-algleichungssystems einen Lösungsansatz der Form
t)cos(ˆ)t( ω= vv (2.4-25)
einführt. Mit
)t cos(ˆ)t( 2 ωω−= vv&& (2.4-26)
erhält man mit f=0 aus (2.4-22) das Eigenwertproblem
( ) 0vMK =ω− ˆ2 (2.4-27)
dessen Lösung die Eigenwerte und Eigenvektoren des Systems liefert.
27 Auch als Rayleigh-Dämpfung bezeichnet.28 Durch diese Annahme erreicht man, daß das Differentialgleichungssystem mitttels der Eigenvektoren modalentkoppelt werden kann.29 Siehe dazu die Vorlesungen Mathematische und Numerische Methoden der Mechanik - Teil II.
62 2. Grundlagen der FEM für die Elastomechanik
2.5. Das Gesamtsystem
Assemblierung der Elementbeiträge zum GesamtsystemDie lokale (d.h. elementweise) Definition der Ansatzfunktionen zur Approximation derunbekannten Feldgrößen, die jeweils nur in einem Element ungleich Null sind, führtdazu, daß das Elastische Gesamtpotential π als Summe der Elastischen Potentiale dereinzelnen finiten Elemente π(e) aufgeschrieben werden kann (siehe Gleichung 2.3-6).
∑π=πe
)e( (2.5-1)
Der hochgestellte Index (e), der auf diesen Sachverhalt verweist, wurde im Abschnitt2.4. zur Vereinfachung weggelassen,30 wird aber hier wieder konsequent eingeführt. .
Wir hatten bisher stillschweigend vorausgesetzt, daß die Knoten eines finiten Elemen-tes von 1 bis n(e) durchnumeriert sind, wenn mit n(e) die Anzahl der Elementknoten be-zeichnet wird. In einer Vernetzung ist jedem Elementknoten eine globale Nummer zu-geordnet, die sich aus der Durchnumerierung aller Knoten eines Netzes ergibt. Im Bild2.5-1 ist das am Beispiel eines Dreieckselementes mit 3 Knoten dargestellt. Den glo-balen Knotennummern L, M und N entsprechen die lokalen Knotennummern 1, 2 und3. Mit dieser Zuordnung der lokalen Numerierung eines Elementes (diese muß natür-lich für jedes finiten Element eindeutig definiert sein) zur globalen Numerierung einerVernetzung, ist die Topologie eines Netzes beschrieben. Wenn man zusätzlich noch dieKnotenkoordinaten aller Knoten kennt, ist das Netz auch geometrisch eindeutig be-schrieben und kann z.B. graphisch dargestellt werden. Jedem Knoten sind ndof Unbe-kannten zugeordnet.31
Die Abbildung zwischen den Unbekannten v im globalen System (dem Gesamtnetz)und dem Einzelelement v(e) kann formal durch eine aus Nullen und Einsen bestehendeAbbildungsmatrix S(e) in der Form 30 Damit gelten die Formeln in Abschnitt 2.3 auch für den Fall, daß im Sinne des Ritzschen Verfahrens eine kon-tinuierliche Ansatzfunktion für das Gesamtgebiet eingeführt wird.31 So hat z.B. im 2D-Fall (Scheibenproblem) jeder Konten L mindestens zwei Unbekannte - die Knotenverschie-bungen in x- und in y-Richtung (hier als uL und vL bezeichnet).
eL
K
M
e1=L
3=K
2=M
11 2 3
45 6 7
2 3
Bild 2.5-1 Globale Knotennummern L, M, K eines Dreieckselementes (e) und zugehörige lokale Kontennummern 1, 2, 3
2.5. Das Gesamtsystem 63
vSv )e()e( = (2.5-2)
beschrieben werden, wobei sich die Größen ohne hochgestellten Index (e) stets auf dasGesamtsystem beziehen. Mit (2.5-2) werden aus dem Gesamtvektor der Unbekannten,die zum e-ten Element gehörenden Unbekannten herausgezogen. Setzt man in Glei-chung (2.5-1) die Elastischen Potentiale der Einzelelemente ein, erhält man
∑ ∑
−=π=π
e e
)e(T)e()e()e(T)e()e(
21
fvvKv
Mit Gleichung (2.5-2) ergibt sich daraus
∑∑ −
=π
e
)e(T)e(T
e
)e()e(T)e(T
21
fSvvSKSv (2.5-3)
Wenn man jetzt die Minimierung von π bezüglich aller Unbekannten, die im Vektor vzusammengefaßt sind, bildet, ergibt sich
0fSvSKSv
=−
=
∂π∂ ∑∑
e
)e()e(
e
)e()e(T)e( (2.5-4)
oder0fKv =− . (2.5-5)
Der erste Summand ist die gesuchte Gesamtsteifigkeitsmatrix K des Systems
∑=e
)e()e(T)e( SKSK , (2.5-6)
und der zweite Summand ist der Systemkraftvektor
∑=e
)e()e( fSf (2.5-7)
Man erkennt, daß sich die Systemmatrix und der Systemkraftvektor additiv aus denElementbeiträgen zusammensetzen. Im Abschnitt 1.4 wurde am Beispiel eines aus dreieindimensionalen finiten Elementen bestehenden Systems gezeigt, wie man über dieMinimierung des Gesamtpotentials zum Gesamtgleichungssystem (2.5-5) gelangt. Dasgleiche Ergebnis läßt sich auch durch die Erfüllung des Kraftgleichgewichts an jedemKnotenpunkt erzielen. Die obige Darstellung ist aber allgemeingültiger. Über die Ma-trix S(e) werden die Komponenten der e-ten Elementsteifigkeitsmatrix und des e-tenElementkraftverkors auf die Positionen in dem Gesamtgleichungssystem aufaddiert,die sich aus der Abbildung der lokalen Numerierung der Knoten (lokale Knotenfrei-heitsgradnummern) auf die globalen Knotennummern (globale Knotenfreiheitsgrad-nummern) ergeben. Für diese Einspeicherung ist also die Tologiematrix der Vernet-zung erforderlich, mit der die Position in der Gesamtsteifigkeitsmatrix berechnet wer-den kann, auf die ein Element aus der Elementsteifigkeitsmatrix aufaddiert werdenmuß. Programmtechnisch läßt sich das über einen Einspeicherungsalgorithmus für dasElementgleichungssystem realisieren. Die Beschreibung über die Gleichungen (2.5-6)und (2.5-7) hat also nur eine symbolische Bedeutung.
64 2. Grundlagen der FEM für die Elastomechanik
Beispiel 2.5-1An dem im Bild 2.5-1 dargestellten Beispiel soll der Einspeicherungsalgorithmus verdeutlicht werden.Wir nehmen an, daß an jedem Knoten 2 Unbekannte, dieVerschiebung in x-Richtung und die Ver-schiebung in y Richtung, vorliegen. Die Elementsteifigkeitsbeziehung für das Dreieckselement (e) mitdrei Knoten hat dann sechs Gleichungen. Die knotenbezogene Struktur lautet
=
)e(3
)e(2
)e(1
)e(3
)e(2
)e(1
)e(33
)e(23
)e(22
)e(13
)e(12
)e(11
fff
vvv
KKKKKK
,
wobei die Indizes 1,2 und 3 die lokalen Knotennummern des Elementes bezeichnen. Da jedem Knotenzwei Unbekannte zugeordnet sind, hat das Elementgleichungssystem die Form
=
)e(3y
)e(3x
)e(2y
)e(2x
)e(1y
)e(1x
)e(3
)e(3
)e(2
)e(2
)e(1
)e(1
)e(33
)e(33
)e(33
)e(23
)e(23
)e(23
)e(23
)e(22
)e(22
)e(22
)e(13
)e(13
)e(13
)e(13
)e(12
)e(12
)e(1211
)e(12
)e(11
)e(11
)e(11
ff
ff
ff
vu
vu
vu
k
kkhsymmetrisc
kk
kk
k
kk
kk
kk
kk
kk
k
kk
22
1211
2221
1211
22
1211
2221
1211
2221
1211
Diese Art der Darstellung ermöglicht es nun, den Einspeicheralgorithmus in einer allgemeingültigenForm anzugeben. Dazu zeigen wir nachfolgend zunächst einen Ausschnitt aus der TopologiematrixMtopo (vergleiche Bild 2.5-1):
=
MMMMMM
MMMKML
732672621
topoM
Wenn wir nun voraussetzten, daß die Gesamtsteifigkeitsmatrix als quadratische Matrix gespeichertwerden soll und alle Knoten die gleiche Anzahl von Unbekannten ndof aufweisen (in unserem Beispielist ndof=2), so ergibt sich für die Einspeicherung des Matrixelementes )e(
IJ ijk , wobei I,J=1,2,3 und
i,j=1,2 annehmen können, folgende Position (k,m) in der Gesamtmatrix K
)e(IJkmkm ij
kKK +=
mit( )( ) jn1)J,e(m
in1)I,e(k
doftopo
doftopo
+⋅−=+⋅−=
MM
Da die Matrix K symmetrisch ist, brauchen nur die Elemente km ≥ aufgebaut zu werden. Um Spei-cherplatz zu sparen, erfolgt die Speicherung von K üblicherweise in einem speziellen Format. Ein ein-fache Variante eines solchen speziellen Formats stellt die folgende Rechteckmatrix dar:
2.5. Das Gesamtsystem 65
44444 844444 76
MMMMMMM
LLMMMMMM
LLLL
ibw
nn
ibw,i1i,iii
ibw22322
ibw11211
0000k0000000000
kkk
kkkkkk
= +K ,
bei der die Hauptiagonale in der ersten Spalte steht. Mit ibw wird die sogenannte Bandweite der Ma-trix bezeichnet.32 Nur innerhalb der Bandweite sind in der Matrix Element ungleich Null vorhanden.Man spricht in solchen Fällen auch von Bandmatrizen. Die Bandweite berechnet sich folgendermaßenaus der Differenz der Knotennummern, wenn man voraussetzt, daß alle Knoten die gleiche Anzahlvon Unbekannten ndof aufweisen.
)1dmax(nibw )e(
)e(dof +⋅= (2.5-8)
d(e) bezeichnet die maximale Differenz der Knotennummern des Elementes (e).
Die Gesamtsteifigkeitsmatrix K ist positiv-semidefinit. Sie wird nach Einarbeitung derRandbedingungen, die mindestens die Starrkörperbewegungen des Systems verhindernmüssen, positiv definit.
RandbedingungenDie Verschiebungsrandbedingungen33 lassen sich auf unterschiedliche Art und Weisein des Gesamtgleichungssystem einarbeiten. Wenn einzelne Verschiebungsfreiheits-grade Null sein müssen, kann das durch Streichen der betreffenden Zeilen und Spaltengeschehen. Als Konsequenz ist ein Umordnen der Matrix erforderlich, was bei derGröße der Gleichungssysteme mit einem erheblichen Rechenaufwand verbunden ist.Alternativ können die Randbedingungen auch mit Hilfe der Penalty-Methode oder mitder Methode der Langrangeschen Multiplikatoren berücksichtigt werden. Dabei er-möglicht es die Penalty-Methode, auch komplizierte Rand-, Zwangs- und allgemeineKoppelbedingungen relativ einfach im Gleichungssystem zu berücksichtigen,34 ohnedaß wie bei der Methode der Lagrangeschen Multiplikatoren das Gleichungssystemvergrößert wird. Diesem Vorteil steht der Nachteil gegenüber, daß die Randbedingun-gen nur näherungsweise und zwar in Abhängigkeit von der Größe der Penalty-Zahl er-füllt werden. Mit steigenden Penalty-Zahl werden die Randbedingungen immer bessererfüllt, in gleichem Maße steigt aber auch die Konditionszahl der Koeffizientenmatrix, 32 Manchmal spricht man auch von halber Bandweite, da ja genaugenommen die symmetrische linke Seite derMatrix auch noch vorhanden ist. Bei unsymmetrischen Matrizen kann man eine linke und rechte Bandweite an-geben, die nicht gleich sein müssen.33 Auch als wesentliche Randbedingungen oder als kinematische Randbedingungen bezeichnet.34 Das ist der Hauptgrund, weshalb in COSAR alle Rand- und Zwangsbedingungen mit Hilfe der Penalty-Methode berücksichtigt werden.
66 2. Grundlagen der FEM für die Elastomechanik
die bei zu großer Penalty-Zahl singulär werden kann.35 Wir werden später nochmalsauf diesen Punkt zurückkommen und verweisen hier auf die Ausführungen im Ab-schnitt 1.4.
Lösung des GleichungssystemsNach der Einarbeitung der Randbedingung kann das Gleichungssystem (2.5-5) gelöstwerden, wobei dafür eine große Anzahl von Lösungsverfahren zur Verfügung steht.36
Wegen der großen Anzahl von Gleichungen37 und der dadurch bedingten hohen Re-chenzeiten für die Lösung des Gleichungssystems kommt es bei der Anwendung di-rekter Lösungsverfahren vor allem darauf an, eine möglichst geringe Bandweite zu er-reichen. Das kann durch eine interne Umnumerierung der Knoten erreicht werden. Da-für steht eine Vielzahl von Verfahren zur Verfügung, die aus einer schlechten Aus-gangsnumerierung, die meist das Ergebnis automatischer Netzgenerierungssoftwareist, eine „brauchbare“ interne Neunumerierung erzeugen. Wenn man beachtet, daß sichdie Rechenzeit bei einem direkten Lösungsverfahren aus
2CPU )ibw(nct ⋅⋅= (2.5-9)
berechnen läßt, wobei c eine vom Computer abhängige Konstante ist und die Koeffizi-entenmatrix das Gleichungssystem das Format (n,ibw) hat (vergleiche die im Beispiel2.5-1 angegebene Rechteckform für K), wird deutlich, welche Bedeutung eine Reduk-tion der Bandweite für die Gleichungslösung hat.38 Neben Verfahren, die auf zykli-schen Tauschalgorithmen zur Reduktion der Bandweite basieren, haben sich vor allemgraphentheoretische Verfahren als besonders effizient erwiesen. Neben den direkten Lösungsverfahren kommen zunehmend auch iterative Ver-fahren zum Einsatz, wobei sich hier vor allem vorkonditionierte Gradientenverfahren(pcg-Verfahren) als besonders leistungsfähig erwiesen haben.
SpannungsberechnungNachdem das Gleichungssystem gelöst wurde, sind die Verschiebungen für jedenKnotenpunkt bekannt. Diese können in den Elementverschiebungsansatz eingesetztwerden, so daß die Verschiebungen auch innerhalb eines Elementes berechnet werdenkönnen. Damit lassen sich nun eine Vielzahl von sekundären Ergebnisgrößen berech-nen. Die Spannungsberechnung kann beispielsweise mit Hilfe der Gleichung (2.4-4)erfolgen und man erhält
)e(0
)e()e()e()e()e( eEvBEs −= (2.5-10)
Dabei soll hier der hochgestellte Index (e) nochmals hervorheben, daß die Spannungenelementbezogene Größen sind. Während die Verschiebung über das ganze Lösungsge-biet mindestens stückweise stetige Funktionen sind, treten üblicherweise Sprünge in
35 Empfehlung für die Wahl der Penalty-Zahl: ii
)i(
p kmin10 ⋅=µ , mit p≈4...6.
36 Siehe dazu die Vorlesung Mathematische und Numerische Methoden der Mechanik - Teil II.37 Bei der Lösung praktische Aufgabenstellungen sind mehrere Millionen Unbekannte keine Seltenheit. Siehe da-zu beispielsweise die Angaben in den Diplomarbeiten, die im Abschnitt 1.1. zitiert wurden.38 Für das Programmsystem COSAR wurde beispielsweise ein schnelles graphentheoretisches Verfahren entwik-kelt und implementiert, das auf dem Cuthill-McKee Verfahren basiert.
2.5. Das Gesamtsystem 67
den Spannungen zwischen benachbarten Elementen auf. Die Spannungen sind also inder Regel zwischen benachbarten Elementen nicht stetig. Die Größe der Spannungs-sprünge ist ein Maß für die Genauigkeit der Näherungslösung. Wir werden darauf imAbschnitt 7 noch genauer eingehen. Meist erfolgt programmintern eine einfache Mit-telwertbildung der Spannungen an den Knoten, wobei in die Mittelwertbildung alleKnotenspannungen aus den angrenzenden Elementen eingehen. Dadurch wird demAnwender des Programmes ein stetiger Verlauf der Spannungen vorgetäuscht, der ei-gentlich so nicht vorhanden ist. In vielen Fällen führt dieser Mittelungsprozeß aller-dings zu geglätteten Spannungen, die genauer als die Einzelspannungen in den Ele-menten sind. Hier ist aber Vorsicht geboten, weil sich gerade dadurch teilweise auchvöllig unbrauchbare oder unsinnige Ergebnisse ergeben können wie wir später nochsehen werden. Meist wird ein einfaches arithmetisches Mittel berechnet, das sich aus
∑=
σ=σm
1
)k()(ij
)k(ij m
1
ll (2.5-11)
berechnen läßt. Dabei bezeichnet (k) den Knoten, an dem die Mittelung erfolgt,(k)
)(ij lσ sind die Spannungen aus dem l-ten Element, das an den Knoten (k) angrenzt, undm ist die Gesamtanzahl der Elemente, die an den Knoten (k) angrenzen.
68 3. Simplex-Elemente
3. Simplex-Elemente
3.1. Das 2-Knoten Stabelement
Das 2-Knoten Stabelement wurde bereits im Abschnitt 1.4 als Einführungsbeispielausführlich behandelt. Hier dient das Element (vor allem aus didaktischen Gründen)dazu, die im Abschnitt 2 dargestellten allgemeinen Grundlagen an einem ersten, sehreinfachen Beispiel anzuwenden. Als neuen Aspekt wollen wir hier erstmals die Nut-zung von dimensionslosen - sogenannten natürlichen Koordinaten - einführen, die sichfür die Ableitung von finiten Elementen im weiteren als nützlich erweisen werden.
Im Bild 3.1-1 ist das Element abgebildet.1 In dem natürlichen Koordinatensystem ξ hatdas Element die Länge 2. Der Zusammenhang zwischen der natürlichen Koordinate ξ,deren Ursprung der Elementmittelpunkt liegt, und der x-Koordinate ergibt sich aus
ξ+= l21
xx 0 (3.1-1)
ξ=
ξξ
=dd2
dd
dxd
dxd
l (3.1-2)
In den natürlichen Koordinaten lassen sich die Formfunktionen N1(x) und N2(x) fürdas 2-Knotenlement folgendermaßen schreiben2
)1(21
N1 ξ−= (3.1-3a)
)1(21
N2 ξ+= (3.1-3b)
Die Matrix der Formfunktionen N hat die Form
ξ+ξ−= 1(
21
)1(21
N (3.1-4)
Der Verschiebungsansatz kann daher folgendermaßen geschrieben werden:
[ ]
ξ+ξ−=
==
2
1
2
121 u
u1(
21
)1(21
uu
NNu Nv
1 Der Index bezeichnet hier stets die Knotennummer.2 Vergleiche Abschnitt 1.4 Gleichung (1.4-11). Die Summe der Formfunktionen ergibt 1, wie es Gleichung (2.3-18) fordert, um Starrkörperverschiebungen korrekt zu erfassen.
xl
21
x=x1 x=x2
= -1 = 1
Bild 3.1-1: 2-Knoten Stabelement mit natürlichen Koordinaten
3.1. Das 2-Knoten Stabelement 69
Für die Berechnung der Elementsteifigkeitsmatrix brauchen wir noch die Matrix B.Für ein Stabelement ist der Zusammenhang zwischen Dehnungen und Verschiebungengegeben durch
BvDNvD ====ε udxdu
(3.1-5)
mit
[ ]111
dd2
dxd
−=
ξ
=
=
llB (3.1-6)
Die Hookesche Matrix besteht hier aus dem Elastizitätsmodul multipliziert mit derStabfläche A, die sich durch die Reduktion des Volumenintegrals auf ein Linieninte-gral ergibt.
[ ]EA=E (3.1-7)
Da wir die natürlichen Koordinaten auch für die Berechnung der Integrale benutzenwollen, müssen wir noch folgende Transformation beachten
ξ= ∫∫−
d2
dx1
1
x
x
2
1
lKK (3.1-8)
Jetzt können die Elementmatrizen berechnet werden, und wir erhalten:3
Elementsteifigkeitsmatrix
[ ]∫∫−−
ξ−⋅
−=ξ=
1
1
1
1
T d1111
21EA
d2
EAl
lBBK
−
−=
1111EA
lK (3.1-9)
Elementkraftvektoren
Volumenlasten infolge Eigengewicht in x-Richtung4
∫∫−−
ξ
ξ+ξ−
ρ=ξ=1
1
1
1
Tp d
211
21
gAd2
All
pNf
ρ=
11
gA21
p lf (3-1-10)
Massenmatrix
ξρ= ∫−
d2
A T1
1
lNNM
3 Dabei nehmen wir an, daß der Elastizätsmodul und die Elementfläche elementweise konstant sind.4 Hier ist gp ρ= .
70 3. Simplex-Elemente
[ ] ξξ+ξ−
ξ+ξ−
ρ= ∫ d)1()1(21
)1()1(
21
A21
lM
ρ=
2112
A61
lM (3.1-11)
Um Rechenzeit zu sparen und zu schnelleren numerischen Algorithmen zu gelangen,wird häufig mit diagonalisierten Massenmatrizen gearbeitet.5 Die einfachste Möglich-keit einer Diagonalisierung besteht darin, eine Zeilensummenbildung vorzunehmen.6
Diese führt hier zu
ρ=
1001
A21
lM (3.1-12)
Im Abschnitt 1.4 wurde bereits gezeigt wie man die auf ein lokales Stabsystem bezo-genen Elementmatrizen und -vektoren auf ein globalen Bezugssystem transformierenkann. Das soll daher hier nur noch einmal kurz rekapituliert werden.7
Stabelement im Raum
Wir bezeichnen jetzt zur Unterscheidung die auf das lokale Stabkoordinatensystem be-zogenen Größen mit einem Querstrich. Für die Transformation der Verschiebungenzwischen dem lokalen und dem globalen Koordinatensystem gilt mit )x,xcos(c xx =
Tvv =
=
=
3
3
2
2
1
1
zxyxxx
zxyxxx
2
1
vuvuvu
ccc000000ccc
vu
(3.1-13)
Analog gilt für den Kraftvektor
Tff = (3.1-14)
Einsetzen in die lokale Steifigkeitsbeziehung fvK = führt zu
fTvK = (3.1-15)
Durch Linksmultiplikation mit TT folgt daraus
fTTvKT TT =
5 Eine solche diagonalisierte Matrix wird im Englischen als lumped matrix bezeichnet.6 Durch die Zeilensummenbildung werden die Gesamtmasse des Elementes und auch seine Trägheitseigen-schaften nicht verändert. Allerdings kann die Zeilensummenbildung zu negativen Hauptdiagonalelementen füh-ren, was bei der Anwendung bestimmter numerischer Verfahren von Nachteil ist. In diesem Fall muß entwederdas numerische Verfahren modifiziert werden oder die Diagonalisierung muß so erfolgen, daß eine positiv defi-nite Diagonalmatrix entsteht.7 Die Ableitung der Elemente in einem lokalen Koordinatensystem und anschließende Transformation auf einallgemeines Bezugssystem ist eine übliche Herangehensweise bei der Ableitung von finiten Elementen.
3.1. Das 2-Knoten Stabelement 71
fKv =mit
TKTK T= (3.1-16)
fTf T= (3.1-17)
Die Steifigkeitsmatrix K ergibt sich aus (3.1-15) durch Ausmultiplizieren.
K
−−−−−−−−−
=
2zx
zxyx2
yx
zxxxyxxx2
xx
2zxzxyxzxxx
2zx
zxyx2
yxyxxxzxyx2
yx
zxxxyxxx2
xxzxxxyxxx2
xx
cccc.symmccccccccccc
cccccccccccccccccc
EAl
(3.1-18)
Nach dem die Verschiebungen im Gesamtsystem berechnet wurden, erfolgt die Rück-transformation mit (3.1-13). Damit können die Stabkräfte und Spannungen wieder imlokalen Koordinatensystem berechnet werden.
72 3. Simplex-Elemente
3.2. Das 2-Knoten Balkenelement (Bernoulli-Balken)Während das 2-Knoten Stabelement das einfachste C0-stetige finite Element ist, dasstellvertretend für 2D- und 3D-Elemente steht, ist das 2-Knoten Balkenelement daseinfachste C1-stetige finite Element, das stellvertretend für Platten und Schalenelementsteht. Wir wollen hier zunächst nur ein ebenes Element mit gerader Stabachse bei klei-nen Verformungen betrachten, so daß die Biegedeformationen des Balken unabhängigvon der Längsdeformation sind.8 Die Längsverformungen haben wir ja gerade im Ab-schnitt 3.1 behandelt. Dieser Anteil kann bei kleinen Verformungen einfach additivder Steifigkeitsbeziehung des Biegeanteils hinzugefügt werden, um ein komplettesebenes Balkenelement zu erhalten.9 Die nachfolgende Ableitung des Balkenelementessetzt voraus, daß die Deformation des Balkens infolge Querkraftschub vernachlässigtwerden kann, was bei schlanken Balkentragwerken stets zulässig ist.10 Derartige Bal-ken werden auch als Bernoulli-Balken bezeichnet.11 Das Deformationsverhalten einesBalkenelementes mit 2-Knoten ist in Bild 3.2-1 abgebildet.
Bild 3.2-1: 2-Knoten Balkenelement in der x-z Ebene
Das Elastische Potential für ein gerades Balkenelement lautet12
( )( ) ( )
∑∑∫∫ ϕ−−−′′=πk
kkk
kk
l
0
l
0
2 MwFdxqwdxwEI21
(3.2-1)
Das entspricht der Gleichung (2.3-1) wenn man für die Dehnung xε=e und für dieSpannung xσ=s einsetzt.13 Es gilt 8 Bei gekrümmter Stabachse sind in den Differentialgleichungen Längskräfte und Momente miteinander ver-knüpft (Längskräfte rufen beispielsweise Biegemomente hervor) und können daher nicht unabhängig voneinan-der betrachtet werden. Man kann natürlich einen gekrümmten Balken auch durch stückweise gerade Balkenele-mente approximieren und auf diesem Weg eine brauchbare Näherungslösung erhalten.9 Analog kann die Steifigkeitsbeziehung für den 2D-Fall auf das komplette 3D-Balkenelement erweitert werden.Dazu sind zusätzlich der Torsionsanteil erforderlich, der analog zum Stabelement aufgebaut ist, sowie der Bie-geanteil in der anderen Balkenebene, der analog zu dem hier behandelten Biegeanteil aufgebaut ist. Es entstehtauf dies Weise ein 3D-Balkenelement mit gerader Stabachse, das 12 Unbekannte (Freiheitsgrade) aufweist.10 Ein Kriterium dafür ist das Verhältnis der Querschnittsabmessung (mit H bezeichnet) zur Gesamtlänge des
Balkens (mit L bezeichnet). Der Querkraftschub kann bei 20HL ≥ meist vernachlässigt werden.
11 Im Unterschied zum Timoshenko-Balken, bei dem der Querkraftschub näherungsweise berücksichtigt wird.12 Siehe Vorlesung Mathematische und Numerische Methoden der Mechanik - Teil I13 Der Balken wird hier nur in der x-z Ebene betrachtet. In die Energieformulierung geht entsprechen der Ber-noulli-Hypothese nur die Biegeformänderungearbeit ein.
x,ul
21
x=x1 x=x2
= -1 = 1
w1
w2
q(x)
z, w
1
2
3.2. Das 2-Knoten Balkenelement (Bernoulli-Balken) 73
zDwdx
wdz
2
2
x −=−=ε (3.2-2)
DwzEE xx −=ε=σ (3.2-3)mit
2
2
22
2
dd4
dxd
Dξ
==l
(3.2-4)
Setzt man (3.2-2) und (3.2-3) in die Formel für das Elastische Potential (2.3-1) ein, er-hält man
,dx)Dw(EI
21
Ddx)Dw(E)dA(21
dV)Dw(Ez21
dV21
2
2
V V
22xx
K
KKK
+=
+=+=+σε=π
∫
∫ ∫∫ ∫
l
l A
2z
wobei ∫=A
2dAzI das Flächenträgheitsmoment ist. Man erkennt, das sich aus der all-
gemeinen Formel (2.3-1) der Spezialfall der Balkenbiegung (3.2-1) ableiten läßt.
Formfunktionen
In dem Elastischen Potential steht die zweite Ableitung der Durchbiegung, so daß daszu entwickelnde finite Element offensichtlich C1 stetig sein muß, d.h an den beidenKnotenpunkten müssen die Durchbiegung w als auch der Biegewinkel ϕ zwischen be-nachbarten Elementen stetig sein.14 Damit ergeben sich insgesamt 4 Unbekannte (w1,ϕ1, w2, ϕ2), so daß für die Approximation der Durchbiegung ein kubischer Polynoman-satz erforderlich ist. Dieser lautet
( ) ( )aP ξ=ξ+ξ+ξ+=ξ 34
2321 aaaaw (3.2-5)
Für den Biegewinkel ergibt sich daraus
( )2432 a3a2a
2ddw2
dxdw
w)( ξ+ξ+=ξ
==′=ξϕll
(3.2-6)
Wir schreiben jetzt den Ansatz (3.2-5) und (3.2-6) an den beiden Knoten 1 und 2 desElementes an und erhalten
Aa=
−
−−
=
ϕ
ϕ
4
3
2
1
2
2
1
1
aaaa
3210111132101111
21w
21w
l
l (3.2-7)
Die Inversion der Matrix A liefert
14 Falls dort nicht ein Gelenk oder eine andere spezielle Verbindung vorliegt.
74 3. Simplex-Elemente
−−
−−−−
=−
25,025,025,025,025,0025,0025,075,025,075,025,05,025,05,0
1A (3.2-8)
Führt man jetzt noch den Vektor v der Unbekannten des Elementes mit
ϕ
ϕ=
2
2
1
1
w
w
v (3.2-9)
ein und bildet entsprechend Gleichung (2.3-13) die Matrix G=P(ξ)A-1, so erhält manfür den Verschiebungsansatz
( ) [ ]
ϕ
ϕ=ξ=ξ=ξ −
2
2
1
1
43211
w
w
NNNN)()(w vNvAP (3.2-10)
mit
( ) ( )ξ+ξ−= 2141
N 21 (3.2-11a)
( )( )ξ−ξ−= 1181
N 22 l (3.2-11b)
( ) ( )ξ−ξ+= 2141
N 23 (3.2-11c)
( )( )ξ+ξ−−= 1181
N 24 l (3.2-11d)
Die Formfunktionen sind im Bild 3.2-2 grafisch dargestellt.
Bild 3.2-2: Grafische Darstellung der Formfunktionen für das 2-Knoten Balkenelement
3.2. Das 2-Knoten Balkenelement (Bernoulli-Balken) 75
Elementsteifigkeitsmatrix
Für die Berechnung der Elementsteifigkeitsmatrix werden die zweiten Ableitungen derFormfunktionen nach der x-Koordinate benötigt.
[ ]43212
2
22
2
BBBB)(dd4
)(dxd
)(D)( =ξξ
=ξ=ξ=ξ NNNBl
(3.2-12)
Daraus ergibt sich
ξ=ξ216
)(Bl
(3.2-13a)
( )ξ+−==ξ 311
)(B2 l (3.2-13b)
ξ−=ξ236
)(Bl
(3.2-13c)
( )ξ+=ξ 311
)(B4 l (3.2-13d)
Die Elementsteifigkeitsmatrix kann dann folgendermaßen berechnet werden
∫−
ξ=1
1
T d2
EIl
BBK (3.2-14)
Für konstantes EI erhält man daraus15
−−−
=
2
22
3
l4.symml612
l2l6l4l612l612
lEI
K (3.2-15)
Kraftvektor
Wir wollen hier nur den Kraftvektor infolge einer linear veränderlichen Linienlast an-geben. Der Elementkraftvektor hat die Form
=
2
2
1
1
MFMF
f (3.2-16)
15 Für einen Balken mit linear veränderlichem Trägheitsmoment (I1 am Knoten 1 und I2 am Knoten 2) ergibt sich
α+α+−α+
α+α+−α+α+α+−α+α+
=
)31()21(2)1(6
)1()2(2)3()21((2)1(6)2(2)1(6
EI
2
22
31
lsymm.l
lllll
lK mit
2
1
II
=α
76 3. Simplex-Elemente
Für eine linear veränderliche Linienlast mit der Intensität 1q am Knoten 1 und 2q amKnoten 2, wird der Verlauf der Funktion )(q ξ längs der Balkenachse durch
2qq
2qq
)(q 1212 ++ξ
−=ξ (3.2-17)
beschrieben. Der äquivalente Knotenkraftvektor berechnet sich nach Gleichung (2.4-9)folgendermaßen:
∫−
ξξξ=1
1
Tq d
2)(q)(
lNf (3.2-18)
( )
( )
( )
( )
+−
+
+
+
=
21
2
21
21
2
21
q
q3q260
q7q320
q2q360
q3q720
l
l
l
l
f (3.2-19)
Für qqq 21 == ergibt sich
−
=
12
21
12
21
l
l
lf (3.2-20)
Man erkennt, daß sich in diesem Fall auch Momente im Kraftvektor ergeben, obwohlnur eine Linienlast wirkt.16
Massenmatrix
Die Massenmatrix berechnet sich nach Gleichung (2.4-19) zu
∫−
ξξξρ=1
1
T d2
)()(Al
GGM (3.2-21)
−−−
ρ=
2
22
4symm221563134135422156
420A
ll
lllll
lM (3.2-22)
16 Die äquivalenten Knotenkräfte und Momente entsprechen den Lagerreaktionen des beidseitig eingespanntenund durch eine Linienlast belasteten Elementes.
3.2. Das 2-Knoten Balkenelement (Bernoulli-Balken) 77
Stabknickung als Beispiel für ein Stabilitätsproblem
Wir wollen nachfolgend noch kurz auf die Berechnung von Knickkräften unter Ver-wendung des oben angegebenen Balkenelementes eingehen. Bei Knickproblemen nachTheorie 2-ter Ordnung müssen wir den Einfluß der Verformung auf das Gleichgewichtberücksichtigen, können aber noch die für kleine Verformungen abgeleitete Differenti-algleichung der Balkenbiegung verwenden. Diese lautet17
( )( )
( )( )
.Min.dxwF21
dxwEI21
l
2N
l
2 →′+′′=π ∫∫ (3.2-23)
Mit FN ist hier die Schnittkraft in Stablängsrichtung bezeichnet. Setzen wir in Glei-chung (3.2-23) den Näherungsansatz (3.2-10) ein, erhalten wir
.Min21
21 TT →+=π σvKvKvv (3.2-24)
mit
( )
ξ=
=
∫
∫
−
σ
d2
F
dxF
1
1
TN
l
TN
lLL
LLK
(3.2-25)
und
[ ]4321 LLLL)(dd2
)(dxd
)( =ξξ
=ξ=ξ NNLl
(3.2-26)
mit
)33(21
)(L 21 ξ+−=ξ
l (3.2-27a)
)321(41
)(L 22 ξ+ξ−−=ξ (3.2-27b)
)33(21
)(L 23 ξ−=ξ
l (3.2-27c)
)321(41
)(L 24 ξ−ξ−−=ξ (3.2-27d)
Die Matrix Kσ wird auch als geometrische Steifigkeitsmatrixbezeichnet, weil sie denEinfluß der verformten Geometrie auf die Biegespannungen berücksichtigt. Für denFall, daß FN elementweise konstant ist, ergibt sich mit FFN = aus (3.2-25)
−−−
−
=σ
2
22
l4.symml336
ll3l4l336l336
30F
lK (3.2-28)
17 Die erste Variation δπ=0 liefert als Eulersche Gleichung 0wF)wEI( N =′′−′′′′ , die wir bereits als Differential-gleichung für Knickstäbe kennengelernt haben (siehe MNMM-Teil I)
78 3. Simplex-Elemente
Die Minimierung von π (3.2-24) bezüglich der eingeführten Unbekannten des Nähe-rungsansatzes führt zu
.0vKKvv
=+=∂π∂
σ (3.2-29)
Für positive Längskräfte liefert (3.2-29) nur die triviale Lösung v=0. Falls FN auch ne-gativ sein kann, schreiben wir
FFN λ−= , (3.2-30)
wobei F eine bekannte Bezugskraft (Einheitskraft) ist. Den Eigenwert λ kann man vorKσ schreiben, und man erhält das Eigenwertproblem
( ) 0vKK =− σ? , (3.2-31)
Für ein Element mit konstanter Längskraft lautet das Eigenwertproblem:18
0=
ϕ
ϕ
−−−
−
λ−
−−−
2
2
1
1
2
22
2
22
3 w
w
4.symm336
34336336
30F
4.symm612
264l612612
EI
ll
lllll
ll
llll
l
l (3.2-32)
Elastisch gebetteter BalkenWenn der Balken auf einer elastischen Unterlage liegt, für die wir die Gültigkeit derWinklerschen Hypothese annehmen können,19 erhalten wir den Einfluß der Bettungs-reaktion, wenn wir zusätzlich in Gleichung (3.2-1) den Lastterm β+= qqq mit
ßw(x)(x)q −=β (3.2-33)
ergänzen. Damit ergibt sich im Elastischen Potential der folgende Zusatzterm 18 Beispiel: Berechnung der kritischen Knickkraft für das im Bild skizzierte System.
Die Randbedingungen lauten: 0,0w,0,0w 1133 =ϕ==ϕ= . Mit die-sen Randbedingungen ergibt sich für eine Diskretisierung mit zwei finitenElemente unter Benutzung der Gleichung (3.2-32) das folgende Eigenwert-problem:
0w
4l43l3333636
42462662612212
2
22222 =
ϕ
++−+−+
λ−
⋅+⋅+−⋅+−⋅+ ∗
llll
llllll
mit EI30
F 2lλ=λ∗ und 1F = ergibt sich daraus 0
1611
22
=+λ−λ ∗∗ ; die Lö-
sungen lauten: 440983,01 =λ∗ und 559017,12 =λ∗ .Der kleinste Einwert ergibt
die kritische Kraft 2ges
2121EI
918,52EI
229,13EI30
FFlll
==λ=λ= ∗ .
19 Die Hypothese von WINKLER besagt, daß die Bettungsreaktion proportional mit der Einsenkung in das Bet-tungsmaterial zunimmt, d.h. es gilt ßw(x)q(x) −= , wenn β die Proportionalitätskonstante bezeichnet, die ausExperimenten bestimmt werden muß.
3.2. Das 2-Knoten Balkenelement (Bernoulli-Balken) 79
( ) .Mindxw21
dxwEI21
rmBettungste
l
0
2l
0
2 →−β+′′=π ∫∫ K43421
K (3.2-34)
Setzt man in (3.2-34) den Verschiebungsansatz ein, ergibt sich
f21
21 T TvvvKKvv −+=π β (3.2-35)
mit
( )∫β=βl
T dxNNK (3.2-36)
Die Minimierung von π bezüglich der Unbekannten v ergibt
fvKK =+ β )( . (3.2-37)
Natürlich ist die Bettungsmatrix nur für diejenigen Elemente aufzubauen, die auf ei-nem Bettungsmaterial aufliegen. An Hand der Formel (3.2-36) erkennt man, daß dieBettungsmatrix offenbar den gleichen Aufbau wie die Massenmatrix hat, wobei man inder Massenmatrix lediglich statt der Dichte die Bettungskonstante einsetzen muß.
80 3. Simplex-Elemente
3.3. Das 2-Knoten Rotationsscheibenelement
Eine Reihe weiterer Probleme lassen sich mit 1D-Modellen lösen. Dazu gehören bei-spielsweise rotationssymmetrische dünne
- Scheiben- Platten- Zylinderschalen.
Allerdings muß stets vorausgesetzt werden, daß sowohl die Geometrie, als auch dieBelastungen und Randbedingungen rotationssymmetrisch sind. Anderenfalls kann un-ter bestimmten Voraussetzungen mittels einer Fourier-Reihenentwicklung eine Zu-rückführung auf die Lösung einer Anzahl quasi-eindimensionaler Probleme (semi-analytische Methode) erfolgen.20
Nachfolgend betrachten wir zunächst rotationssymmetrische dünne Scheiben (sieheBild 3.3-1), die dadurch charakterisiert sind, daß sowohl die Geometrie als auch dieLagerungen und Belastungen und Materialeigenschaften rotationssymmetrisch in Um-fangsrichtung und konstant in Dickenrichtung sind, wobei angenommen wird, daß dieDicke klein ist gegenüber dem Außenradius.21 Die Mittelfläche muß natürlich ebensein. Diese Voraussetzungen führen dazu, daß die Spannungen und Verformungen un-abhängig von der Umfangsrichtung sind und nur vom Radius r abhängen. DerartigeTragwerke können damit als quasi eindimensionale Modelle behandelt werden, derenVerformungs- und Spannungszustand durch die Verschiebung u(r) vollständig be-schrieben ist.
Bild 3.3-1 Rotationssymmetrische dünne Scheibe
Der Zusammenhang zwischen der Radialverschiebung u(r) und den Dehnungen ergibtsich aus
)r(u
r1drd
)r(u
r
udr
dur
==
=
εε
=ϕ
De (3.3-1)
20 Vergleiche dazu auch die Vorlesung Flächentragwerke.21 Ein typischer Anwendungsfall dafür sind beispielweise Verdichterlaufräder in Turbinen.
r
r
r, u(r)
z
b(r)
3.3. Das 2-Knoten Rotationsscheibenelement 81
mit
=
r1drd
D (3.3-2)
Der ebenfalls zweidimensionale Spannungszustand berechnet sich über das HookescheGesetz aus den Dehnungen. Für isotropes Material gilt
εε
ν
ν
ν−=
σσ
==ϕϕ
r2
r
11
1E
Ees (3.3-3)
Durch Einsetzen von (3.3-1) und (3.3-3) wird das Elastische Potential (2.3-1) durchdie Verschiebungen u(r) ausgedrückt. Man erkennt, daß in Gleichung (2.3-1) nur ersteAbleitungen der Radialverschiebungen eingehen, so daß die finiten Elemente nur C0-
stetig sein müssen. Das leistet bereits ein linearer Ansatz für u(r), für den wir analogzum 2-Knoten Stabelement erhalten
( ) ( ) ( ) ( )[ ]
ξξ=ξ=ξ
2
121 u
uNNvNu (3.3-4)
mit
( )ξ−= 121
N1 , ( )ξ+= 121
N2 (3.3-5)
Für den Zusammenhang zwischen r und der lokalen Koordinate ξ gilt
( ) 0rr21
r +ξ∆=ξ (3.3-6)
( ) 0rr2
rr
2r
∆−
∆=ξ (3.3-7)
Damit ergibt sich für D
+ξ∆
ξ∆=
mr2r2dd
r2
D (3.3-8)
r, u(r)21
r=r1 r=r2
= -1 = 1
0
rz
Bild 3.3-2. Finites Kreisringelement mit zwei Knotenkreisen
82 3. Simplex-Elemente
Anwenden von D auf die Matrix der Formfunktionen N liefert
+ξ∆ξ+
+ξ∆ξ−
∆∆−
==
00 r2r1
r2r1
r1
r1
DNB (3.3-9)
Jetzt kann entsprechend Gleichung (2.4-5) die Steifigkeitsmatrix berechnet werden.Für die Berechnung der Volumenintegrals ist zu beachten, daß wir hier über den Radi-us, den Kreisumfang und die Dicke integrieren müssen. Für das Volumenintegral er-halten wir demzufolge
( ) ( ) ξ∆
+ξ∆ξπ=π= ∫∫∫
−
d2r
rr21
b2drrrb2dV1
10
r
rV
2
1
LLL (3.3-10)
Für K ergibt sich aus (2.4-5) mit (3.3-10) und einer konstante Elementdicke b
ξ+ξ∆ξξ∆π= ∫−
d)rr21
)(()(br1
1m
T EBBK (3.3-11)
( )
( )
∆ν++
−∆ν+−
∆ν−π
=2
1
221
1
221
2
1
222
22r1
rr
lnr.symm
rr
lnrrr1rr
lnr
r)1(bE2
K (3.3-12)
Für die Knotenkräfte infolge einer Rotation um die z-Achse mit der Winkelgeschwin-digkeit ω ergibt sich22
( )ξ∆+ξ∆ξπωρ=ωρ= ∫∫
−
dr21
)rr21
()(2bdVr1
1
2m
T2
V
2Tp GGf (3.3-13)
+−
+−
∆ωρπ
=4
1321
42
422
31
412
pr
121
rr31
r41
r121
rr31
r41
rb2
f (3.3-14)
22 Der Volumenlastvektor ergibt sich in diesem Fall zu 2rp ωρ= .
3.4. Das 2-Knoten Rotationplattenelement 83
3.4. Das 2-Knoten RotationsplattenelementVon rotationssymmetrischen dünnen Platten sprechen wir, wenn im Unterschied zumAbschnitt 3.3 die Belastung ausschließlich senkrecht zur Plattenmittelfläche erfolgt. Esgelten dann die Bild 3.4-1 angegebenen Spannungsdefinitionen.
Bild 3.4-1 Rotationssymmetrische dünne Platte
Nach der Kirchhoffschen Plattentheorie23 ergeben sich die Dehnungen in radialer undtangentialer Richtung zu
)r(wz
drdw
r1dr
wd
z2
2
r De −=
−=
εε
=ϕ
(3.4-1)
mit
=
drd
r1drd
2
2
D (3.4-2)
Die Spannungen lassen sich mit der Hookeschen Matrix durch die Dehnungen aus-drücken. Für isotropes Material gilt
HeEes2
r2
r
1E
11
1E
ν−=
εε
ν
ν
ν−==
σσ
=ϕϕ
(3.4-3)
mit
ν
ν=
11
H (3.4-4)
Die Matrix H wurde eingeführt, um eine Analogie zu der in der Plattentheorie übli-chen Schreibweisen zu erzielen. Damit kann der Ausdruck (2.3-1) für das ElastischePotential benutz werden. 23 Die wesentlichste Annahme der Kirchhoffsche Plattentheorie ist, daß eine Normale zur Plattenmittelflächeauch nach der Verformung normal auf der Mittelfläche steht. Schubspannungen τrz infolge von Querkräften ha-ben damit keinen Einfluß auf die Verformungen, was (analog zum Bernoulli-Balken) nur bei dünnen Kreisplat-ten, mit einer Dicke klein im Verhältnis zum Radius, in guter Näherung erfüllt ist. Kleine Verformungen und li-neares Materialverhalten werden hier generell vorausgesetzt. Die Berücksichtigung von Querkraftschubeinflüs-sen erfolgt im Rahmen der Mindlinschen Plattentheorie.
r
z, w(r)
b(r)
rz
rz
84 3. Simplex-Elemente
Wir wollen annehmen, daß wir für die Approximation der Durchbiegung w(r) einenAnsatz der Form
Nv=)r(w (3.4-5)
gemacht haben. Die Anwendung von D auf die Matrix der Formfunktionen liefert wieüblich
DNB = (3.4-6)
Die Berechnung der Steifigkeitsmatrix nach Gleichung (2.4-5) ergibt dann
∫ ∫ ∫∫π
+
−
ϕν−
=ν−
=2
0
2b
2b
r
r
T2
2
V
T2
22
1
dzdrdr1
EzdV
1E
z HBBHBBK (3.4-7)
Da B nur von r nicht aber von z und ϕ abhängt, kann die Integration über z und ϕ aus-geführt werden. Wenn auch hier die Elementdicke b und der E-Modul als konstant an-genommen wird, erhält man
∫
π
ν−=
2
1
r
r
T3
2 drr22b
31
21
EHBBK
∫+
−
ξ∆+ξ∆πν−
=1
102
1T2
3
dr21
)rr(2)1(12
EbHBBK (3.4-8)
Wenn man die Plattensteifigkeit
)1(12Eb
k 2
3
ν−= (3.4-9)
einführt, erhält man schließlich
∫+
−
ξ+ξ∆∆π=1
102
1T d)rr(rk HBBK (3.4-10)
In B gehen zweite Ableitungen der Durchbiegung nach der r-Koordinate ein, so daßdie Ansatzfunktionen eines finiten Elementes mindestens C1 stetig sein müssen. Dasbedeutet, daß der Näherungsansatz für die Durchbiegung w als Unbekannte minde-stens die Durchbiegung und die Biegewinkel an den beiden Knotenkreisen 1 und 2 desfiniten Elementes aufweisen muß, d.h
ϕ
ϕ=
2
2
1
1
w
w
v . (3.4-11)
3.4. Das 2-Knoten Rotationplattenelement 85
Der Verschiebungsansatz kann also analog zum Balkenelement gewählt werden (esgelten die in Bild 3.2-1 eingeführten Bezeichnungen), und man erhält in der lokalenKoordinate ξ für den Ansatz
( ) [ ]
ϕ
ϕξξξξ=ξ=ξ
2
2
1
1
4321 w
w
)(N)(N)(N)(N)(w vN (3.4-12)
mit
( ) ( )ξ+ξ−= 2141
N 21 (3.4-13a)
( )ξ−ξ−∆= 1)1(r81
N 22 (3.4-13b)
( ) ( )ξ−ξ+= 2141
N 23 (3.4-13c)
)1)(1(r81
N 24 ξ+ξ−∆−= (3.4-13d)
Die Matrix B hat die Form
)(
dd
r2
rr1
dd
r4
)(
drd
r1drd
021
2
2
22
2
ξ
ξ∆+ξ∆
ξ∆=ξ
= NNB (3.4-14)
Einsetzen der Ansatzfunktionen in (3.4-14) ergibt
( ) ( )
+ξ∆ξ−ξ−
−+ξ∆∆
ξ−+ξ∆
ξ+ξ−−+ξ∆∆
ξ+−
ξ+∆
ξ∆
−ξ+−∆
ξ∆=
)rr(4321
)rr(r233
)rr(4321
)rr(r233
31r2
1r6
31r2
1r6
021
2
021
2
021
2
021
2
22B (3.4-15)
Damit läßt sich die Steifigkeitsmatrix nach Gleichung (3.4-10) berechnen.
86 3. Simplex-Elemente
Für die Berechnung ist eine numerische Integration zweckmäßig, z.B. mit einer Gauß-schen Integrationsformel mit N Stützstellen.24
Die äquivalenten Knotenkräfte infolge von Flächenlasten q(ξ) berechnen sich nachGleichung (2.4-9) aus
∫∫−
ξ+ξ∆ξξ∆π==1
102
1T
O
Tq d)rr)((q)(rdOq NNf , (3.4-16)
wobei hier die Integration über die gesamte Oberfläche eines durch die beiden Kno-tenkreise 1 und 2 begrenzten Elementes erfolgen muß.
24 Die Gaußschen Quadraturformeln integrieren Funktionen entsprechend
∑∫=−
ξ=ξξ=N
1I
)N(I
)N(I
1
1
)(fwd)(fI ,
wobei die Gewichte und Stützstellen tabelliert vorliegen (siehe dazu auch die Vorlesung Mathematischen undNumerische Methoden der Mechanik - Teil II). Für die Stützstellen von N=1...4 gilt
N ( )NIw ( )N
Iξ1 2 02 1 ±0,57735026923 0,555...5
0,888...8±0,77459666920
4 0,3478548451 0,6521541549
±0,8611363116±0,3399810435
Für die Integration der Plattensteifigkeitsmatrix ist eine 4-Punkte Gauß-Integration zweckmäßig.
3.5. Das 3-KnotenDreieckselement für Scheibenberechnungen 87
3.5. Das 3-Knoten Dreieckselement für Scheibenberechnungen
In diesem Abschnitten behandeln wir ein ebenes Dreieckselement mit 3 Knoten für dieBerechnung von Scheibenproblemen. Das Modell einer Scheibe setzt eine ebene
Mittelfläche, eine geringe Dicke, die kleingegenüber den Flächenabmessungen ist,und Belastungen ausschließlich in derMittelfläche voraus. Wir können dannannehmen, daß die Spannungsverteilungüber die Dicke näherungsweise konstantist. Der Spannungsvektor hat dann nurdrei Komponenten, die beidenNormalspannungen σx, σy und dieSchubspannung τxy.
25 Analog gibt es zweiDehnungen εx, εy und die Gleitung γxy.
Der Zusammenhang zwischen denVerzerrungen und den Verschiebungenu(x,y) und v(x,y) wird im Falle kleinerVerschiebungen durch Gleichung (2.1-3)ausgedrückt. Für den Zusammenhang zwi-schen Spannungen und Verzerrungen gilthier das Hookesche Materialgesetz entspre-chend Gleichung (2.1-11) mit (2.1-12) fürden isotropen ebenen Spannungszustand.26
Die Verzerrungen ergeben sich aus den er-sten Ableitungen der Verschiebungen nachden (x,y)-Koordinaten. In das ElastischePotential gehen daher nur die ersten Ablei-tungen der Verschiebungen ein, so daß es of-fenbar recht, wenn auf dem gesamten Randeines zweidimensionalen finiten Elementesdie Verschiebungen zwischen benachbartenElementen stetig sind. Das einfachste finite Element ist das im Bild 3.5-1 abgebildeteDreieckselement mit 3 Knoten. Wenn es uns gelingt, einen Verschiebungsansatz ab-zuleiten, der durch die unbekannten Verschiebungen an den Knoten des Dreiecks aus-gedrückt wird, ist diese Stetigkeitsbedingung stets erfüllt wie das Bild 3.5-2 zeigt.
25 Wir verwenden hier zunächst wegen der besseren Übersichtlichkeit die Bezeichnungen x, y für die Koordina-tenrichtungen und u, v für die Verschiebungen in Richtung dieser Koordinaten.26 Manchmal ist für eine Berechnung auch der ebene Formänderungszustand erforderlich. Dazu braucht man nur
in den Gleichungen für den ebenen Spannungszustand 21
EdurchersetzenE
ν− und
ν−νν
1durchersetzen . Die
noch fehlende Spannungskomponente zσ läßt sich mit 0z =ε aus dem Hookeschen Gesetz berechnen. Man er-
hält )( yxz σ+σν=σ .
x3
y, v
12
3
u2
v2
u3
v3
u1
v1
x1 x2x, u
y2
y3
y1
Bild 3.5-1 Dreieckselement mit 3 Knotenu
u1 u2
u3
12
3
x
Bild 3.5-2 Verlauf einer linearenFunktion im Dreieck, die durch diedrei Funktionswerte an den Knoten-punkten aufgespannt wird
88 3. Simplex-Elemente
Die sechs Knotenverschiebungen fassen wir in dem Vektor v folgendermaßen zusam-men:
[ ]332211T vuvuvu=v (3.5-1)
Im Dreieckselement wird ein linearer Ansatz in der Form
yaxaa)y,x(u 321 ++= (3.5-2a)yaxaa)y,x(v 654 ++= (3.5-2b)
eingeführt. In Matrizenschreibweise lautet dieser Ansatz
==
=
6
2
1
a
aa
yx1000000yx1
vu
MPau (3.5-3)
Auswerten des Ansatzes an den Knoten liefert die A-Matrix. Die Inversion von A undbilden von:
1N −= PA (3.5-4)
ergibt
==
3
3
2
2
1
1
321
321
vuvuvu
N0N0N00N0N0N
Nvu (3.5-5)
mit:
( ) ( ) ( )[ ]A21
yxxxyyyxyxN 233223321 −+−+−= (3.5-6a)
( ) ( ) ( )[ ]A21
yxxxyyyxyxN 311331132 −+−+−= (3.5-6b)
( ) ( ) ( )[ ]A21
yxxxyyyxyxN 122112213 −+−+−= (3.5-6c)
und( ) ( ) ( )122131132332 yxyxyxyxyxyxA2 −+−+−= (3.5-7)
A ist der Flächeninhalt des Dreiecks. Die Matrix B=DN läßt sich damit folgenderma-ßen schreiben:
( ) ( ) ( )
( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( )
−−−−−−−−−
−−−=
21121313313223
123123
211332
yyxxyyxxyyxx
xx0xx0xx0
0yy0yy0yy
A21
B (3.5-8)
3.5. Das 3-KnotenDreieckselement für Scheibenberechnungen 89
Da B konstant ist, läßt sich K problemlos berechnen. Wenn auch die Elementdickeund die Hooksche Matrix elementweise konstant sind, ergibt sich
EBBEBBEBBK T
A
T
V
T bAdAbdV === ∫∫ (3.5-9)
Bei einer Linienlast, die senkrecht auf einer Elementkante (siehe Bild 3.5-3) steht, er-gibt sich für die beiden Kräfte F1 und F2
27
( )211 qq261
F += l (3.5-10a)
( )212 q2q61
F += l (3.5-10b)
Die Kräfte F1 und F2 stehen ebenfalls senk-recht auf der Seitenkante. Vor der Einspeiche-rung in das globale Gleichungssystem müssendiese Kräfte noch auf das globale Koordina-tensystem transformiert werden. Falls eineLinienlast tangential zur Elementkante wirkt,ergeben sich analog Knotenkräfte in Richtungdie Elementkante, die ebenfalls noch auf dasglobale Koordinatensystem transformiertwerden müssen.
Natürliche Koordinaten im Dreieck (Dreieckskoordinaten)Bei eindimensionalen Problemen haben wir durch die Einführung der dimensionslosennatürlichen Koordinate ξ ∈[-1,1] die Ableitung von finiten Elementen vereinfacht. BeiDreieckselementen führt die Benutzung von Dreieckskoordinaten (natürliche Koordi-naten im Dreieck) zu ähnlichen Vorteilen (siehe Bild 3.5-4). Die Dreieckskoordinatenhaben ihren Ursprung jeweils auf den Seitenkanten und nehmen am gegenüberliegen-den Knoten den Wert 1 an, d.h. es gilt ζL∈[0,1], L=1,2,3. Ein Punkt P im Dreieck hatdamit die folgenden Koordinaten (siehe Bild 3.5-4).
AA1
P1 =ς ,AA2
P2 =ς ,A
A3P3 =ς (3.5-11)
27 Für die Berechnung reicht es aus, nur die Formfunktionen entlang der Seitenkante zu betrachten. Diese ist of-fensichtlich linear und kann als 22
112
1 us1us1su )()()( ++−= geschrieben werden, wobei hier s die natürliche
Koordinate längs der Seitenkante ist, und u1, u2 sind die Verschiebungen der Knoten normal zur Seitenkante (po-sitiv in Richtung auf das Elementinnere). Dann berechnen sich die Knotenkräfte aus
+
+=++−
+−
=
+−
=
= ∫∫
−− 21
1122
112
11
1 2121
2
1
1 2121
2
1
q2qqq2
6ds}q)s1(q)s1({
)s1()s1(
2ds)s(q
)s1()s1(
FF lllF .
s=1
1
2
F1
F2
q2
q1
l
s
s=-1
Bild 3.5-3 Linienlast auf einer Ele- mentkante
90 3. Simplex-Elemente
Die drei Dreieckskoordinaten sind nicht voneinander unabhängig. Es gilt
1A
AAA 321P3P2P1 =
++=ς+ς+ς (3.5-12)
Zwischen den Koordinaten x, y und ζ1, ζ2, ζ3 besteht ein linearer Zusammenhang.
332211 xxxx ς+ς+ς= (3.5-13a)
332211 yyyy ς+ς+ς= (3.5-13b) 321?1 ς+ς+= (3.5-13c)
Das Gleichungssystem (3.5-13) kann nach ζ1, ζ2, ζ3 umgestellt werden, und man erhält
( )( )( )
−−−−−−−−−
=
ςςς
1yx
yxyx)xx()yy(yxyx)xx()yy(yxyx)xx()yy(
A21
12211221
31133113
23322332
3
2
1
(3.5-14)
Für die Ableitungen ergibt sich: 28
ς∂∂ς∂∂
−−−−
=
ς∂∂ς∂∂
∂ς∂
∂ς∂
∂ς∂
∂ς∂
=
∂∂
∂∂
2
1
3123
)1332
2
1
21
21
)xx()xx(yy()yy(
A21
yy
xx
y
x (3.5-15)
28 Wobei hier vorausgesetzt wurde, daß die abzuleitenden Funktionen im lokalen Koordinatensystem nur nochvon 1ς und 2ς abhängig sind, d.h. die 3ς Koordinate schon mittels Gleichung (3.5-13c) eliminiert wurde.
y, v
12
3
x, u
P
A
AA1
3
2P
1P
3P
2
Bild 3.5-4 Natürliche Koordinaten im Dreieck
3.5. Das 3-KnotenDreieckselement für Scheibenberechnungen 91
Mit den Dreieckskoordinaten lassen sich sehr systematisch Dreieckselemente entwi k-keln. So lauten beispielsweise die Ansatzfunktionen für das 3-Knoten Dreiecks-element (siehe Bild 3.5-4):29
11N ς= (3.5-16a)
22N ς= (3.5-16b)
2133 1N ς−ς−=ς= (3.5-16c)
Mit Hilfe der Dreieckskoordinaten lassen sich Integrationsformeln für die Integrationüber Dreiecksflächen herleiten. So gilt beispielsweise für die Integration in den natür-lichen Koordinaten30
( )!2pnm!p!n!m
A2dAIA
p3
n2
m1pnm +++
=ςςς= ∫ (3.5-17)
Für die Integration über (x, y) muß das Integral zunächst auf die natürlichen Dreiecks-koordinaten transformiert werden, um die Formel (3.5-17) anwenden zu können. Fürdie Lösung das Integrals
∫=A
srrs dAyxI (3.5-18)
ergibt sich
∑∑= =
−−
−−+−−+
++=
r
0i
s
0j
j2
js1
i2
ir1rs yyxx
)!js()!ir(!j!i)!ji()!jisr(
A2)!2sr(
!s!rI , (3.5-19)
wenn man ohne Einschränkung der Allgemeinheit annimmt, daß das lokale Koordina-tensystem )y,x( genau im Knoten 3 des Dreieckselementes liegt, d.h. für die Koordi-naten des Knoten 3 gilt 0yx 33 == .31
29 Für ein Dreieckselement mit 6 Knotenpunkten (quadratischer Ansatz) ergibt sich
y, v
12
3
4
56
( )( )( )
316
325
214
333
222
111
4N4N4N
12N12N
12N
ςς=ςς=ςς=
−ςς=−ςς=
−ςς=
30 Die Herleitung ergibt sich mittels partieller Integration.31 Für eine beliebige Lage des Koordinatensystem gilt ( ) ( )∫∫ ++== dAyyxxdAyxI n
3m
3nm
mn . Dieses
Integral kann durch Auflösung des Integranden mittels der binomischen Formel berechnet werden, wobei dabeiIntegrale der Form (3.5-19) zu berechnen sind.
92 3. Simplex-Elemente
3.6. Das 4-Knoten Rechteckelement für Scheibenberechnungen
Nachfolgend behandeln wir ein Rechteckelement mit 4 Knotenpunkten (siehe Bild3.6-1) und nehmen zunächst an, daß die Seitenkanten parallel zu den Koordinatenach-
sen orientiert sind. Dieses Element isteine Erweiterung des eindimensionalenElementes (siehe Abschnitt 3.1), sodaß wir die Ansatzfunktionen für die-ses Element offensichtlich durch mul-tiplikative Verknüpfung der eindimen-sionalen Ansatzfunktionen aus Glei-chung (3.1-4) ableiten können. In demlokalen (ξ, η) Koordinatensystem istdas Element ein Quadrat der Kanten-länge 2. Da wir jetzt vier Knoten ha-ben, können wir die Verschiebungen injeder Koordinatenrichtung durch vierPolynomanteile approximieren, z.B.durch (1, x, y, xy). Beim Dreiecksele-
ment mit drei Knoten bestand der Ansatz aus den Polynomanteilen (1, x, y), durch diegerade ein vollständiges lineares Polynom beschrieben wird. Wegen des zusätzlichenKnotens haben wir hier den zusätzlichen Anteil (xy) im Ansatz, der die Genauigkeitder Viereckselement gegenüber den Dreieckselementen deutlich verbessert.32 Wirwollen nachfolgend die Ansatzfunktionen unter Nutzung der eindimensionalen An-satzfunktionen (3.1-4) direkt angeben.
( ) ( ) ( )( )η−ξ−=
η−
ξ−=ηξ 11
41
121
121
),(N1 (3.6-1a)
( ) ( ) ( )( )η−ξ+=
η−
ξ+=ηξ 11
41
121
121
),(N2 (3.6-1a)
( ) ( ) ( )( )η+ξ+=
η+
ξ+=ηξ 11
41
121
121
),(N3 (3.6-1a)
( ) ( ) ( )( )η+ξ−=
η+
ξ−=ηξ 11
41
121
121
),(N4 (3.6-1d)
Die Gleichungen (3.6-1) lassen sich auch in kompakter Form als
( )( )ηη+ξξ+=ηξ LLL 1141
),(N , L=1,...4 (3.6-2)
schreiben, wobei (ξL, ηL) die lokalen Koordinaten des Knotens L sind (ξL, ηL∈[-1,1]).Die Ansatzfunktionen (3.6-1) enthalten genau die Polynomterme (1, ξ, η, ξη) wie mandurch Ausmultiplizieren zeigen kann. Die Summe der Formfunktionen ergibt 32 Das ist übrigens der Grund für die Bevorzugung von Viereckselementen gegenüber Dreieckselementen undHexaederelementen gegenüber Teraederelementen, wo der Unterschied infolge des deutlich höheren Grades derPolynomansätze noch stärker ist.
y, v
1 2u2
v2
3u3
v3
u1
v1
x, u
4 u4
v4
b 0
a
Bild 3.6-1 Rechteckelement mit 4 Knoten
3.6. Das 4-Knoten Rechteckelement für Scheibenberechnungen 93
1),(N4
1LL =ηξ∑
=, (3.6-3)
so daß auch die Bedingung (2.3-18) erfüllt ist, die eine notwendige Bedingung dafürist, daß das Element bei Starrkörperbewegungen korrekt reagiert, d.h. zum Beispielkeine Spannungen liefert. Damit ergibt sich der Verschiebungsansatz des Elementes zu
=ηξ=
=
4
4
3
3
2
2
1
1
4321
4321
vuvuvuvu
N0N0N0N00N0N0N0N
),(vu
vNu (3.6-4)
Zu dem gleichen Ansatz gelangt man übrigens auch, wenn man von dem Polynoman-satz
ξηηξ
ξηηξ=ηξ=
ηξ
ηξ=ηξ
8
1
a
a
10000
00001),(
),(v
),(u),( MaPu (3.6-5)
ausgeht, diesen Ansatz an den 4 Knotenpunkten auswertet, die so entstehende MatrixA invertiert und N(ξ,η)=A-1P(ξ,η) berechnet (siehe die allgemeine Darstellung derAbleitung des Verschiebungsansatzes in Abschnitt 2.3). Aus dem Zusammenhangzwischen den Verzerrungen und Verschiebungen Gleichung (2.1-3) erhält man dieDifferentiationsmatrix D, deren Anwendung auf die Matrix der Formfunktionen nachGleichung (2.4-2) die Matrix B ergibt. Um B berechnen zu können, müssen wir dieAbleitungen nach x und y noch durch die Ableitungen nach den lokalen Koordinaten ξund η ersetzen. Der Zusammenhang zwischen den Koordinaten folgt aus Bild 3.6-1 zu
0xa21
x +ξ= (3.6-6a)
0yb21
y +η= (3.6-6b)
Daraus folgt für die Ableitungen
ξ∂∂
=ξ∂∂
∂ξ∂
=∂∂
a2
xx (3.6-7a)
η∂∂
=η∂∂
∂η∂
=∂∂
b2
yy (3.6-7b)
Die Gleichungen (3.6-7) lassen sich in Matrizenschreibweise als
94 3. Simplex-Elemente
η∂∂ξ∂
∂
=
∂∂∂∂
R
y
x (3.6-8)
schreiben, wobei sich die Matrix R folgendermaßen berechnet:
=
∂η∂
∂ξ∂
∂η∂
∂ξ∂
=
b2
0
0a2
yy
xxR . (3.6-9)
Die Matrix R ist die Inverse der bekannten Jacobischen Matrix J.33 Die besonders ein-fache Form ergibt sich hier dadurch, daß die Elementkanten und die Achsen des loka-len Koordinatensystems parallel zum globalen Koordinatensystem liegen (im Ab-schnitt 5. wird ein allgemeinerer Fall betrachtet). Wir können jetzt die Matrix B ange-ben, die wir für die Berechnung der Steifigkeitsmatrix benötigen. Wir wenden dazu dieDifferentiationsmatrix D
ξ∂∂
η∂∂
η∂∂
ξ∂∂
=
∂∂
∂∂
∂∂
∂∂
=
a2
b2
b2
0
0a2
xy
y0
0x
D (3.6-10)
auf die Matrix der Ansatzfunktionen N in Gleichung (3.6-4) an und erhalten
[ ]4321),(),( BBBBDNB =ηξ=ηξ (3.6-11)
mit
=
xLyL
yL
xL
L
BBB0
0BB (3.6-12)
und
)1(a2
1)1(
41
a2),(N
a2
B LLLLL
xL ηη+ξ=ηη+ξ=ξ∂
ηξ∂= (3.6-13a)
)1(b2
1)1(
41
b2),(N
b2
B LLLLL
yL ξξ+η=ξξ+η=η∂
ηξ∂= (3.6-13b)
33 Die JACOBIsche Matrix J=R-1 lautet hier
=
=
∂ξ∂
∂ξ∂−
∂∂
−∂∂
=b00a
21
a20
0b2
ab41
xx
x?
y?
det1
RJ ,
und damit ergibt sich dA=dxdy=detJdξdη=¼abdξdη.
3.6. Das 4-Knoten Rechteckelement für Scheibenberechnungen 95
Damit kann die Steifigkeitsmatrix nach Gleichung (2.4-5) berechnet werden. Die Inte-gration über die Dickenrichtung liefert die Dicke h(ξ,η) als Funktion der lokalen Ko-ordinaten. Wenn die Dicke als Funktion bekannt ist, kann sie durch die lokalen Koor-dinaten ausgedrückt und in die Integration einbezogen werden. Häufig wird eine ver-änderliche Dicke durch die Dickenwerte hL an den Knotenpunkten L des Elementnet-zes beschrieben. Dann kann man die Verschiebungsansatzfunktionen auch benutzen,um die Dickenverteilung im Element zu approximieren, so daß man
∑=
ηξ=ηξ4
1LLL h),(N),(h (3.6-15)
erhält. Die Steifigkeitsmatrix hat folgende Form
ηξ
ηξ=
ηξηξ=
∫ ∫
∫ ∫
− −
− −
dd
.symm
),(hab41
dd2b
2a),(h
1
1
4T4
4T33
T3
4T23
T22
T2
4T13
T12
T11
T1
1
1
1
1
1
1
T
EBBEBBEBBEBBEBBEBBEBBEBBEBBEBB
BEBK
(3.6-16)
Jede der 2×2 Submatrizen KIJ gehört zu einer Knotenkombination (I,J) mit I,J=1,...,4.Damit kann man die Matrix K im Rahmen eines Rechenprogramms systematisch ausdiesen Submatrizen aufbauen. Für den Fall einer orthotropen Hookschen Matrix E mit
=
33
2212
1211
E000EE0EE
E (3.6-17)
ergibt sich für eine solche Submatrix
ηξ
++++
ηξ=
ηξηξ=
∫ ∫
∫ ∫
− −
− −
ddBEBBEBBEBBEBBEBBEBBEBBEB
),(hab41
dd),(hab41
xJ33xIyJ22yIyJ33xIxJ12yI
xJ33yIyJ12xIyJ33yIxJ11xI1
1
1
1
ITI
1
1
1
1IJ EBBK
(3.6-18)
In diesem Fall kann die Integration noch exakt ausgeführt werden. Wenn man sich dieFunktionen in Gleichung (3.6-13) ansieht, erkennt man, daß in (3.6-18) nur vier Typenvon Integralen auftreten, wenn man von konstanten Faktoren absieht. Das sind
∫−
=ξ±1
1
2
38
)1( , ∫−
=ξ−1
1
2
34
)1( , 2d1
1
=ξ∫−
, ∫−
=ξ±1
1
2)1( (3.6-19)
Damit kann problemlos die Steifigkeitsmatrix berechnet34 und programmiert werden.
34 Für das erste Element der Steifigkeitsmatrix ergibt sich beispielsweise )E1E(b31K 3311
11IJ α
+α= , mit ab=α .
96 3. Simplex-Elemente
3.7. Das 4-Knoten Tetraederelement
In den folgenden beiden Abschnitten soll noch kurz auf dreidimensionale finite Ele-mente eingegangen werden. Das indiesem Abschnitt behandelte 4-Knoten Tetraederelement ist eineErweiterung des Dreieckselemen-tes mit drei Knotenpunkten (sieheAbschnitt 3.5.). Das Bild 3.7-1zeigt die Geometrie und Knoten-anordnung für dieses Element. DieBerechnung kann in Analogie zumDreieckselement erfolgen. Es istoffensichtlich, daß wir in jeder derdrei Koordinatenrichtungen einenlinearen Näherungsansatz für dieApproximation des Verschie-bungszustandes einführen müssen.
Der Ansatz kann in Analogie zur Gleichung (3.5-2) aufgeschrieben werden und enthälthier die Polynomterme (1, x, y, z). Die Auswertung dieses Ansatzes an den vier Kno-tenpunkten führt zur Matrix A, deren Inversion schließlich analog zur Gleichung (3.5-4) zur Matrix der Formfunktionen G führt. Die Ableitung der Elementsteifigkeitsma-trix und der Lastvektoren kann schließlich auf dem üblichen Wege erfolgen. Wir wo l-len hier gleich den eleganteren Lösungsweg unter Nutzung der lokalen Tetraederkoor-dinaten beschreiten. Die Tetraederkoordinaten ζL ∈[0,1], L=1,...,4, sind analog zu denDreiekcskoordinaten definiert. Sie nehmen jeweils am Knoten L den Wert Eins an,und sind auf der dem Knoten gegenüberliegenden Seite Null. Natürlich sind nur dreidieser Koordinaten unabhängig voneinander, und die vierte Koordinate ζ4 kann durchdie übrigen drei Koordinaten ausgedrückt werden, wobei dazu die Gleichung
43211 ς+ς+ς+ς= (3.7-1)
genutzt wird.35 Für den Zusammenhang zwischen den Tetraederkoordinaten und denkartesischen Koordinaten gilt analog zu Gleichung (3.5-13)
ςςςς
=
ςςςς
=
4
3
2
1
4
3
2
1
4321
4321
4321
1111zzzzyyyyxxxx
1zyx
H (3.7-2)
In der Koeffizientenmatrix stehen die Koordinaten der vier Knotenpunkte des Te-traeders. Die Inversion von (3.7-2) ergibt 36
35 Diese Gleichung läßt sich analog zum Dreieck dadurch erhalten, daß man von einem Punkt P im Tetraeder dievier Teilvolumina VL mit L=1,...4 berechnet. Diese Teilvolumina ergeben sich mit den Tetraederkoordinaten zu
LL VV ς= . Die Addition liefert ∑ ∑= =
ς=ς=4
1L
4
1LLL VVV und daraus folgt ∑
=
ς=4
1LL1 .
36 Für die Determinante von H ergibt sich detH=6V, wobei V das Volumen des Tetraeders ist.
y, v
z, w
x, u
1
2
3
4
Bild 3.7-1 Tetraederelement mit 4 Knotenpunkten
3.7. Das 4-Knoten Tetraederelement 97
=
ςςςς
1zyx
4
3
2
1
J (3.7-3)
mit
== −
4441
14
1
JJ
J
L
MJHJ (3.7-4)
Damit können die Ableitungen mittels
ς∂∂ς∂∂ς∂∂
=
∂∂∂∂∂∂
3
2
1
T
z
y
x
J (3.7-5)
transformiert werden, wobei auch hier wieder vorausgesetzt wird, daß die zu differen-zierenden Funktionen nur noch von den lokalen Koordinaten ζ1, ζ2, ζ3 abhängig sind.Die Matrix J ist konstant und hängt nur von der Geometrie des Elementes ab. Für dieAnsatzfunktionen in den Tetraederkoordinaten erhält man die besonders einfacheForm37
LLN ς= , L=1,2,3,4 (3.7-6)
Wenn man die Koordinate ζ4 ersetzt, gilt
3214 1N ς−ς−ς−= (3.7-7)
Der Verschiebungsansatz lautet damit
37 Für ein Tetraederelement mit 10 Knoten (quadratischer Verschiebungsansatz) gilt:
1
2
3
4
5
6
8
10
7 9
Eckknoten (L=1,2,3,4):)12(N LLL −ςς=
Seitenmittenknoten (L=5,6,...,10):
4310
326
215
4N
4N4N
ςς=
ςς=ςς=
L
98 3. Simplex-Elemente
[ ]
ς−ς−ς−ςςς==
=
4
3
2
1
321321 )1(wvu
u
uuuu
IIIINv , (3.7-8)
mit
=
100010001
I (3.7-9)
und
=
L
L
L
L
wvu
u (3.7-10)
Damit stehen alle Grundlagen bereit, um die Steifigkeitsmatrix K (12×12) für diesesElement zu berechnen. Der Zusammenhang zwischen Verzerrungen und Verschiebun-gen wird im 3D Fall durch Gleichung (2.1-4) beschrieben. Damit ergibt sich die Diffe-rentiationsmatrix zu
ς∂∂ς∂∂ς∂∂
=
ς∂∂ς∂∂ς∂∂
+++++++++
=
∂∂
∂∂
∂∂
∂∂∂∂
∂∂
∂∂
∂∂
∂∂
=
3
2
1
3
2
1
333123211311
333223221312
323122211211
332313
322212
312111
JJJJJJJJJJJJJJJJJJ
JJJJJJ
JJJ
x0
z
yz0
0xy
z00
0y
0
00x
RD (3.7-8)
Die Ableitung der Formfunktionen nach den lokalen Tetraederkoordinaten liefert
−−
−=
ς∂∂ς∂∂ς∂∂
=100100000000
010000010000001000000001
3
2
1
NS (3.7-9)
Jetzt kann die Matrix B, die für dieses Element offensichtlich konstant ist, berechnetwerden aus
3.7. Das 4-Knoten Tetraederelement 99
RSDNB == (3.7-10)
Dann ergibt sich für die Steifigkeitsmatrix
VdV T
V
T EBBBEBK == ∫ , (3.7-11)
wobei das Volumen aus HdetV 61= berechnet werden kann, wenn man die Knoten-
koordinaten des Tetraederelementes kennt. Hier ist also der Integrand in (3.7-11) kon-stant, so daß nur mit dem Volumen zu multiplizieren ist. Bei der Verwendung höhererAnsatzfunktionen ist das nicht mehr der Fall. Für die dann erforderlichen Integrationenist folgende Formel hilfreich:38
)!3qpnm(!q!p!n!Vm6
dVJ q4
V
p3
n2
m1mnpq ++++
=ςςςς= ∫ (3.7-12)
38 Die Ableitung kann wieder unter Beachtung der veränderlichen Integrationsgrenzen und Anwendung der par-tiellen Integration abgeleitet werden.
100 3. Simplex-Elemente
3.8 Das 8-Knoten Hexaederelement39
Das 8-Knoten Hexaederelemnt ist eine Erweiterung des 4-Knoten Rechteckelementes(siehe Abschnitt 3.6). Die Ableitung der Formfunktionen und die Berechnung der Stei-
figkeitsmatrix und der Kraftvektorenerfolgt analog zu der Darstellung inAbschnitt 3.6. Das Element hat 8Knotenpunkten, so daß in jeder derdrei Koordinatenrichtung ein Ver-schiebungsansatz mit den Termen
(1, x, y, z, xy, yz, zx, xyz)verwirklicht werden kann. Nebendem vollständigen linearen Ansatz(1,x,y,z) enthält dieser Ansatz nochvier zusätzliche Terme, durch die die-ses Elemnt wesentlich genauer als dasTetraederelement mit 4 Knotenpunk-ten ist. Natürlich hätten auch andere
Terme in dem Ansatz benutzt werden können. Es sollte aber stets das Ziel verfolgtwerden, die Ansatzfunktionen so zu konstruieren, daß keine Richtung bevorzugt wird,und außerdem sollten die Ansatzfunktionen keine Polynomanteile mit unnötig hohenPotenzen enthalten, die nicht Teil eines vollständigen Polynoms sind. Die angegebe-nen Terme des Ansatzes sichern übrigens, daß die Verschiebungen in jeder der dreiKoordinatenrichtungen durch eine lineare Funktion approximiert werden. Damit kannman wie beim Rechteckelement auch hier die Ansatzfunktionen aus einer Kombinati-on der linearen Ansatzfunktionen entwickeln. In den lokalen Koordinaten erhält mandann
( ) ( ) ( )( )
( ) ( ) ( )( )
( ) ( ) ( )( )
( ) ( ) ( )( )
( ) ( ) ( )( )
( ) ( ) ( )( )
( ) ( ) ( )( )
( ) ( ) ( )( ) )1(1181
)1(21
121
121
),,(N
)1(1181)1(
211
211
21),,(N
)1(1181
)1(21
121
121
),,(N
)1(1181)1(
211
211
21),,(N
)1(1181
)1(21
121
121
),,(N
)1(1181
)1(21
121
121
),,(N
)1(1181)1(
211
211
21),,(N
)1(1181
)1(21
121
121
),,(N
8
7
6
5
4
3
2
1
ς+η+ξ−=
ς+
η+
ξ−=ςηξ
ς+η+ξ+=
ς+
η+
ξ+=ςηξ
ς+η−ξ+=
ς+
η−
ξ+=ςηξ
ς+η−ξ−=
ς+
η−
ξ−=ςηξ
ς−η+ξ−=
ς−
η+
ξ−=ςηξ
ς−η+ξ+=
ς−
η+
ξ+=ςηξ
ς−η−ξ+=
ς−
η−
ξ+=ςηξ
ς−η−ξ−=
ς−
η−
ξ−=ςηξ
(3.8-1)
Diese Ansatzfunktionen enthalten in den lokalen Koordinaten genau die oben angege-benen Polynomterme (1, ξ, η, ζ, ξη, ηζ, ζξ, ξηζ). Die Summe der Formfunktionen
39 Im englischen Sprachgebrauch als brick element bezeichnet.
1
2
345
6
7
8
x
y
z
a b
c
Bild3.8-1 Hexaederelement mit 8 Knoten
3.8. Das 8-Knoten Hexaederelement 101
(3.8-1) ergibt auch hier wieder den Wert 1, so daß eine Starrkörpertranslation zu kei-nen Spannungen führt. Die acht Ansatzfunktionen lassen sich kompakt in der Form
)1)(1)(1(81
),,(N LLLL ςς+ηη+ξξ+=ςηξ (3.8-2)
schreiben, wobei die (ξL, ηL, ζL) die lokalen Koordinaten der Knotenpunkte L=1,...,8sind mit (ξL, ηL, ζL∈[-1,1]). Analog zu Gleichung (3.6-4) läßt sich der Verschiebungs-ansatz in der Form
[ ]
=ςηξ=
=
8
8
8
1
1
1
8L1
wvu
wvu
NNN),,(wvu
MKK IIIvNu , (3.8-3)
mit
),,(N1
11
N LL ςηξ
=I (3.8-4)
schreiben. Die weiteren Ableitungen folgen analog zur Darstellung in Abschnitt 3.6und sollen hier nicht weiter vertieft werden.
Top Related