2. Grundlagen der FEM für die Elastomechanik 2bis3-8.pdf · 2. Grundlagen der FEM für die...

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.


[ ]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.


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.


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ν−ν+

ν=λ

ν+==µ


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= .


ν−ν−

ν−ν−

ν−ν−

−νν−νν

−νν

ν−ν+ν−

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


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.


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.


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.


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.


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.


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)


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.


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.


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.


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.


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:


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.


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.


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


[ ]∫∫−−

ξ−⋅

−=ξ=

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 ρ= .


[ ] ξξ+ξ−

ξ+ξ−

ρ= ∫ 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.


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.


−−

−−−−

=−

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



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

=α


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

qq

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.


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)


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ß.


( ) .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ß.


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


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


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.


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


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)


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


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


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.


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


η∂∂ξ∂

∂

=

∂∂∂∂

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=α .


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


[ ]

ς−ς−ς−ςςς==

=

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.


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.

2. Grundlagen der FEM für die Elastomechanik 2bis3-8.pdf · 2. Grundlagen der FEM für die...

Documents

Transcript of 2. Grundlagen der FEM für die Elastomechanik 2bis3-8.pdf · 2. Grundlagen der FEM für die...